From: Ken Kelley <KKIII_at_indiana.edu>

Date: Tue 11 Oct 2005 - 14:38:09 EST

Date: Tue 11 Oct 2005 - 14:38:09 EST

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@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.htmlReceived on Tue Oct 11 14:43:57 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:40:41 EST
*