Sometimes you can recast the equation so you can iterate with a vector. For
instance, if you can write the equation as p = fun(x,p), then you can plug p
into the rhs, get new p values, plug it in again and so on until
convergence. It doesn't always converge, but there may be alternative
recastings that do converge. Here's an example where it does converge
(except for x=0)
lam1 <- 0
lam2 <- 1
lam3 <- 2
lam4 <- 2
fun2 <- function(p) (lam1+(p^lam3-(1-p)^lam4)/lam2) # original equation:
x=fun2(p)
fun <- function(p) (lam2*(x-lam1)+(1-p)^lam4)^(1/lam3) # recasting of the
equation p=fun(x,p)
x <- seq(0,1,.01)
p <- fun(1) # initialize
for(i in 1:1000) p <- fun(p)
round(fun2(p)-x,6)
You can speed things up by looking for convergence and only iterating on
parts that haven't converged yet, e.g.
fun <- function(x,p) (lam2*(x-lam1)+(1-p)^lam4)^(1/lam3) # recasting of the
equation
nc <- rep(T,length(x))
for (i in 1:1000) {
prev <- p
p[nc] <- fun(x[nc],p[nc])
nc <- abs(p-prev) > 1e-6
}
Nick Ellis
CSIRO Marine Research mailto:Nick.Ellis@csiro.au
PO Box 120 ph +61 (07) 3826 7260
Cleveland QLD 4163 fax +61 (07) 3826 7222
Australia http://www.marine.csiro.au
> -----Original Message-----
> From: Steve Su [mailto:s.su@qut.edu.au]
> Sent: Saturday, March 22, 2003 8:44 PM
> To: s-news@lists.biostat.wustl.edu
> Subject: [S] Finding the inverse of a function (not necessarily
> monotonic) through numerical methods.
>
>
> Dear All,
>
> 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.
>
> Thank you for your attention, all comments welcome.
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> **************************************************************
> **************
> **********
>
> Steve Su (s.su@qut.edu.au)
> PhD student.
>
> School of Accountancy
> School of Mathematical Sciences
> Queensland University of Technology
>
> Postal Address: Steve Su, School of Accountancy, QUT, PO Box 2434,
> Brisbane,
> Queensland, Australia, 4000.
>
>
> Phone: +61 7 3864 2017
> Fax: +61 7 3864 1812
> Mobile: 0421 840 586
> .
> _--_|\
> / QUT
> \_.--._/
> v
>
> **************************************************************
> **************
> **********
>
> --------------------------------------------------------------------
> This message was distributed by s-news@lists.biostat.wustl.edu. To
> unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
> the BODY of the message: unsubscribe s-news
>
|