[R] problem in using optim

[Berend Hasselman]

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.