evaluating function over seq of values
The reason you are having trouble with using an *apply function is
that f_like does
not have an argument 'N', so the N it uses is the N from the
environment in which
f_like was defined, .GlobalEnv, not one you might set in *apply's FUN argument.
Hence, make N an argument to f_like and use it in *apply. I like
vapply since it gives
you error checking and predictable results.
f_like2 <- function(par, N) {
p1 <- par[1];
p2 <- par[2];
p3 <- par[3];
p4 <- par[4];
lfactorial(N)-lfactorial(N-79) +
(30*log(p1)+(N-30)*log(1-p1)) +
(15*log(p2)+(N-15)*log(1-p2)) +
(22*log(p3)+(N-22)*log(1-p3)) +
(45*log(p4)+(N-45)*log(1-p4)) }
N <- seq(80, 200, 1)
like_mod <- vapply(N,
FUN = function(Ni) {
optim(c(0.2,0.2,0.2,0.2), function(par) f_like2(par, N=Ni),
method = "L-BFGS-B", lower=c(0.005,0.005,0.005,0.005),
upper=c(0.990,0.995,0.995,0.995),
hessian = TRUE,control=list(fnscale=-1))$value
},
FUN.VALUE=0.0)
plot(N, like_mod)
datNew <- cbind(count = seq_along(N), N = N, like_mod = like_mod) #
like your 'dat'
Bill Dunlap
TIBCO Software
wdunlap tibco.com
On Mon, Feb 13, 2017 at 8:41 AM, Evan Cooch <evan.cooch at gmail.com> wrote:
The MWE (below) shows what I'm hoping to get some help with:
step 1\ specify the likelihood expression I want to evaluate using a
brute-force grid search (not that I would do this in practice, but it is a
convenient approach to explain the idea in a class...).
step 2\ want to evaluate the likelihood at each of a sequence of values of N
(for this example, seq(80,200,1)).
step 3\ take all of the values of the likelihood at each value for N, and
(say) plot them
I'm sure there is a clever way to vectorize all this, but my token attempts
at wrestling sapply into submission have failed me here. In my MWE, I use a
simple loop, which has the advantages of working, and being completely
transparent as to what it is doing. For teaching purposes, this is perhaps
fine, but I'm curious as to how I could accomplish the same thing avoiding
the loop.
Thanks in advance...
-----<MWE code below>-----
# ML estimation by simple grid search
rm(list=ls())
library("optimx")
# set up likelihood function
f_like <- function(par) {
p1 <- par[1];
p2 <- par[2];
p3 <- par[3];
p4 <- par[4];
lfactorial(N)-lfactorial(N-79) +
(30*log(p1)+(N-30)*log(1-p1)) +
(15*log(p2)+(N-15)*log(1-p2)) +
(22*log(p3)+(N-22)*log(1-p3)) +
(45*log(p4)+(N-45)*log(1-p4)) }
# do the otimization using optimx nested in a loop (works, but guessing
there is an
# easier way using lapply or some such...)
count <- 1;
dat <- matrix(c(0,0,0),length(seq(80,200,1)),3)
for (N in seq(80,200,1)) {
results_optx <- optim(c(0.2,0.2,0.2,0.2), f_like,
method = "L-BFGS-B", lower=c(0.005,0.005,0.005,0.005),
upper=c(0.990,0.995,0.995,0.995),
hessian = TRUE,control=list(fnscale=-1))
like_mod <- results_optx$value;
dat[count,1] <- count;
dat[count,2] <- N;
dat[count,3] <- like_mod;
count=count+1;
}
plot(dat[,2],dat[,3],type="l",bty="n", xlim=c(79,205), yaxs =
"i",main="likelihood profile",xlab="N", ylab="Likelihood")
______________________________________________ R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.