R help - Mar 2008 - problem with optim and integrate

Dear all,

I want to min "integrate( (p1*dnorm+p2*dnorm+p3*dnorm)^(1.3))" for p,
mu,
and sigma.

So, I have to estimate 8 parameters(p3=1-p1-p2).

I got this warning-"Error in integrate(numint, lower = -Inf, upper = Inf) :
non-finite function value."

My questions are

How could I fix it? I tried to divide into several intervals and sum up, but
I got same message.


My code is 

-----------------------------------------------------------------------------------

alpha=1.3; theta0=c(0.5,0,sqrt(10),0.25,-.3,sqrt(0.05),.3,sqrt(0.05))

fc = function(theta){

  numint = function(z)
  {
	    
(theta[1]*dnorm(z,mean=theta[2],sd=theta[3])+theta[4]*dnorm(z,mean=theta[5],sd=theta[6])
             
+(1-(theta[1]+theta[4]))*dnorm(z,mean=theta[7],sd=theta[8]))^alpha
  }


  inte = integrate(numint,lower=-Inf, upper=Inf)$value
  inte
}

optim(theta0, fc, method="L-BFGS-B",hessian=T,
lower=c(rep(c(0,-Inf,0),2),-Inf,0),upper=c(rep(c(1,Inf,4),2),Inf,4))

-----------------------------------------------------------------------------------

Thank you in advance.

Kathryn lord
-- 
View this message in context:
http://www.nabble.com/problem-with-optim-and-integrate-tp16149860p16149860.html
Sent from the R help mailing list archive at Nabble.com.

kathie wrote:
> Dear all,
> 
> I want to min "integrate( (p1*dnorm+p2*dnorm+p3*dnorm)^(1.3))"
for p, mu,
> and sigma.
> 
> So, I have to estimate 8 parameters(p3=1-p1-p2).
> 
> I got this warning-"Error in integrate(numint, lower = -Inf, upper =
Inf) :
> non-finite function value."
> 
> My questions are
> 
> How could I fix it? I tried to divide into several intervals and sum up,
but
> I got same message.
Kathryn,

at first, please make your code readable and clarify it. I tried as follows:

fc <- function(theta, alpha){
   numint <- function(z, theta){
     (theta[1] * dnorm(z, mean = theta[2], sd = theta[3]) +
      theta[4] * dnorm(z, mean = theta[5], sd = theta[6]) +
      (1 - (theta[1] + theta[4])) *
           dnorm(z, mean = theta[7], sd = theta[8]))^alpha
   }
   integrate(numint, lower = -Inf, upper = Inf, theta = theta)$value
}

theta0 <- c(0.5, 0, sqrt(10), 0.25, -0.3, sqrt(0.05), 0.3, sqrt(0.05))

optim(theta0, fc, method = "L-BFGS-B", hessian = TRUE,
       lower = c(rep(c(0,-Inf,0), 2), -Inf, 0),
       upper = c(rep(c(1,Inf,4), 2), Inf, 4),
       alpha = 1.3)


Using some debugging tools (such as a simple browser() call or setting 
options(error=recover)), you will find that numint() can be NaN, because 
(1-(theta[1] + theta[4])) can be negative and hence the sum can be 
negative and a negative value to the power of alpha is not defined.

You may need to penalize for such cases.

Best,
Uwe Ligges


> 
> My code is 
> 
>
-----------------------------------------------------------------------------------
> 
> alpha=1.3; theta0=c(0.5,0,sqrt(10),0.25,-.3,sqrt(0.05),.3,sqrt(0.05))
> 
> fc = function(theta){
> 
>   numint = function(z)
>   {
> 	    
>
(theta[1]*dnorm(z,mean=theta[2],sd=theta[3])+theta[4]*dnorm(z,mean=theta[5],sd=theta[6])
>              
> +(1-(theta[1]+theta[4]))*dnorm(z,mean=theta[7],sd=theta[8]))^alpha
>   }
> 
> 
>   inte = integrate(numint,lower=-Inf, upper=Inf)$value
>   inte
> }
> 
> optim(theta0, fc, method="L-BFGS-B",hessian=T,
> lower=c(rep(c(0,-Inf,0),2),-Inf,0),upper=c(rep(c(1,Inf,4),2),Inf,4))
> 
>
-----------------------------------------------------------------------------------
> 
> Thank you in advance.
> 
> Kathryn lord