[R] Sometimes having problems finding a minimum using optim(), optimize(), and nlm() (while searching for noncentral F parameters)

Gabor Grothendieck ggrothendieck at gmail.com
Tue Oct 11 07:27:43 CEST 2005


I haven't look at your code but here a couple of things to try:

1. try using the square of the difference rather than the absolute value
as your objective so that your objective is differentiable.

2. your objective function may be relatively flat in which case it will be
difficult to get a precise answer.   plot your objective function to see and
try transforming the variable being optimized, e.g. 1/lambda, and then
plotting that to see if its less flat in the region of interest.


On 10/11/05, Ken Kelley <KKIII at indiana.edu> wrote:
> Hi everyone.
>
> I have a problem that I have been unable to determine either the best
> way to proceed and why the methods I'm trying to use sometimes fail. I'm
> using the pf() function in an optimization function to find a
> noncentrality parameter that leads to a specific value at a specified
> quantile. My goal is to have a general function that returns the
> noncentrality parameter that leads to a given value at a defined
> quantile. For example, with 5 and 200 degrees of freedom, what
> noncentrality parameter has at its .975 quantile a value of 4 (it is
> 3.0725 by the way)? The code I've written, using three different methods
> works great at times, but at other times it fails (sometimes all
> sometimes not). It isn't even that the functions I'm trying to write
> fail, but the reason they sometimes fail and sometimes do not is what is
> really bothering me; I simply don't understand why the functions at
> times stop the iterative process of minimization and return what the
> function believes to be a successful convergence value (e.g., optim()
> sometimes returns a 0 stating successful convergence when it clearly is
> not).
>
> I'm using three function [optim(), optimize(), and nlm()] to try and
> accomplish the same goal (which was stated above). I believe that they
> should all return the same value, and at times they do just that, but at
> other times the methods return inappropriate results. I'll paste my code
> that illustrates an example where all is well and one where things fail.
>
> Is there are easier way to do what I'm trying to accomplish? The analog
> in SAS of what I'm trying to come up with is FNONCT.
>
> #Begin code
> ##################################################################
> # Define necessary values.
> F.value <- 4
> tol <- 1e-8
> df.1 <- 5
> df.2 <- 200
> alpha.lower <- .025
> maxit<-1000
>
> # The function to be minimized. Here we are looking for the noncentral
> # value, 'Lambda', that has at its .975 quantile 'F.value'.
> Low.Lim.NC.F <- function(Lambda, alpha.lower, F.value, df.1, df.2)
> {
> abs(pf(q=F.value, df1=df.1, df2=df.2, ncp=Lambda) - (1-alpha.lower))
> # This will be near zero when an appropriate Lambda value is found.
> # The Lambda value that leads to a solution of zero is the noncentrality
> # value that has at its .975 quantile a value of 'F.value'.
> }
>
> # Use the quantile from a central F distribution as a minimum.
> LL.0 <- qf(p=alpha.lower, df1=df.1, df2=df.2)
>
> optim(par=LL.0, fn=Low.Lim.NC.F,
> method="L-BFGS-B", # Others return the same result usually.
> lower=LL.0, upper = Inf, control = list(maxit=maxit, reltol=1e-10),
> hessian = FALSE, alpha.lower=alpha.lower, F.value=F.value, df.1=df.1,
> df.2=df.2)
>
> # Try to accomplish the same task with a different R function.
> optimize(f=Low.Lim.NC.F, lower=LL.0, upper=50, maximum=FALSE, tol=tol,
> alpha.lower=alpha.lower, F.value=F.value, df.1=df.1, df.2=df.2)
>
> # Try to accomplish the same task with a different R function.
> nlm(f=Low.Lim.NC.F, p=LL.0, fscale=1,
> print.level = 0, ndigit=12, gradtol = 1e-6,
> stepmax = max(1000 * sqrt(sum((LL.0/10)^2)), 1000),
> steptol = 1e-6, iterlim = 1000, check.analyticals = TRUE,
> alpha.lower=alpha.lower, F.value=F.value, df.1=df.1, df.2=df.2)
>
> # The answer in each case is 3.0725. Thus, a noncentral F with
> # 5 and 200 df with a noncentrality parameter 3.0725 has at its .975
> # quantile a value of 4 (this has been verified in another software).
>
> # But, suppose we triple the F.value to 12 and rerun the code.
>
> F.value <- 12
>
> # The function to be minimized. Here we are looking for the noncentral
> # value, 'Lambda', that has at its .975 quantile 'F.value'.
> Low.Lim.NC.F <- function(Lambda, alpha.lower, F.value, df.1, df.2)
> {
> abs(pf(q=F.value, df1=df.1, df2=df.2, ncp=Lambda) - (1-alpha.lower))
> }
>
> # Use the quantile from a central F distribution as a minimum.
> LL.0 <- qf(p=alpha.lower, df1=df.1, df2=df.2)
>
> optim(par=LL.0, fn=Low.Lim.NC.F,
> method="L-BFGS-B",  # Others return the same result usually.
>  lower=LL.0, upper = Inf, control = list(maxit=maxit, reltol=1e-10),
> hessian = FALSE, alpha.lower=alpha.lower, F.value=F.value, df.1=df.1,
> df.2=df.2)
>
> # Try to accomplish the same task with a different R function.
> optimize(f=Low.Lim.NC.F, lower=LL.0, upper=500, maximum=FALSE, tol=tol,
> alpha.lower=alpha.lower, F.value=F.value, df.1=df.1, df.2=df.2)
>
> # Try to accomplish the same task with a different R function.
> nlm(f=Low.Lim.NC.F, p=LL.0, fscale=1,
> gradtol = 1e-6, stepmax = max(1000 * sqrt(sum((LL.0/10)^2)), 1000),
> steptol = 1e-6, iterlim = 1000, check.analyticals = TRUE,
> alpha.lower=alpha.lower, F.value=F.value, df.1=df.1, df.2=df.2)
>
> # Now only optimize() works and optim() and nlm() both return the
> # same (wrong) answer. Why would the function stop when the
> # minimized value was .025 (when it should stop when the value is
> # very close to zero)?
>
> # But, optimize() isn't always the answer either, because if the
> # upper limit is too large, the function will fail.
> # For example, changing the upper limit of optimize in
> # this example to 1000 leads to a failure.
> ##################################################################
> #End Code
>
> Am I going about this the best, or even a reasonable, way? Are there
> other functions I'm missing that would be more appropriate given what
> I'm trying to do? Any help would most certainly be appreciated.
>
> Thanks,
> Ken
>
> --
> Ken Kelley, Ph.D.
> Inquiry Methodology Program
> Indiana University
> 201 North Rose Avenue, Room 4004
> Bloomington, Indiana 47405
> http://www.indiana.edu/~kenkel
>
> ______________________________________________
> 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
>




More information about the R-help mailing list