R help - Sep 2007 - Root finding problem

Hello,

I have a problem finding a root of a function, which I define like this:

tuki <- function(u, x, a, lambda){
if((lambda%%1>0) & u<0) {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.

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:

ex.: uniroot(tuki, c(-1,1), x=1, a=5, lambda=0.5)

The results are the same as I get in MATLAB.

2) nlm() does not find the correct value!!!
ex.: nlm(tuki, p=0.1, x=1, a=5, lambda=0.5)

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

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")

There is noone here to ask about this problem, and I just don't know
what to do. :-( Again, what I want is to find the root of the
above-mentioned function, for arbitrary values of x, a, and lambda,
and just specifying one value for u as a starting guess. Could anyone
please give me a suggestion of how to do that? Thanks in advance.

In MATLAB I do the following:

function u=tukeysolve(a,lambda,xgrid)
u=zeros(length(xgrid),1);
old=0.001; % just a guess
for i=1:length(xgrid)
  x=xgrid(i);
  u(i) =
fsolve(@tuki,old,optimset('Display','Off','tolf',5e-8),x,a,lambda);
  old=u(i); % use previous value as new start value
end

function f=tuki(u,x,a,lambda)
f = x - (a * u^lambda - (1-u)^lambda)/lambda;

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.