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-boun...@r-project.org
> [mailto:r-help-boun...@r-project.org] On Behalf Of ch_dinesh
> Sent: Monday, May 30, 2011 10:31 AM
> To: r-help@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@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.
>
______________________________________________
R-help@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.
