0
$\begingroup$

Considering binomial expansion

$1 = (p+(1-p))^n = \sum_{k=0}^n p^n (1-p)^{n-k},$

which is true for any real value of $p$, I'm getting diverging results when I numerically evaluate the above even with very high precision. The code I use is a simple one, say for $n=50$ and $p=10$:

SetPrecision[Sum[Binomial[n, k] p^k (1. - p)^(n - k) /. 
{n -> 50., p -> 10.}, {k,0., 50.}], 1000] // N

For the above computation, each $p^k$ requires at most 50 digits of precision, and same with $(1-p)^k$ and $\binom{n}{k}$ which is at most $2^n$, so in total a couple of hundred digits should be plenty but even with 1000 digits I'm getting $-1.48779 \times 10^{49}$ instead of 1. What am I doing wrong, and how can I change my code to get the correct result 1?

$\endgroup$

closed as off-topic by Daniel Lichtblau, m_goldberg, garej, MarcoB, Bob Hanlon Jun 25 '17 at 14:45

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Daniel Lichtblau, m_goldberg, garej, MarcoB, Bob Hanlon
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 3
    $\begingroup$ The computation as written is done at machine precision, before SetPrecision is applied. Compare with bsum[n_, p_] := Sum[Binomial[n, k] p^k (1 - p)^(n - k), {k, 0, n}] and evaluating bsum[50, 10], bsum[50, 10.] and bsum[50, SetPrecision[10, 80]]. $\endgroup$ – ilian Jun 23 '17 at 20:04
  • $\begingroup$ Thanks, but your example evaluates the sum symbolically to 1 and then applies the parameters (doing nothing). I gave this example to make a point for a much more complex one (that cannot be done symbolically but suffers from the same precision issue as the one I gave). Could you write an example that indeed forces Mathematica to use numerical computation but also gives the correct answer 1.0000? Thanks! $\endgroup$ – MCH Jun 23 '17 at 21:05
  • 2
    $\begingroup$ I think the result of bsum[50, 10.] shows there is no such thing as first evaluating the sum symbolically to 1... $\endgroup$ – ilian Jun 23 '17 at 21:11
  • $\begingroup$ Yes, seems so. I'll try this trick on the more complex calculation and see if it still gives nonsense. Thanks again. $\endgroup$ – MCH Jun 23 '17 at 21:18
  • 1
    $\begingroup$ Maybe you want something like this: Sum[Binomial[n, k] p^k (1.`100 - p)^(n - k) /. {n -> 50.`100, p -> 10.`100}, {k, 0, 50.`100}] -- Or Hold[Sum[Binomial[n, k] p^k (1. - p)^(n - k) /. {n -> 50., p -> 10.}, {k, 0., 50.}]] /. {x_Real /; x != 0 :> SetPrecision[x, 100], z_ /; z == 0 -> 0} // ReleaseHold $\endgroup$ – Michael E2 Jun 24 '17 at 0:34
0
$\begingroup$

You can only get high procession results when every number in the computation is a high-precison inexact number or is exact. If even one machine precision number, such as 1., gets included, the whole computation gets reduced to machine precision. So try this

With[{pr = 50}, 
  N[Sum[Binomial[n, k] p^k (1 - p)^(n - k) /. {n -> 50, p -> 10}, {k, 0, 50}], pr]]

or better yet

With[{pr = 50}, N[Sum[Binomial[50, k] 10^k (1 - 10)^(50 - k), {k, 0, 50}], pr]]

Both give

1.0000000000000000000000000000000000000000000000000

$\endgroup$
  • $\begingroup$ Thanks. To make sure that Mathematica doesn't "cheat" by symbolically evaluating the sum to the exact value 1 and then numerically estimate just 1, I tried to change 1-p to something like SetPrecision[1, 500] - p. Seems to work. $\endgroup$ – MCH Jun 24 '17 at 18:44

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