[R] (no subject)

peter dalgaard pdalgd at gmail.com
Sun Jan 3 14:41:14 CET 2016


Your code could do with some revision: If I read the curly braces correectly, you simulate 1000 data sets and define the same function 1000 times, then use the function on the last simulated data set. Also, inside the function, you define gama and beta, but do not actually use them.

However, none of that would cause your symptoms. More likely something is wrong in in your function P so that it does not actually compute the negative log likelihood.

I'm not familiar with this particular distribution, and I think it would be your job and not mine to check against the theory. However, looking at the function as written, I see two divergent terms as theta[2] approaches zero: 

-n/2*log(theta[2]) goes to plus infinity

sum(-0.5/theta[1]^2*(x/theta[2])) goes to minus infinity,
if theta[1] > 0

The latter term is faster so it seems that your function has no lower limit as theta[2] approaches zero. And optim looks for a minimum.

This suggests a sign error somewhere. However, is there not a dgbs function in the gbs package, to compute the density of the distribution? If so, it would be much easier to minimize 

-sum(dgbs(x, gama, beta, log=TRUE))  

or, in case the log argument doesn't work

-sum(log(dgbs(x, gama, beta))

-pd

> On 03 Jan 2016, at 01:02 , Muhammad Kashif <mkashif at uaf.edu.pk> wrote:
> 
> Dear i optimized the gama and beta value using MLE via simmulation of birnbaum saunders distribution. if i run this code it generate very small value of beta. Could any one help me in this regard. i use gbs package to generate data.
> 
> gama=1.0
> beta=1.3
> n=25
> iterCount=1000
> for(i in 1:iterCount){
> x<-rgbs(n,gama,beta)
> P<-function(theta,x){
> n<-length(x)
> gama<-theta[1]
> beta<-theta[2]
> -n*log(5.013)+ n*(theta[1])^-2-n*log(theta[1])-n/2*log(theta[2])+ sum(log(theta[2]+x)-0.5/theta[1]^2*(x/theta[2]+theta[2]/x))}}
> P.out<-optim(theta<-c(gama,beta),ll.wd,x=x,method = "Nelder-Mead",hessian=FALSE)
> gamhat<-P.out$par[1]
> betahat<-P.out$par[2]
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com



More information about the R-help mailing list