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

[Ken Kelley]

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