0

I am working on a maximum likelihood estimator, and one of the parameters is estimated using the digamma function. I'm trying to use uniroot to solve the equation but am unable to do so. Here's my code:

dig = function(alpha){
    digamma(2 + alpha) - digamma(alpha) - (1/(2+alpha)) + (2/(2+alpha))
}

curve(dig, from = 0, to = 10)
uniroot(dig, lower = 0, upper = 10)

This produces the following error:

Error in uniroot(dig, lower = 0, upper = 10) : f.lower = f(lower) is NA
In addition: Warning messages:
1: In digamma(alpha) : NaNs produced
2: In digamma(alpha) : NaNs produced

The first error sort of makes sense based on the curve, but the second has me stuck. It's entirely possible that I'm misunderstanding how to find the roots of the digamma function, or that there's a numerical package (maybe rootsolve?) in R that could help. Not sure what I'm missing here- any tips would be appreciated. Thanks!

sqlck
  • 87
  • 1
  • 7

1 Answers1

1

Consider the following

curve(dig, from = 0.01, to = 10)
uniroot(dig, lower = 0.01, upper = 10)

enter image description here

Error in uniroot(dig, lower = 0.01, upper = 10) : f() values at end points not of opposite sign

MichaelChirico
  • 33,841
  • 14
  • 113
  • 198
  • MichaelChirico, thanks for this example- I should have checked my function more closely. I find that dig has roots in the intervals (-1, 0) and (0, +Inf). Is there a way to pass open intervals to uniroot? – sqlck Feb 12 '17 at 20:43
  • Just run it on the intervals separately – MichaelChirico Feb 12 '17 at 20:43
  • @sqlck Can you explain why dig has a root on (0, +inf)? I'm not understanding it. I do see the one on (-1, 0). – Robert Dodier Feb 12 '17 at 21:08
  • @RobertDodier I'm not saying it does; I'm saying the root-search problem on two intervals can be split into two separate root search problems – MichaelChirico Feb 12 '17 at 21:09
  • @RobertDodier, you're right- sorry about that. There is a root on (-1,0), and two roots on (0, +inf). I used this code: `curve(dig, from = -1.5, to = 10); abline(h = 0, lty = 3);` – sqlck Feb 12 '17 at 21:28