Reputation: 85
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
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()
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
Upvotes: 2