[R] [dfoptim] 'Error in fn(ginv(par), ...) : object 'alpha' not found'

[Simon Zehnder]

Hi Carlos,

your problem is a wrong definition of your Likelihood function. You call symbols in the code (alpha, beta) which have no value assigned to. When L the long calculation in the last lines is assigned to L alpha and beta do not exist. The code below corrects it. But you have a problem with a divergent integral when calling integrate. A problem you can surely fix <as you know what your function is doing.

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_zero)-log(alpha_zero+T))-x*log(alpha_zero+T)+s*(log(beta_zero)-log(beta_zero+T)))+exp(lgamma(r+x)-lgamma(r)+r*log(alpha_zero)+log(s)+s*log(beta_zero)+log(integral$V1))
 f <- -sum(log(L))
 return (f)
}


Best

Simon


On Sep 3, 2013, at 1:28 PM, Carlos Nasher <carlos.nasher at googlemail.com> wrote:

> 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.nasher at gmail.com
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at 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.