evaluating function over seq of values

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

evaluating function over seq of values

Evan Cooch
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")

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: evaluating function over seq of values

Bert Gunter-2
The apply() family of functions **are** loops (at the interpreted
level). They are **not vectorized** (looping at the C level). Their
typical virtue is in code clarity and (sometimes) the utiity of the
return structure, not greater efficiency. Sometimes they are a bit
faster, sometimes a bit slower, than *well-written* explicit loops.

Note that in your likelihood function, N can be a vector of values, so
you can compute the likelihood for all values of N and just access the
value you want via subscripting rather than repeatedly computing it
for different N's.  If all you want to do is maximize over a grid of
values, I don't know why you need optim() at all -- but I probably
just don't get what you're doing.

Cheers,
Bert




Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Mon, Feb 13, 2017 at 8:41 AM, Evan Cooch <[hidden email]> 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")
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: evaluating function over seq of values

Evan Cooch

> Note that in your likelihood function, N can be a vector of values, so
> you can compute the likelihood for all values of N and just access the
> value you want via subscripting rather than repeatedly computing it
> for different N's.

OK -- that is the part I'm stuck at - pointers to how to do precisely
that appreciated (I don't use R a lot -- and it shows ;-)

>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: evaluating function over seq of values

R help mailing list-2
In reply to this post by Evan Cooch
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 <[hidden email]> 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")
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.