# Re: [R] problem with optim and integrate

From: Uwe Ligges <ligges_at_statistik.tu-dortmund.de>
Date: Fri, 21 Mar 2008 19:55:00 +0100

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,

fc <- function(theta, alpha){

numint <- function(z, theta){

```     (theta * dnorm(z, mean = theta, sd = theta) +
theta * dnorm(z, mean = theta, sd = theta) +
(1 - (theta + theta)) *
dnorm(z, mean = theta, sd = theta))^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 + theta)) 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*dnorm(z,mean=theta,sd=theta)+theta*dnorm(z,mean=theta,sd=theta)
>
> +(1-(theta+theta))*dnorm(z,mean=theta,sd=theta))^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))
>
> -----------------------------------------------------------------------------------
>