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.
