Dear R helpers, I have problems to properly define a Likelihood function. Thanks to your help my basic model is running quite well, but I have problems to get the enhanced version (now incorporating covariates) running.
Within my likelihood function I define a variable 'alpha'. When I want to optimize the function I get the error message: 'Error in fn(ginv(par), ...) : object 'alpha' not found' I think it's actually not a problem with the optimization function (nmkb), but with the Likelihood function itself. I do not understand why 'alpha' is a missing object. 'alpha' should be part of the dataframe 'data' (as 'beta' should be too), like 'x', 'tx', ''T. But it obviously isn't. Here's a minimum example which reproduces my problem: ################################################################## library(plyr) library(dfoptim) ### Sample data ### x <- c(3, 0, 2, 5, 1, 0, 0, 1, 0, 2) tx <- c(24.57, 0.00, 26.86, 34.57, 2.14, 0.00, 0.00, 8.57, 0.00, 14.29) T <- c(33.29, 30.71, 31.29, 34.57, 36.00, 35.43, 31.14, 33.86, 35.71, 35.86) IS <- c(54.97, 13.97, 122.33, 110.84, 30.72, 14.96, 30.72, 20.74, 29.16, 83.00) data <- data.frame(x=x, tx=tx, T=T) rm(x, tx, T) ### Likelihood function ### Likelihood_cov <- function(params, x, tx, T, IS) { r <- params[1] alpha_zero <- params[2] s <- params[3] beta_zero <- params[4] gamma_1 <- params[5] gamma_2 <- params[6] data$alpha <- alpha_zero*exp(-gamma_1*IS) data$beta <- beta_zero*exp(-gamma_2*IS) f <- function(x, tx, T, alpha, beta) { g <- function(y) (y + alpha)^(-( r + x))*(y + beta)^(-(s + 1)) integrate(g, tx, T)$value } integral <- mdply(data, f) L <- exp(lgamma(r+x)-lgamma(r)+r*(log(alpha)-log(alpha+T))-x*log(alpha+T)+s*(log(beta)-log(beta+T)))+exp(lgamma(r+x)-lgamma(r)+r*log(alpha)+log(s)+s*log(beta)+log(integral$V1)) f <- -sum(log(L)) return (f) } ### ML optimization ### params <- c(0.2, 5, 0.2, 5, -0.02, -0.02) fit <- nmkb(par=params, fn=Likelihood_cov, lower=c(0.0001, 0.0001, 0.0001, 0.0001, -Inf, -Inf), upper=c(Inf, Inf, Inf, Inf, Inf, Inf), x=data$x, tx=data$tx, T=data$T, IS=IS) ################################################################## Maybe you could give me a hint were the flaw in my code is. Many thanks in advance. Carlos ----------------------------------------------------------------- Carlos Nasher Buchenstr. 12 22299 Hamburg tel: +49 (0)40 67952962 mobil: +49 (0)175 9386725 mail: carlos.nas...@gmail.com [[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.