not sure just what you want, but here are some snippets
newton <-
function(lfun, est, tol = 1e-7, niter = 500) {
cscore <- lfun$score(est)
if (abs(cscore) < tol)
return(est)
for (i in 1:niter) {
new <- est - cscore / lfun$d2(est)
cscore <- lfun$score(new)
if (abs(cscore) < tol)
return(new)
est <- new
}
stop("exceeded allowed number of iterations")
}
with the likelihood function
Rmklike <-
function(data) {
n <- length(data)
sumx <- sum(data)
lfun <- function(mu) n * log(mu) - mu * sumx
score <- function(mu) n / mu - sumx
d2 <- function(mu) -n / mu^2
list(lfun = lfun, score = score, d2 = d2)
}
Leeds, Mark (IED) wrote:> In their paper, "Lexical Scope and Statistical Computing", the
authors (
> Gentleman and Ihaka ) go to great length explaining why R's use of
> lexical scoping creates advantages when doing statistical computations.
>
> If anyone has or is familiar with this paper, could they provide the
> main program code for how the "newton" function would be called
in their
> example on page 500 of the paper. The authors are extremely clear in
> their writing and the paper is quite an eye opener for me but it seems
> like lfun somehow needs to be initialized so that it "grabs" the
> environment of "Rmklike".
Rmlike is called to create the likelihood function, and since that
function is defined in the body of Rmlike, it has the evaluation
environment by default. And any return value will have the correct
environment.
> I'm not sure how one would go about doing this so I am wondering what
> the main program that calls "newton" would be
> if there was one. Thanks.
>
So,
data=rexp(10, rate=.3)
lf = Rmklike(data)
newton(lf, .1)
seems to do the trick (and surprisingly still works, as written some
unbelievably long time ago)
Robert
>
>
> Mark
> --------------------------------------------------------
>
> This is not an offer (or solicitation of an offer) to buy/se...{{dropped}}
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>
--
Robert Gentleman, PhD
Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M2-B876
PO Box 19024
Seattle, Washington 98109-1024
206-667-7700
rgentlem at fhcrc.org