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