Skip to content

Error in minimizing an integrand using optim

2 messages · ch_dinesh, William Dunlap

ch_dinesh
# 
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.
William Dunlap
# 
You can put print statements into your integrand to see
what is happening:
+     on.exit(print(cbind(x, fx))) ;
+     fx <-
(e*pnorm((qnorm(p)+sqrt(R)*x)/sqrt(1-R))*dnorm(x))^2/dnorm(x+mu)
+     fx
+ }
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
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
-using-optim-tp3561226p3561226.html