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. ______________________________________________ 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.