1
$\begingroup$

Assuming $0<q<1$. 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?

$\endgroup$
  • 2
    $\begingroup$ Mor me: B[25, 0.6] give me: 11.09 not: 5.61124*10^-7 ? $\endgroup$ – Mariusz Iwaniuk Aug 14 '20 at 9:50
  • $\begingroup$ @MariuszIwaniuk It looks like there are some characters missing in the definition, note the big empty spaces in e.g. n 1. $\endgroup$ – Lukas Lang Aug 14 '20 at 11:32
  • $\begingroup$ 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 $\endgroup$ – Lukas Lang Aug 14 '20 at 11:37
  • $\begingroup$ @MariuszIwaniuk I'm sorry, you are right. I just edited the post with the correction. $\endgroup$ – Slash020 Aug 14 '20 at 13:54
  • 1
    $\begingroup$ 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. $\endgroup$ – Michael E2 Aug 14 '20 at 14:11
1
$\begingroup$

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  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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