[R] Double Integral Minimization Problem

[MVika]

Hello all,

I am trying to minimize a function which contains a double integral, using
"nlminb" for the minimization and "adapt" for the integral. The integral is
over two variables (thita and radiusb) 
and the 3 free parameters I want to derive from the minimization are
counts0, index and radius_eff.

I have used both tasks in the past successfully but this is the first time
I've tried using the functions where the free parameters are contained
within the integral. My code is copied below:

############## Begin of the code
alpha=3.99
vita=4.4
ellip=0.2

majaxis_pix=c(3.0  6.5 10.0 13.5 17.0 20.5 24.0 27.5 31.0 34.5 38.0 41.5
45.0 48.5 52.0 55.5 59.0)
counts=c(2479.117 718.061 298.918 157.963 98.869  63.883 44.524 31.918
23.500 17.877 13.915 11.032 8.881  7.245  5.978  4.972  4.175112)
countserr=c(80.085 12.181  4.247  3.148  1.963 1.279  0.448  0.242 0.097 
0.055 0.034  0.022  0.015  0.011  0.0082 0.006  0.0043)


intensity=function(x,counts0,index,radius_eff){
	
thita=x[1]
radiusb=x[2]
	
counts2=(1-ellip)*(vita-1)/(pi*alpha^2)*counts0*radiusb*exp(-(2*index-0.324)*(radiusb/radius_eff)^(1/index))*(1+(((majaxis_pix^2)+(radiusb^2)-2*majaxis_pix*radiusb*cos(-thita)+(ellip^2-2*ellip)*((radiusb*sin(thita))^2))/(alpha^2)))^(-vita)
return(counts2)
}
	
	
sersic=function(p)
{
	
counts0=p[1]
index=p[2]
radius_eff=p[3]

value1=adapt(2,c(0,0),c(2*pi,200),functn=intensity,minpts=1000,maxpts=NULL,eps=0.01,counts0=counts0,index=index,radius_eff=radius_eff)

test=value1$value


f=sum(((counts-test)/countserr)^2)
return(f)
}


out<-nlminb(c(16000.0, 3.0,10.0), sersic) 
##############End of the code

Any suggestions are welcome,
Thank you,
Marina





-- 
View this message in context: http://n4.nabble.com/Double-Integral-Minimization-Problem-tp1475150p1475150.html
Sent from the R help mailing list archive at Nabble.com.