[R] 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.
##
[R] Sometimes having problems finding a minimum using optim(), optimize(), and nlm() (while searching for noncentral F parameters)
Ken Kelley <[EMAIL PROTECTED]> 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,5), 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: [EMAIL PROTECTED] 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
__
R-help@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
Re: [R] Sometimes having problems finding a minimum using optim(), optimize(), and nlm() (while searching for noncentral F parameters)
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 <[EMAIL PROTECTED]> 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 wit
