From: Gabor Grothendieck <ggrothendieck_at_gmail.com>

Date: Tue 11 Oct 2005 - 15:27:43 EST

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 Received on Tue Oct 11 15:38:51 2005

Date: Tue 11 Oct 2005 - 15:27:43 EST

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

- try using the square of the difference rather than the absolute value as your objective so that your objective is differentiable.
- 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@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@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
**>
*

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 Received on Tue Oct 11 15:38:51 2005

*
This archive was generated by hypermail 2.1.8
: Sun 23 Oct 2005 - 18:38:50 EST
*