[R] Error in minimizing an integrand using optim

[William Dunlap]

You can put print statements into your integrand to see
what is happening:

> integrand<-function(x){
+     on.exit(print(cbind(x, fx))) ;
+     fx <-
(e*pnorm((qnorm(p)+sqrt(R)*x)/sqrt(1-R))*dnorm(x))^2/dnorm(x+mu)
+     fx
+ }
> mu<-c(-1)
> output<-optim(par=mu, fx2, method="BFGS")
               x           fx
 [1,]   4.090232 3.923504e-05
 [2,] 236.155401          NaN
 [3,]   3.094523 8.917953e-04
 [4,]  41.389072          NaN
 [5,]   3.116343 8.462806e-04
 [6,]  16.890184 4.150785e-68
 [7,]   3.162696 7.553260e-04
 [8,]   9.828110 4.534241e-24
 [9,]   3.238647 6.225113e-04
[10,]   6.922168 2.484860e-12
[11,]   3.351197 4.599135e-04
[12,]   5.456358 5.109337e-08
[13,]   3.512864 2.878964e-04
[14,]   4.614799 4.190869e-06
[15,]   3.746156 1.366194e-04
Error in integrate(integrand, lower = qnorm(0.999), upper = Inf) : 
  non-finite function value
> integrand(20:41)
       x            fx
 [1,] 20  2.270161e-94
 [2,] 21 1.044060e-103
 [3,] 22 1.766442e-113
 [4,] 23 1.099459e-123
 [5,] 24 2.517470e-134
 [6,] 25 2.120582e-145
 [7,] 26 6.571300e-157
 [8,] 27 7.491228e-169
 [9,] 28  0.000000e+00
[10,] 29  0.000000e+00
[11,] 30  0.000000e+00
[12,] 31  0.000000e+00
[13,] 32  0.000000e+00
[14,] 33  0.000000e+00
[15,] 34  0.000000e+00
[16,] 35  0.000000e+00
[17,] 36  0.000000e+00
[18,] 37  0.000000e+00
[19,] 38  0.000000e+00
[20,] 39  0.000000e+00
[21,] 40           NaN
[22,] 41           NaN
...

For x>=40 you are overflowing and/or underflowing (making a value bigger
than 10^308 or smaller than 10^-308), leading to a 0/0 or Inf/Inf
situation, which leads to a NaN (not a number).  For x>=28 you get
exactly 0.

You can
  (a) modify your integrand to return the appropriate limit (0?)
  when x is so large that you are calculating 0/0
  (b) do some algebra to choose an upper limit that avoids
  the NaN problem (if the NaN's really represent limits of 0)
  (c) do some algebra to reparameterize things so you don't
  try to compute 0/0 for large x.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com  

> -----Original Message-----
> From: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] On Behalf Of ch_dinesh
> Sent: Monday, May 30, 2011 10:31 AM
> To: r-help at r-project.org
> Subject: [R] Error in minimizing an integrand using optim
> 
> Hi,
> 
> Am not sure if my code itself is correct. Here's what am trying to do:
> Minimize integration of a function of gaussian distributed 
> variable 'x' over
> the interval qnorm(0.999) to Inf by changing value of 
> parameter 'mu'. mu is
> the shift in mean of 'x'.
> 
> Code:
> # x follows gaussian distribution
> # fx2 to be minimized by changing values of mu
> # integration to be done over the interval (qnorm(0.999),Inf)
> 
> p<-0.009 #constant
> R<-0.25 # constant
> e<-11 #constant
> 
> integrand<-function(x){
> 		
> (e*pnorm((qnorm(p)+sqrt(R)*x)/sqrt(1-R))*dnorm(x))^2/dnorm(x+mu)}
> 
> fx2<-function(mu)	{
> 	integrate(integrand, lower = qnorm(0.999), upper=Inf)$value}
> 
> mu<-c(-1) #initial value for mu, which needs to be optimized 
> such that fx2
> is minimized
> output<-optim(par=mu, fx2, method="BFGS")
> 
> I get the following error message:
> Error in integrate(integrand, lower = qnorm(0.999), upper = Inf) : 
>   non-finite function value
> 
> If upper is changed to 10, error doesn't appear. However, mu 
> retains its
> value and is not optimized.
> 
> Pls. help.
> Thanks
> Dinesh
> 
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Error-in-minimizing-an-integrand
-using-optim-tp3561226p3561226.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> ______________________________________________
> 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.
>