[R] Numerical integration again

[casperyc]

Hi all,

Here is an integration function

require(pracma)   # for 'quadinf'

myint=function(j) { 
    quadinf(function(x)
(1/(1+exp(-x)))^j*(1-1/(1+exp(-x)))^(k-j)*dnorm(x,mu,casigma),-Inf,Inf)
}

in any optimization routine. It works fine most of the time but failed with
some particular sets of values, say one of the following:

k=20
mu=-1.978295
casigma=0.008326927

> sapply(0:k,myint)
Error in integrate(f, xa, xb, subdivisions = 512, rel.tol = tol) : 
the integral is probably divergent

I did some investigation on the it and it turns out that it fails when j=7,
ie

> myint(7)
Error in integrate(f, xa, xb, subdivisions = 512, rel.tol = tol) : 
  the integral is probably divergent

All other values from 0 to 20 works fine.

I have even plotted the graph (when j=0:20), most only has a single 'peak',
close to 0 however.


> j=7
> myf=function(x){ 
> (1/(1+exp(-x)))^j*(1-1/(1+exp(-x)))^(k-j)*dnorm(x,mu,casigma) }
> plot(-3:3,myf(-3:3))
> myf(-2)
[1] 1.053024e-07
I just wonder why it does not work when j=7. It could at least give me a
value 1.053024e-07, instead of an error.

> integrate(myf,-Inf,Inf)
0 with absolute error < 0

works fine though.

I am using 'quadinf' to avoid some large parameters that will cause the
function to return an error or return an inaccurate value, discussed in 
http://r.789695.n4.nabble.com/R-numerical-integration-td4500095.html
R-numerical-integration 

Any comments or suggestions?

Thanks!

casper

-----
######################
PhD candidate in Statistics
Big R Fan
Big LEGO Fan
Big sTaTs Fan
######################

--
View this message in context: http://r.789695.n4.nabble.com/Numerical-integration-again-tp4569031p4569031.html
Sent from the R help mailing list archive at Nabble.com.