[R] Double Integral Minimization Problem

[Ravi Varadhan]

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.