# Round-off error

Assuming $$0. I built these two functions

A[k_,q_]:=Sum[PDF[BinomialDistribution[i, q], k]*PDF[ZipfDistribution[n], i], {i, k, Infinity}]


and

B[k_,q_]:=Piecewise[{
{(1/Zeta[n + 1])*PolyLog[n + 1, 1 - q], k == 0},
{(1/(Zeta[n + 1]*k!))*(q/(1 - q))^k*Sum[CoefficientList[Product[x - j, {j, 1, k - 1}],x][[-i]]* PolyLog[n + i - k, 1 - q], {i, 1, k}], k > 0}}]


Functions A and B are mathematically equivalent for $$k\geq1$$. (Function B is faster than function A.)

When I compute $$A[25,0.6]$$ and $$B[25,0.6]$$ with $$n=3$$, I start to get a little difference between the results $$(5.61124*10^{-7},5.61132*10^{-7})$$. So, the question is, why is this difference? Does anyone know how to avoid this round-off problem?

• Mor me: B[25, 0.6] give me: 11.09 not: 5.61124*10^-7 ? – Mariusz Iwaniuk Aug 14 '20 at 9:50
• @MariuszIwaniuk It looks like there are some characters missing in the definition, note the big empty spaces in e.g. n 1. – Lukas Lang Aug 14 '20 at 11:32
• It's normal for different forms of the same mathematical expression to give different numerical approximations, especially when including things like infinite sums. For your particular case, you can use exact numbers: A[25, 6/10] == B[25, 6/10] // FullSimplify and N[A[25, 6/10], 20]. In general, see the many questions on this site about increasing precision of the various Mathematica functions, such as Sum – Lukas Lang Aug 14 '20 at 11:37
• @MariuszIwaniuk I'm sorry, you are right. I just edited the post with the correction. – Slash020 Aug 14 '20 at 13:54
• Exact solvers (like Sum) sometimes have difficulty with inexact input (e.g. 0.6 instead of 6/10). The number 5.61124*10^-7 is the result for A[] returned by NSum. Possibly Sum failed and called NSum as a backup. – Michael E2 Aug 14 '20 at 14:11

This calculates the symbolic infinite sum A[] for a given integer n and caches the result. The result is used to compute subsequent values of A[k, q, n].

ClearAll[A, B];
A[k_, q_, n_Integer?Positive] := (
DownValues[A] = Prepend[
DownValues[A],
HoldPattern[A[kk_, qq_, n]] :> Evaluate@
Sum[
PDF[BinomialDistribution[i, qq], kk]*
PDF[ZipfDistribution[n], i],
{i, kk, Infinity}]];
A[k, q, n])

B[k_, q_, n_Integer?Positive] :=
Piecewise[{{(1/Zeta[n + 1])*PolyLog[n + 1, 1 - q],
k == 0}, {(1/(Zeta[n + 1]*k!))*(q/(1 - q))^k*
Sum[CoefficientList[Product[x - j, {j, 1, k - 1}], x][[-i]]*
PolyLog[n + i - k, 1 - q], {i, 1, k}], k > 0}}];


Difference on floating-point input:

A[25, 0.6, 3] - B[25, 0.6, 3]

(*  3.27923*10^-12  *)


The exact diffference is zero (it takes a few seconds to evaluate the HypergeometricPFQ that appear in A[k, q, 3] on exact input):

Simplify[A[25, 6/10, 3] - B[25, 6/10, 3]]

(* 0  *)

N[A[25, 6/10, 3] - B[25, 6/10, 3], 20]


N::meprec: Internal precision limit \$MaxExtraPrecision = 50. reached while evaluating....

(*  0.*10^-66  *)
`