[R] optim() argument scoping: passing parameter values into user's subfunction

Bert Gunter gunter.berton at gene.com
Wed Sep 9 18:29:45 CEST 2009


Inline below.

Bert Gunter
Genentech Nonclinical Biostatistics

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Ben Bolker
Sent: Tuesday, September 08, 2009 7:29 PM
To: r-help at r-project.org
Subject: Re: [R] optim() argument scoping: passing parameter values into
user's subfunction




Szumiloski, John wrote:
> 
> Dear useRs,
> 
> I have a complicated function to be optimized with optim(), and whose
> parameters are passed to another function within its evaluation.  This
> function allows for the parameters to enter as arguments to various
> probability distribution functions.
> 
> However, I am violating some scoping convention, as somewhere within the
> hierarchy of calls a variable is not visible.  I will give a genericized
> example here.
> 
>> myFxn <- function(parms, Y, phi, <other args>) {<body of function>}
> ### I want to optimize this over its first argument
> 
>>  optim(par=<numeric(>2)>, fn=myFxn, ### end of named args, next are
> all in "..."
>              Y=<data>,  <other args>, 
>              phi=function(r) pnorm(r, mean=parms[2], sd=parms[1])
>              )

This is complicated, but I believe this is what's happening (corrections
from wiser folks would be appreciated if I misstate the situation below):

parms is visible only within the scope of myFxn; by lexical scoping, when
phi is called within myFxn, it goes looking for parms within the environment
it was _defined_, which is optim's, and doesn't find it. Hence the error
message.  Ben's solution below works because it explicitly makes myFxn phi's
closure, thereby giving it access to parms. 

An alternative is to pass parms as a named argument to myFxn. Then it will
look for parms in the _calling_ environment, which will be that of myFxn,
where it will find parms. That is, modify Ben's "original" version to:

myFxn2 <- function(parms, Y,phi,X) {
 ## phi(1) ## not necessary
  sum((Y-(parms[1]+parms[2]*X))^2)
}

set.seed(1001)
X = 1:10
Y = 1+2*X+rnorm(10)
  
optim(par=c(1,2), fn=myFxn2,
      phi = function(r,parms) pnorm(r, mean=parms[2], sd=parms[1]),
      Y=Y, X=X)


This works just fine.

Given the complexity of the original query, I'm not sure this helps, but
maybe ...

------------------


