Steve Su wrote:
>Splus 6 on Windows currently have an optimize function which can help one
>find the inverse of a function quickly, e.g:
>
>function is x=lambda1 + (p^lambda3 - (1 - p)^lambda4)/lambda2
>
># Define the Objective function:
>qgl.rs.f<-
>function(p, lambda1, lambda2, lambda3, lambda4,x)
>{
>quants <- abs(lambda1 + (p^lambda3 - (1 - p)^lambda4)/lambda2-x)
>quants
>}
>
># Using optimize:
>
>optimize(qgl.rs.f,lambda1=0, lambda2=-1, lambda3=-1.5,
>lambda4=3,x=-5,interval=c(0,1))$minimum
>
>This works very nicely but it is fairly slow if you want to evalue for
>x<-seq(-10,3,length=1000), the syntax is as follows:
>
># vectorised optimize
>
>apply(matrix(seq(-10,3,length=1000)),1,function(k)
>optimize(qgl.rs.f,lambda=0, lambda2=-1, lambda3=-1.5,
>lambda4=3,x=k,interval=c(0,1))$minimum)
>
>I am wondering if there is a quicker way to achieve this result? I am not
>entirely sure how Splus find the inverse of a normal density for example, as
>the programs are coded in .Internal but if there is a function which can be
>used generally in Splus then it will be useful to know this.
(1) It would be more appropriate to use uniroot() than optimize()
to solve an equation.
optimize() can converge happily to an incorrect value, when
the function is not monotone, or when a solution is not possible.
For example, if you try using optimize() to solve cos(theta)=99, it
would tell you the solution is theta=0.
This happens in the example above, for values of x greater than -1,
where no solution is possible in the range 0 < p < 1.
(2) There is a "solveNonlinear" function defined in the help file
for nlmin() (in newer versions of S+). It is intended for solving
multivariate problems, rather than vectorized univariate problems.
It can also converge to a false solution, as noted in the latest help
file.
(3) There is a vectorized function, inverseFunction(),
in the beta resampling library (see www.insightful.com/Hesterberg/bootstrap)
It uses the following syntax:
f _ function(p, lambda1, lambda2, lambda3, lambda4){
lambda1 + (p^lambda3 - (1 - p)^lambda4)/lambda2
}
inverseFunction(f,
seq(-10, -1, length=1000), # values to solve for
initial = c(.2, .8), # couple of values within the domain
lambda1 = 0, lambda2 = -1, lambda3 = -1.5, lambda4 = 3)
# Need to specify initial values, because f is undefined
# at one of the two default initial values, 0.
# Make the problem well-defined, by omitting impossible values greater
# than -1.
(4) The method used to find the inverse of the normal cdf is problem-specific.
Tim Hesterberg
========================================================
| Tim Hesterberg Research Scientist |
| timh@insightful.com Insightful Corp. |
| (206)802-2319 1700 Westlake Ave. N, Suite 500 |
| (206)802-2500 (fax) Seattle, WA 98109-3044, U.S.A. |
| www.insightful.com/Hesterberg |
========================================================
|