# Bug in integral related to beta distribution

I've encountered a problem, which I'm pretty sure is a bug. I've contacted Wolfram's support and they were less than helpful. I'd like to bring the issue up here to get either a confirmation that it is a bug, or a convincing explanation why it is not.

The issue begins when I try to evaluate the following integral:

Integrate[p^1.0 (1 - p)^100.0, {p, 0, 1}]

I get the nonsensical result -1.38685*10^12, whereas the correct result is 0.0000970685.

The Wolfram support representative claimed this is nothing but numerical error caused by using machine-precision floating point, and as such is expected. I claim that this is not a satisfactory answer. Here is a Mathematica notebook that attacks the problem from different angles and, in my opinion, clearly demonstrates there is something else going on:

"Numerical error" shouldn't be used as a fully general excuse for the result of every floating point calculation. Some calculations are numerically unstable, and some are not. If I input 1.0 + 2.0 into Mathematica and got the nonsensical result -3567.321, I wouldn't be satisfied with the explanation "well this is a machine precision calculation, those aren't reliable". I'd demand to get the answer 3.0 (or something very close to it) since there is nothing about the calculation that should trigger wild numerical errors.

It's the same here, the calculation is fairly simple, there is nothing in it that should mess things up. This can be clearly observed from the fact that both numerical integration (with machine precision) works, and symbolic evaluation & substitution works. Whatever method Mathematica chooses to use when presented with the integral directly, it should work.

Furthermore, even if we use custom precision instead of machine precision, we get a nonsensical result and a bogus report of precision. My understanding is that custom precision calculations are tracked for their remaining precision. Any operation that reduces precision (e.g. subtracting two similar numbers) is accounted for, and the final result reflects the remaining precision.

But here we get a result which differs by sign and several orders of magnitude, but still is reported as having high precision:

Integrate[p^SetPrecision[1, 20] (1 - p)^100, {p, 0, 1}]

Gives -1.5682115931362842167*10^8, and

Precision[Integrate[p^SetPrecision[1, 20] (1 - p)^100, {p, 0, 1}]]

Gives 19.936.

Add to that the fact that some results are just weird. For example, adding an irrelevant assumption fixed the problem.

Integrate[p^1.0 (1 - p)^100.0, {p, 0, 1}, Assumptions -> {a > 0}]

Correctly gives 0.0000970685. If this was truly just numerical error, adding an irrelevant assumption wouldn't fix it. There's obviously some bug here which is triggered in some circumstances but not in others.

UPDATE (in response to @michael-e2):

I've thought briefly that maybe it was trying to expand the polynomial and this causes problems, but I didn't think it would try doing that for an inexact exponent. Anyway, this still doesn't explain why we get false precision. For example, suppose I tried

Integrate[Expand[p (1 - p)^100], p] /. p -> SetPrecision[1, 20]

I will get the error "No significant digits are available for display" and the precision of the result will be 0. It recognizes there is not enough precision to numerically calculate it in this form. That's the expected behavior. Increasing the precision works:

Integrate[Expand[p (1 - p)^100], p] /. p -> SetPrecision[1, 100]

Gives a correct result. But in the direct evaluation you get false precision:

Integrate[p^SetPrecision[1, 20] (1 - p)^100, {p, 0, 1}]

And even if that's somehow explained... The point remains that with my original command, there are ways to calculate it that work and ways to calculate it that don't. In some cases it chooses the method that doesn't work and gives a wrong result. This means that when I used Mathematica in practice, my work was disrupted by having to figure out what the problem is and figuring out workarounds. Wouldn't this be considered a bug? Something that Wolfram might want to improve, finding approaches that don't fail so spectacularly?

UPDATE 2:

I originally encountered the problem while working on a slightly older version of Mathematica (some 10.X). Before contacting Wolfram I upgraded to the latest version at the time (11.0) and the results were exactly the same, so I assumed they weren't going to change any time soon.

Taking a cue from @michael-e2, I've upgraded now to 11.1 which has since been released. As he reported, the numerical results are a bit different, but most notably - the results when working with custom precision are now as expected. That is, they do not report a precision which does not exist, and when attempted with a low accuracy (20 or 30 digits) they report 0 precision.

This is an improved behavior of the new version, and I wonder if my report contributed in any way to this change (I doubt it). Together with Michael's discussion of some transformations that Integrate might be trying to do which aren't numerically stable, I guess I can no longer claim this is a hard bug.

I still consider it a soft bug - failing the original command is behavior I did not expect and which could be improved. Michael suggests that Integrate, being a symbolic solver, does not mix well with numeric arguments. But the fact that it is symbolic is exactly why I find this so surprising - what I expected it to do is to solve the general integral symbolically, and then substituting the values numerically.

Of course I can manually have Integrate evaluate the symbolic integral and then substitute, but I'd like Mathematica to work for me, not the other way around.

• Hi; M just does not seem to like that ^1.0 jazz. Integrate[p (1 - p)^100.0, {p, 0, 1}] yields 0.0000970685 as expected. – bobbym Apr 16 '17 at 12:50
• It is generally a bad idea to use inexact numbers with Integrate, or any other functions meant for symbolic computation ... – Szabolcs Apr 16 '17 at 13:07
• (1) If you expand the integrand as a polynomial in the approximate term p^1.0, i.e., with PolynomialReduce[(1 - p)^100. u, p^1. - u, {p^1.`}, CoefficientDomain -> Rationals], you'll certainly get numerical issues. (2) The use of Assumptions throws the integration into an improper solver which rewrites the (1-p)^100.0 in terms of MeijerG[{{}, {101.}}, {{0}, {}}, p], which it can integrate without numerical issues. -- As a symbolic solver, Integrate will make symbolic transformations not guaranteed to be numerically stable. (I'm on V11.1 and get slightly different results.) – Michael E2 Apr 16 '17 at 14:18
• Many symbolic transformations of polynomials are numerically unstable with inexact coefficients. Your tests assume certain ones are done. Unless you happen to hit upon the exact ones used by Integrate, you're likely to get wildly differing results, no? As I implied before, it's generally unwise to use a symbolic solver on a numerical problem. – Michael E2 Apr 16 '17 at 15:32
• What version are you using? I get "no sig. digs." with Integrate[p^SetPrecision[1, 20] (1 - p)^100, {p, 0, 1}]. – Michael E2 Apr 16 '17 at 15:32