R-beta: S Compatibility (again)

Luke Tierney luke at stat.umn.edu
Mon Apr 13 17:54:14 CEST 1998


Peter Dalgaard BSA wrote:
> 
> Peter Dalgaard BSA <p.dalgaard at biostat.ku.dk> writes:
> 
> > > Here is a cute example of what can be done in S but not in R.
> > > Make a function for the pdf of an order statistic.
> > > 
> > > > pdf.order <- function(n, r, pfun, dfun) {
> > >   con <- round(exp(lgamma(n + 1) - lgamma(r) - lgamma(n - r + 1)))
> > >   substitute(
> > >     function(x) {
> > >       Fx <- p(x)
> > >       K*Fx^r1*(1 - Fx)^nr*f(x)
> > >     }, 
> > >     list(p = substitute(pfun), f = substitute(dfun), 
> > >          r1 = r-1, nr = n-r, K = con)
> > >   )
> > > }
> > > > pdf.order(9, 5, pnorm, dnorm)
> > > function(x)
> > > {
> > >         Fx <- pnorm(x)
> > >         630 * Fx^4 * (1 - Fx)^4 * dnorm(x)
> > > }
> > > 
> > > The substitute()s to get unevaluated arguments do work but the
> > > one to modify the function definitely does not.  substitute() is
> > > a very different kind of function in R from what it is in S.
> > 
> > I see the problem. I suspect that there's a workaround, though. Will
> > look at it. 
> 
> [a couple of hours later...]
> 
> > pdf.order                    
> function (n, r, pfun, dfun) 
> {
>         con <- round(exp(lgamma(n + 1) - lgamma(r) - lgamma(n - 
>                 r + 1)))
>         eval(eval(expression(substitute(function(x) {
>                 Fx <- p(x)
>                 K * Fx^r1 * (1 - Fx)^nr * f(x)
>         })), list(p = substitute(pfun), f = substitute(dfun), 
>                 r1 = r - 1, nr = n - r, K = con)), .GlobalEnv)
> }
> > pdf.order(9, 5, pnorm, dnorm)
> function (x) 
> {
>         Fx <- pnorm(x)
>         630 * Fx^4 * (1 - Fx)^4 * dnorm(x)
> }
> > pdf.order(9, 5, pnorm, dnorm)(0)
> [1] 0.981772
> 
> This could have been much easier if substitute allowed a named list
> like eval() does. I see no reason why it shouldn't...
> 
> The necessity of the outer eval is more fundamental: R is sometimes
> more stringent in distinguishing between objects and expressions that
> evaluate to them. So without this necessary eval (<hehe..>) you'd get
> an expression (or more precisely, a "call" object) that evaluates to a
> function, not the function itself:
> 
> I.e.
> 
> S:
> > mode(substitute(function()a+b))
> [1] "function"
> 
> R:
> > mode(substitute(function()a+b))
> [1] "call"
> > mode(eval(substitute(function()a+b)))
> [1] "function"
> 
> You also need to eval() it in .GlobalEnv {or sys.frame(sys.parent(2))}
> , otherwise the result will carry with it the environment of
> pdf.order, which is a direct consequence of R's scoping rules.
> 

Maybe we should think about starting an obfuscated S/R code contest?
:-)

The need to make heavy use of substitute or eval is, in my view, not a
strengh of a language. The right tool for this job, in my view at
least, is lexical scoping and closures. R privides these directly, so
a much simpler definition is

pdf.order<-
function (n, r, pfun, dfun) 
{
  con <- round(exp(lgamma(n + 1) - lgamma(r) - lgamma(n - r + 1)))
  function(x) {
    Fx <- pfun(x)
    con * Fx^(r - 1) * (1 - Fx)^(n - r) * dfun(x)
  }
}

> pdf.order(9, 5, pnorm, dnorm)(0)
[1] 0.981772


The free variables con, n, r, pfun and dfun in the returned function refer
to the variables in the defining environment (this is lexical scope).

Since the free variables in the value function are not modified you
can do this in S if you abstact out the closure creation operation
into a function MC,

pdf.order<-
function (n, r, pfun, dfun) 
{
  con <- round(exp(lgamma(n + 1) - lgamma(r) - lgamma(n - r + 1)))
  MC(function(x) {
       Fx <- pfun(x)
       con * Fx^(r - 1) * (1 - Fx)^(n - r) * dfun(x)
     },
     list("con"=con, "pfun"=pfun, "dfun"=dfun, "r"=r, "n"=n))
}

> pdf.order(9, 5, pnorm, dnorm)(0)
[1] 0.981772

This isn't as clean as in R but it is clearer in its intent, to me at
least, than the eval/substitute stuff.

The definition of MC (make closure) I use is not based on substitute.
Substitute, no matter how sophisticated, is in my view fundamentally
the wrong tool for this job since it does not understand syntax.  (It
is a reasonable tool in this case since the function is very small,
but in general it is not). Both S and R have syntactic constructs that
can result in shadowing of variables -- substitute cannot understand
these and will mess them up. An alternative is to add bindings for
free variables as additional arguments to the function; in this
example this produces

> pdf.order(9, 5, pnorm, dnorm)
function(x, con = 630,
            pfun = function(q, mean = 0, sd = 1) { ... },
            dfun = function(x, mean = 0, sd = 1) { ... },
            r = 5,
            n = 9)
{
        Fx <- pfun(x)
        con * Fx^(r - 1) * (1 - Fx)^(n - r) * dfun(x)
}

This isn't perfect (real lexical scope is much better) but, in my
opinion, it is a better approach than using substitute.  Here is the
definition of MC I use to implement this.  It is fairly convoluted,
but once you have it and understand conceptually what it does, which
is quite simple, you don't need to look at it again.

> MC
function(f, env = NULL)
{
        env <- as.list(env)
        if(mode(f) != "function")
                stop(paste("not a function:", f))
        if(length(env) > 0 && any(names(env) == ""))
                stop(paste("all arguments are not named:", env))
        fargs <- if(length(f) > 1) f[1:(length(f) - 1)] else NULL
        fbody <- f[length(f)]
        cf <- c(fargs, env, fbody)
        mode(cf) <- "function"
        return(cf)
}

Just my 2c worth :-)

luke


-- 
Luke Tierney
University of Minnesota                      Phone:           612-625-7843
School of Statistics                         Fax:             612-624-8868
206 Church Street                            email:      luke at stat.umn.edu
Minneapolis, MN 55455 USA                    WWW:  http://www.stat.umn.edu
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list