Hi,
Is your code reproducible?
I can't even find the package "adapt" on the CRAN repository. I
am not sure what exactlt happened to that package, but do remember seeing
something about it relatively recently in R-help. Are you using an older
version of adapt?
Ravi.
____________________________________________________________________
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvaradhan at jhmi.edu
----- Original Message -----
From: MVika <mv56 at st-andrews.ac.uk>
Date: Tuesday, February 9, 2010 5:41 pm
Subject: [R] Double Integral Minimization Problem
To: r-help at r-project.org
> 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:
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
>
> PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.