[R] Root finding problem

[bbolker]

tuki <- function(u, x, a, lambda){
  u <- u+0i
  f <- Re(x-(a*(u)^lambda-(1-(u))^lambda)/lambda)
  f
}

## What I want to do is to find the root, but without specifying the
## interval within which to search for it. I can do it easily in MATLAB
## with fsolve() or fzero() functions.

## BB: how do these functions pick an interval?

## I compared results of running root-finding in MATLAB and R:

## 1) I found (by hand) a suitable interval of u-values for a=5,
## lambda=0.5 and ran uniroot for x from 1 to 5:

curve(tuki(u=x,x=1,a=5,lambda=0.5),from=-1,to=1)
u1 = uniroot(tuki, c(-1,1), x=1, a=5, lambda=0.5)
abline(h=0,col=2)
abline(v=u1$root,col=2)

## The results are the same as I get in MATLAB.

## 2) nlm() does not find the correct value!!!

## BB: nlm is not solving the same problem !!!
##  it's looking for a MINIMUM, not a ROOT.

nlm(tuki, p=0.1, x=1, a=5, lambda=0.5)

## to use nlm to find a minimum:

tukisq <- function(u, x, a, lambda){
  tuki(u,x,a,lambda)^2
}

n1 = nlm(tukisq, p=0.1, x=1, a=5, lambda=0.5)
abline(v=n1$estimate,col=4)

n1$estimate-u1$root
## works fine.

## 3) if I change lambda to 1.5, while keeping a=5, root finding in
## MATLAB returns the following results for x from 1 to 5:

##    0.5134
##    0.7330
##    0.9345
##    1.1289 - 0.0058i
##    1.3085 - 0.0199i

curve(tuki(u=x,x=1,a=5,lambda=1.5),from=-1,to=2,ylim=c(-2,8))
curve(tuki(u=x,x=5,a=5,lambda=1.5),from=-1,to=2,add=TRUE,lty=2)
u2 = uniroot(tuki, c(-1,1), x=1, a=5, lambda=1.5)
abline(h=0,col=2)
abline(v=u2$root,col=2)
u3 = uniroot(tuki, c(-1,2), x=5, a=5, lambda=1.5)
abline(v=u3$root,col=4)

## BB: this seems to work fine.  Did you forget to extend the
##  x range to bracket the root?

## With correctly chosen interval, uniroot() finds correct values only
## for x = 1:3, not for x=4 or x=5. (Obviously, I return real value from
## tuki, but without that uniroot() does not work, returning "Error in
## f(lower, ...) * f(upper, ...) > 0 :  invalid comparison with complex
## values")


## BB: I haven't been through the rest of your code,
## but the approach above again seems to work fine.

n3 = nlm(tukisq, p=0.1, x=5, a=5, lambda=1.5)
abline(v=n3$estimate,col=6)

## also works fine.

-- 
View this message in context: http://www.nabble.com/Root-finding-problem-tf4508016.html#a12865617
Sent from the R help mailing list archive at Nabble.com.