That was great! I tried as what you said, not as many warnings as before,
and it seems that it converges.
Thank you Berend. and thank you all!
Carol
On Sat, May 8, 2010 at 5:00 PM, Berend Hasselman <b...@xs4all.nl> wrote:
>
>
> Carol Gao wrote:
> >
> > I have generated some random variable from generalised gamma
> distribution,
> > and now I need to use them to test on my optim function. Please see
> below:
> >
> > *##Need to install package VGAM first before calling rggamma function
> > library(VGAM)
> > x <- rggamma(500,scale=gamma(1.5)/gamma(3.5),d=0.5,k=1.5)
> > mlogl <- function(param){
> > n <- length(x)
> > psi <- numeric(n)
> > psi[1] <- 1.0
> > a0 <- param[1]; a1 <-param[2]; b1 <- param[3]; alpha <- param[4]; kappa
> > <-
> > param[5];
> >
> > for (i in 2:n) {
> > psi[i] <- a0 + a1*x[i-1] + b1*psi[i-1]
> > lam <- gamma(kappa)/gamma(kappa+(1/alpha))
> > ll <-
> >
> n*log(alpha/gamma(kappa))+kappa*alpha*sum(log(x))-kappa*alpha*sum(log(psi*lam))-lam^(-alpha)*sum((x/psi)^alpha)
> > }
> > return(-ll)
> > }
> > fittedparam<-optim(c(0.1,0.05,0.5,0.65,1.8),mlogl,x)*
> >
>
> The for loop in your problem is wrongly structured.
> You only need to compute lam once, but you're doing that (n-1) times.
> Same goes for the computation of ll.
> Since psi was initialised to all zeros by the numeric(n), psi*lam will
> still
> contain zero's while i < n.
>
> So move the calculation of lam and ll outside the for loop and you
> shouldn't
> get so many messages about NaN
>
> for (i in 2:n) { psi[i] <- a0 + a1*x[i-1] + b1*psi[i-1] }
>
> lam <- gamma(kappa)/gamma(kappa+(1/alpha))
> ll <-
>
> n*log(alpha/gamma(kappa))+kappa*alpha*sum(log(x))-kappa*alpha*sum(log(psi*lam))-lam^(-alpha)*sum((x/psi)^alpha)
>
> /Berend
> --
> View this message in context:
> http://r.789695.n4.nabble.com/problem-in-using-optim-tp2133938p2135929.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help@r-project.org mailing list
> 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.
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
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.
