Charles Yan
Charles Yan

Reputation: 85

"uniroot" doesn't work

Here is the equation to be solved:

sum( (x-miu)/(1+l*(x-miu)) ) == 0

x is a vector

x<- c(0.490239414, -0.041047069, -0.062440582, -0.020759616, -0.084667527,  0.006101447,  0.985401602, -0.665158752,  0.153003354,  0.515112122)

miu=0.1 

I tried to use "uniroot"

a<-uniroot(function(l){sum((x-miu)/(1+l*(x-miu)))},c(-20,20))$root

I found "a" is 8.280825. However, when I want to substitute l into a to have a check, I found

sum((x-miu)/(1+a*(x-miu))) 

is -11257.84, instead of 0.

Upvotes: 0

Views: 731

Answers (1)

IRTFM
IRTFM

Reputation: 263411

I can reproduce your difficulties and will show you how I have attempted to overcome them. First off, you do not want x to be the name for to be those data values. You want them assigned to some other name since R will choose x as the default formal parameter to vary for curve and uniroot, so I chose to reverse x with l and further to make it ll since the letter l is so similar to the numeral 1 so:

ll<- c(0.490239414, -0.041047069, -0.062440582, -0.020759616, -0.084667527,  0.006101447,  0.985401602, -0.665158752,  0.153003354,  0.515112122)
miu=0.1

I tried to plot the function with initial failure:

curve(function(x){sum((ll-miu)/(1+ll*(x-miu)))}, -20,20)
Error in curve(function(x) { : 
  'expr' did not evaluate to an object of length 'n'

That mathematical expression on the top of your question wasn't being evaluated as a summation since the sum function isn't "vectorized". So I decided to use the Vectorize function to make a function that would behave in the manner we expect for calculation ( or at least for plotting) purposes:

Vfunc <- Vectorize( function(x){sum((ll-miu)/(1+x*(ll-miu)))} )
png(); curve(Vfunc, -20,20); abline(h=0, lty=3, col="red"); dev.off()  

enter image description here

So looking at the curve lets you see its asymptotic behavior near each of the points at which one would expect the 'll'-values in the denominator to create a degeneracy at the solutions to (1+x*(ll-miu) == 0, and therefore to decide on appropriate ranges over which to seek solution within the non-degenerate intervals.

> Vfunc(-15)
[1] -0.2580872   # Ooops, not far enough to the left, try another point.
> Vfunc(-18)
[1] 0.7168581  # better
> uniroot(Vfunc, c(-18,-5) )   
$root
[1] -16.72707

$f.root
[1] -3.477722e-06

$iter
[1] 7

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05

> png(); curve(Vfunc, -20,20); abline(h=0, lty=3, col="red"); abline(v= -16.72707, col="blue"); dev.off()

> Vfunc(-16.72707)
[1] -4.13423e-06

enter image description here

Upvotes: 2

Related Questions