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.