> Error in pnorm(r, mean = parms[2], sd = parms[1]) : 
>   object 'parms' not found
>> debugger()        ## options(error=expression(dump.frames())) in
> .RProfile
> Message:  Error in pnorm(r, mean = parms[2], sd = parms[1]) : 
>   object 'parms' not found
> Available environments had calls:
> 1: optim(par = <numeric(>2)>, fn = myFxn, Y=<data>, <other args>, ..
> 2: function (par) 
> 3: fn(par, ...)
> 4: ifelse(<logical vector>, phi(Y), <other stuff within myFxn>
> 5: phi(Y)
> 6: pnorm(r, mean = parms[2], sd = parms[1])
> 
> Now, when using the debugger in environments 1 and 2, I can see the
> value of 'par' correctly.  I cannot access environment 4 as it just
> returns the original error message.  Trying to access environment 3
> gives the (to me) cryptic  'Error in get(.obj, envir =
> dump[[.selection]]) :  argument "..." is missing, with no default' and
> returns to the top level without debugging.  
> 
> I will try to explain to the best of my ability what I think is
> happening here.  Environments 2 and 3 are from the first lines of
> optim(), where it is building an internal function to evaluate the
> candidate parameter values.  When accessing environment 3, it seems like
> when it fills out the "..." argument of fn(), it is passing
> phi=function(r) pnorm(r, mean=parms[2], sd=parms[1]) but upon trying to
> evaluate the variable 'parms', it cannot see it in the search path.
> When actually running the original call, 'parms' is apparently not
> evaluated yet, but is once the pnorm call is hit.  It appears the
> 'parms' variable is being evaluated before the fn(par) is evaluated into
> myFxn(parms=par).   
> 
> A point which is probably, but not certainly, irrelevant:  myFxn() has
> "..." as its final argument, so as to pass tuning arguments to
> integrate().  The function being integrated contains phi(), as well as
> other stuff, of a dummy variable.  The calls in the debugging tree,
> however, are *not* those involving integrate().  
> 
> It would probably be possible to include some
> substitute/eval/as.name/etc. constructions within the function myFxn, in
> order to avoid this problem, but as myFxn already involves numeric
> integrations whose integrand involves the optimized parameters
> themselves, computing on the language at each step of the optimization
> seems like a bad idea.  
> 
> My question: is there a straightforward and efficient way of
> constructing this function and optimizing it with optim(), allowing for
> the argument phi to pass an arbitrary distribution function whose
> parameters are the global ones being optimized?  In particular, in the
> third environment of the debugger tree, is there a way to force the
> fn(par, ...) ----> myFxn(parms=par, ...) evaluation before the "..." get
> evaluated?
> 
> Thanks, John
> 
> John  Szumiloski,  Ph.D.
> 
> Senior Biometrician
> Biometrics Research
> WP53B-120
> Merck Research Laboratories
> P.O. Box 0004
> West Point, PA 19486-0004
> 
>> (215) 652-7346 (PH)
>> (215) 993-1835 (FAX)
>> 
>> 
> #################################################
> ####  obligatory session info
> ####
> ####  Windows XP 2002 sp3, customized by corporate IT
> 
>> version
>                _                                          
> platform       i386-pc-mingw32                            
> arch           i386                                       
> os             mingw32                                    
> system         i386, mingw32                              
> status         Patched                                    
> major          2                                          
> minor          9.2                                        
> year           2009                                       
> month          09                                         
> day            05                                         
> svn rev        49600                                      
> language       R                                          
> version.string R version 2.9.2 Patched (2009-09-05 r49600)
> 
>> search()
>  [1] ".GlobalEnv"        "package:geometry"  "package:datasets"
> "<mylib1>"           "<mylib2>"          
>  [6] "package:VGAM"      "package:stats4"    "package:Design"
> "package:Hmisc"     "package:boot"     
> [11] "package:splines"   "package:MASS"      "package:nnet"
> "package:utils"     "package:stats"    
> [16] "package:graphics"  "package:grDevices" "<mylib3>"
> "<mylib4>"           "<mylib5>"           
> [21] "<mylib6>"              "package:methods"   "Autoloads"
> "package:base"     
> 
> Notice:  This e-mail message, together with any attachme...{{dropped:15}}
> 
> ______________________________________________
> R-help at 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.
> 
> 


Will the following tweak work for you?  It's not very realistic
but it might solve the problem.

## new version
myFxn <- function(parms, Y,
                  phi=function(r) pnorm(r, mean=parms[2], sd=parms[1]),
                  X) {
  phi(1)
  sum((Y-(parms[1]+parms[2]*X))^2)
}

set.seed(1001)
X = 1:10
Y = 1+2*X+rnorm(10)
  
optim(par=c(1,2), fn=myFxn, 
      Y=Y, X=X)

### original

myFxn2 <- function(parms, Y,
                  phi,
                  X) {
  phi(1)
  sum((Y-(parms[1]+parms[2]*X))^2)
}

set.seed(1001)
X = 1:10
Y = 1+2*X+rnorm(10)
  
optim(par=c(1,2), fn=myFxn2,
      phi = function(r) pnorm(r, mean=parms[2], sd=parms[1]),
      Y=Y, X=X)


-- 
View this message in context:
http://www.nabble.com/optim%28%29-argument-scoping%3A-passing-parameter-valu
es-into-user%27s-subfunction-tp25350121p25357520.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
R-help at 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.




More information about the R-help mailing list