R help - Oct 2005 - Sometimes having problems finding a minimum using optim(), optimize(), and nlm() (while searching for noncentral F parameters)

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

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
>

Ken Kelley <KKIII at Indiana.Edu> wrote
>
>
> 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.
[SNIP]
> 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.
>
Perhaps uniroot() might be better
 Low.Lim.NC.F <- function(Lambda, alpha.lower, F.value, df.1, df.2) {
    pf(q=F.value, df1=df.1, df2=df.2, ncp=Lambda) - (1-alpha.lower)
 }
 uniroot(Low.Lim.NC.F, interval=c(0,50000), alpha.lower=0.025, df.1=5,
 df.2=200, F.value=12)

If you plot your example, the gradients are much steeper from above than
below, so nlm() works when the starting value is say 50.  In addition,
underflow in pf for more extreme values of lambda affect these more
general search algorithms.


| David Duffy (MBBS PhD)                                         ,-_|\
| email: davidD at qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v