5
$\begingroup$

I'm trying to compute this function, which sums to 1 for $0 < \alpha < 1$ and $k \rightarrow \infty$

  f =(1 - α) * Sum[(α*(ν - k))^ν/Exp[Log[ν!]]/ E^(α*(ν - k)), {ν, 0, k}]

It works when this is substituted...

 f/. {α -> 0.5, k -> 20} = 1.0

But when I try higher values of $\alpha$ (i.e., 0.95), then there appears to be factorial overflow/underflow. For example,

 f/. {α -> 0.95, k -> 40} = 1200.43

which is incorrect. The function should sum to 1.

How do you handle this using Mathematica?

$\endgroup$

2 Answers 2

9
$\begingroup$

Try this:

$Version
(*"12.1.0 for Microsoft Windows (64-bit) (March 14, 2020)"*)

Block[{$MaxExtraPrecision = 5000}, N[f /. {\[Alpha] -> 99/100, k -> 4000}, 100]]

(*0.99999999999999999999999999999999998629286165036357732187047671477498799163903665856827728579195409935*)

Increase $MaxExtraPrecision if you receive overflow/underflow.

$\endgroup$
1
  • $\begingroup$ Nice solution! Thanks $\endgroup$
    – PiE
    Commented Apr 4, 2020 at 19:53
5
$\begingroup$

Why do you believe this series sums to $1$? That claim is clearly false for $k = 1$, $2$, or $3$.

$$f(\alpha, k) = (1-\alpha ) \sum _{\nu =0}^k \frac{e^{-\alpha (\nu -k)} (\alpha (\nu -k))^{\nu }}{\nu !}$$

f[\[Alpha]_, k_] := Simplify[ 
  (1 - \[Alpha])*
  Sum[ (\[Alpha]*(\[Nu] - k))^
    \[Nu]/Exp[Log[\[Nu]!]]/
    E^(\[Alpha]*(\[Nu] - k)), 
    {\[Nu], 0, k}
  ] 
]

(Why "Exp[Log[\[Nu]!]]" when "\[Nu]!" would suffice?)

f[\[Alpha], 1]
(*  -(E^\[Alpha]*(-1 + \[Alpha]))  *)

which is clearly not $1$. More support:

Plot[f[\[Alpha], 1], {\[Alpha], 0, 1}]

f(alpha,1)

f[\[Alpha], 2]
Plot[f[\[Alpha],2], {\[Alpha], 0, 1}]
f[\[Alpha], 3]
Plot[f[\[Alpha],3], {\[Alpha], 0, 1}]

(*  E^\[Alpha]*(-1 + \[Alpha])*(-E^\[Alpha] + \[Alpha])  *)

f(alpha,2)

(*  -(E^\[Alpha]*(-1 + \[Alpha])*(2*E^(2*\[Alpha]) - 4*E^\[Alpha]*\[Alpha] + \[Alpha]^2))/2  *)

f(alpha,3)

The correct question is not "Why isn't this $1$?", because it isn't. The correct question is "Why did you ever get $1$ and how can we get correct non-$1$ answers?"

When you give the input "0.5" you force the number to only have one decimal of precision. So your input cripples the precision of the computation. Let's use better input and then try using more and more internal digits of precision.

f[0.5, 20]
N[f[1/2, 20]]
N[f[1/2, 20], 10]
N[f[1/2, 20], 20]
N[f[1/2, 20], 40]

(*  1.  *)
(*  1.000000000  *)
(*  0.99999999999192806326  *)
(*  0.9999999999919280632566119872342247496372  *)

So your first example is not $1$. There's no recovering from low-precision input.

What about your second example?

Plot[f[\[Alpha], 40],{\[Alpha],0,1}]

f(alpha, 40)

Plot internally numerically computes with relatively low precision and usually, that's fine. But $f(\alpha, 40)$ is computed by catastrophic cancellation in an alternating series of terms. For $\alpha = 19/20 = 0.95$ and $k = 40$, the terms range over $67$ orders of magnitude, so rather a lot of internal precision is required to compute these values.

With[{k = 40, \[Alpha] = 19/20},
  Table[N[
    (E^(-\[Alpha] (-k + \[Nu])) 
    (\[Alpha] (-k + \[Nu]))^\[Nu])/\[Nu]!, 
    20], {\[Nu], 0, k}]
]

(*  {3.1855931757113756220*10^16, -4.5645583886373022450*10^17, 
     3.1046682821270707438*10^18, -1.3337541279910105247*10^19, 
     4.0622295739967888631*10^19, -9.3340156653846143830*10^19, 
     1.6811033091204642577*10^20, -2.4342450042655433475*10^20, 
     2.8841669458885993646*10^20, -2.8312430370213576032*10^20, 
     2.3231556136930412425*10^20, -1.6032265394721459332*10^20, 
     9.3427750895301949039*10^19, -4.6079626182780457491*10^19, 
     1.9249643753359022086*10^19, -6.8068870388418987405*10^18, 
     2.0335373130901240145*10^18, -5.1161163977674949641*10^17, 
     1.0790630006201254205*10^17, -1.8966858018675564914*10^16, 
     2.7576716749377697991*10^15, -3.2861911204615024150*10^14, 
     3.1738286541005634875*10^13, -2.4509224372254223950*10^12, 
     1.4887321436217583874*10^11, -6.9730422544208907612*10^9, 
     2.4583303769637718372*10^8, -6.3323063242935942845*10^6, 
     114855.58252487235285, -1400.2672731086127563, 
     10.810497030509802187, -0.048881621973443540706, 
     0.00011654704546238566397, -1.2661551118342602259*10^-7, 
     5.0705963061585947764*10^-11, -5.4068270460994059133*10^-15, 
     8.9535107693414320584*10^-20, -8.4773911192564840061*10^-26, 
     5.0035227874237341214*10^-34, -1.7148071562535009841*10^-47, 0 }  *)

Plot[f[\[Alpha], 40], {\[Alpha], 0, 1}, 
  WorkingPrecision -> 100, PlotRange -> All]

f(alpha, 40)

So, $f(\alpha, 40)$ also isn't constantly $1$. (It's also not $1200.43$.)

f[0.95, 40]
(*  1200.43  *)

N[f[19/20, 40]]
N[f[19/20, 40], 20]
N[f[19/20, 40], 40]
N[f[19/20, 40], 80]

(*  -36172.3
    0.98347497258026635978
    0.983474972580266359782310302180233513950298347497258026635978231030218023351395023669509725484400437731716599125986976662  *)

Also, using more internal precision cannot rescue a no-precision input.

N[f[0.95, 40], 80]
(*  1200.43  *)

If you are committed to using decimal notation for input, you need to specify more precision.

N[f[0.95`80, 40], 80]
(*  0.9834749725802663597823103021802335139502366950972548440044  *)

The "`80" asserts that this is $0.95$ with $80$ digits of precision and consequently, the computation is able to give a result that is not numerical garbage.

$\endgroup$
3
  • 1
    $\begingroup$ The series is the ss probabilities of an M/D/1 queue so they should sum to 1 for $k \rightarrow \infty$. $\endgroup$
    – PiE
    Commented Apr 5, 2020 at 8:55
  • 1
    $\begingroup$ @PiE : Twice you write the sum is $1$ in your question, unconditionally on $k$. I've also more thoroughly answered your question and explained why you were having numerical precision trouble. $\endgroup$ Commented Apr 5, 2020 at 8:59
  • $\begingroup$ Yes. You had a great answer. Thank you. $\endgroup$
    – PiE
    Commented Apr 5, 2020 at 9:44

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