# Possible bug / numerical issues with HypergeometricU — any suggestions for a fast workaround?

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

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.

• Looks like a bug has been identified... – Igor Rivin Sep 18 '14 at 12:04
• 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. – Daniel Lichtblau Sep 18 '14 at 17:27
• 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'. – Confused-cius Sep 19 '14 at 1:10
• 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. – Daniel Lichtblau Sep 19 '14 at 15:53