Besag

I think I will learn a lot from this: Markov Chain Monte Carlo for Statistical Inference (2000) (2001) by Julian Besag

[R] books about MCMC to use MCMC R packages?

---

Computing the mean in R (2.6.0):

 case REALSXP:
     PROTECT(ans = allocVector(REALSXP, 1));
     for (i = 0; i < n; i++) s += REAL(x)[i];
     s /= n;
     if(R_FINITE((double)s)) {
  for (i = 0; i < n; i++) t += (REAL(x)[i] - s);
  s += t/n;
     }
     REAL(ans)[0] = s;
     break;

where t and s were previously initialized to 0.0. (Again I wonder if it would make a noticeable difference to cache the result of REAL(x).) So, in exact math, at the end of the loop t = sum(x) - sum(x) = 0. What does this accomplish in finite precision math?
> y <- seq(-1e6, 1e6, len=1000)
> mean(y)
[1] -1.433398e-11
> sum(y) / 1000
[1] -1.44355e-11
Hm.

And what does Recall (main/eval.c) look like on the inside?

SEXP attribute_hidden do_recall(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    RCNTXT *cptr;
    SEXP s, ans ;
    cptr = R_GlobalContext;
    /* get the args supplied */
    while (cptr != NULL) {
 if (cptr->callflag == CTXT_RETURN && cptr->cloenv == rho)
     break;
 cptr = cptr->nextcontext;
    }
    args = cptr->promargs;
    /* get the env recall was called from */
    s = R_GlobalContext->sysparent;
    while (cptr != NULL) {
 if (cptr->callflag == CTXT_RETURN && cptr->cloenv == s)
     break;
 cptr = cptr->nextcontext;
    }
    if (cptr == NULL)
 error(_("'Recall' called from outside a closure"));

    /* If the function has been recorded in the context, use it
       otherwise search for it by name or evaluate the expression
       originally used to get it.
    */
    if (cptr->callfun != R_NilValue)
 PROTECT(s = cptr->callfun);
    else if( TYPEOF(CAR(cptr->call)) == SYMSXP)
 PROTECT(s = findFun(CAR(cptr->call), cptr->sysparent));
    else
 PROTECT(s = eval(CAR(cptr->call), cptr->sysparent));
    ans = applyClosure(cptr->call, s, args, cptr->sysparent, R_BaseEnv);
    UNPROTECT(1);
    return ans;
}

Renaming a function screws up a recursive function defined in terms of its own name:
> fact <- function(n) if (n > 1) n * fact(n - 1) else 1
> fact(3)
[1] 6
> ffact <- fact
> fact <- function(x) 0 ##Whoops!
> ffact(3)
[1] 0
But not one defined in terms of Recall:
> fact <- function(n) if (n > 1) n * Recall(n - 1) else 1
> fact(3)
[1] 6
> ffact <- fact
> fact <- function(x) 0 ##Whoops!
> ffact(3)
[1] 6
It occurs to me that I've never written a recursive function in R before. It looks like Recall gets a modest amount of use in the base code though.