5
$\begingroup$

I've encountered some problematic behaviour with HypergeometricU. I have a probability distribution on the positive integers that takes the following form after integrating out a few degrees of freedom:

prob[α_, γ_, T_, k_] := ((T)^k Gamma[k + α] Gamma[α + γ]
    HypergeometricU[k + α, 1 + k - γ, T]) / (k! Gamma[α] Gamma[γ])

For many ranges of parameters, this works nicely, but if you plug in something like α = 0.145, γ = 1.71, and T = 53, you'll see that this (incorrectly) evaluates to zero for 14 <= k <= 128. This is somewhat "annoying" (not my original choice of words), given that I would like to perform maximum likelihood / Bayesian analysis on a set of data, but keep running up against Indeterminate outputs whenever my log-likelihood function collides with these erroneous results. (And as luck would have it, the maximum likelihood values for my data set look like they fall in the "broken" zone.)

Now, I can go back to the integral definition of HypergeometricU and attempt to perform that numerically, say with something naive like:

myU[a_?NumericQ, b_?NumericQ, z_?NumericQ] := 
    Re[1/Gamma[a] NIntegrate[Exp[-z t] t^(a - 1) (1 + t)^(b - a - 1), {t, 0, ∞}]]

but this makes my log-likelihood function run s...l...o...w, taking on the order of 30s per call. In addition, the functional form of the integral can be problematic for numerical integration -- as a gets large, the function becomes a very tall, but very narrow spike. (I've tried several different options in NIntegrate, but none quiet the error messages that abound whenever I try applying NMaximize to my log-likelihood.)

Any suggestions from the community? Or are many fun hours reading Abramowitz and Stegun, looking up, checking, rederiving, and ultimately coding up appropriate asymptotic forms in my future? (I can do this, but part of my reason for using Mathematica is so that I don't have to do this ...)

Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

I think it is a working precision problem because you work with big numbers (for Factorial, Gamma and HypergeometricU these numbers are big).

Therefore, you can simply increase the precision

Nprob[α_, γ_, T_, k_] := prob @@ SetPrecision[{α, γ, T, k}, 100]

prob[0.145, 1.71, 53, 100]
Nprob[0.145, 1.71, 53, 100]

0.
0.00024978428152868438468357163500888181009355521204244571578248072896\
   52983233226993080280022314532640

However, it is a problem that the precision of HypergeometricU is poor with MachinePrecision.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Looks like a bug has been identified... $\endgroup$
    – Igor Rivin
    Sep 18, 2014 at 12:04
  • 3
    $\begingroup$ No, this is not a bug. If one gives approximate numbers to begin with, then numerical error is impossible to avoid. A way around this, in some cases, is to do an exact computation first and numericize later, e.g. In[9]:= N[prob[145/1000, 171/100, 53, 17], 20] Out[9]= 0.0048112713266325461380. But the high precision input approach in the response is, most often, the better way to go. $\endgroup$
    – Daniel Lichtblau
    Sep 18, 2014 at 17:27
  • $\begingroup$ Increasing the precision sorted out my problem, although my calls to the log-likelihood function still took a little over 7 seconds. However, the factor of four speedup was sufficient for my purposes. I would tend to agree with @DanielLichtblau that this is not a bug, per se, although some form of warning about the loss of precision would be appreciated -- it currently fails "silently". As I was cursing at this at 2AM last night, it looked like somebody had screwed up on stitching code for two asymptotic approximations together -- I didn't think to try scaling up `WorkingPrecision'. $\endgroup$
    – Confused-cius
    Sep 19, 2014 at 1:10
  • 1
    $\begingroup$ Regarding warnings and loss of precision, the problem is this: the numeric code is not able to discern loss of precision when working with machine numbers. It can only assess that when working with software floats (aka bignums), because only in software arithmetic are error estimates done. That said, whether there is a good warning system for HypergeometricU is something I do not know. Most of the special functions do a reasonable job at getting appropriate error estimates but I gather there are some edge cases. $\endgroup$
    – Daniel Lichtblau
    Sep 19, 2014 at 15:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.