2
$\begingroup$

I have this discrete probability function that I want to utilize with Mathematica's ProbabilityDistribution[] function

The distribution is:

Pn[0, \[Alpha]_] := 1 - \[Alpha];
Pn[1, \[Alpha]_] := (1 - \[Alpha]) (E^\[Alpha] - 1);
Pn[n_, \[Alpha]_] := (1 - \[Alpha])*
   Sum[E^(j*\[Alpha])*(-1)^(n - 
        j)*(((j*\[Alpha] + n - j)*(j*\[Alpha])^(n - j - 1))/(n - j)!), 
         {j, 0, n}]; 

f[k_] := Pn[k - 1, \[Alpha]]/Sum[Pn[j, \[Alpha]], {j, 0, k - 1}] - 
  Pn[k, \[Alpha]]/Sum[Pn[j, \[Alpha]], {j, 0, k}]

and I would like to use other statistics functions like MomentGeneratingFunction[], ProductDistribution[], etc.

When I type in this:

pmf = ProbabilityDistribution[f[k], {k, 1, 20}];

It just hangs and does nothing...I've waited for over an hour sometimes...I don't even know if I can do this.

Is there a better way to do this? I.e., create a pdf using the function above

$\endgroup$
4
  • 1
    $\begingroup$ Your f[k] seems quite computationally demanding to me. Have you tried calculating a few values for it on its own, to check whether it is the cause of the hang? Is f[k] already normalized? $\endgroup$
    – MarcoB
    Feb 19, 2022 at 14:45
  • $\begingroup$ Yes I have - Sum[f[k], {k, 1, 20}] /. {\[Alpha] -> 0.5} computes quickly... $\endgroup$
    – PiE
    Feb 19, 2022 at 14:49
  • $\begingroup$ I am unclear precisely what you mean by normalization, but the pmf does sum to (almost) 1 - even for {k, 1, 20}. Theoretically, the support is $[1,\infty]$ $\endgroup$
    – PiE
    Feb 19, 2022 at 14:52
  • $\begingroup$ Evaluating Sum[f[k] /. α -> 7/8, {k, 1, 20}] // Simplify // N[#, 20]& shows that either {k, 1, 20} is not a valid range for the distribution or that the distribution needs to be normalized. Note also that ProbabilityDistribution is not a "pmf" but rather a distribution, call it dist[α]. Then the "pmf" is given by PDF[dist[α], x]. $\endgroup$
    – Bob Hanlon
    Feb 19, 2022 at 15:41

1 Answer 1

5
$\begingroup$
Clear["Global`*"]

Start by preventing Pn from evaluating unless n has a nonnegative integer value.

Pn[0, α_] := 1 - α;
Pn[1, α_] := (1 - α) (E^α - 1);
Pn[n_Integer?NonNegative, α_] := (1 - α)*
  Sum[E^(j*α)*(-1)^(n - j)*(((j*α + n - j)*
    (j*α)^(n - j - 1))/(n - j)!), {j, 0, n}]

f[k_] := Pn[k - 1, α]/Sum[Pn[j, α], {j, 0, k - 1}] - 
  Pn[k, α]/Sum[Pn[j, α], {j, 0, k}]

Checking the total probability,

Sum[f[k] /. α -> 7/8, {k, 1, 20}] // Simplify // N[#, 20] &

(* 0.99852321012132571990 *)

Since the total probability for {k, 1, 20} does not equal 1 (and the discrepancy worsens for larger α), the distribution must be normalized.

dist[α_] = ProbabilityDistribution[
   f[k], {k, 1, 20, 1},
   Assumptions -> {0 < α < 1},
   Method -> "Normalize"];

Verifying that the total probability is 1

dist[α] /. ProbabilityDistribution -> Sum // Simplify

(* 1 *)

For symbolic values of α computations are very complicated, e.g.,

LeafCount[Mean[dist[α]] // Simplify] // AbsoluteTiming

(* {1.5569, 3616} *)

LeafCount[MomentGeneratingFunction[dist[α], t] // Simplify] // 
  AbsoluteTiming

(* {31.216, 5629} *)

For numeric values of α the computations are more responsive

Mean[dist[#]] & /@ Range[0.1, 0.9, 0.1] // AbsoluteTiming

(* {1.21928, {1.10051, 1.20419, 1.31467, 1.43658, 1.57661, 1.74568, 1.96436, 
  2.27823, 2.77035}} *)

Variance[dist[#]] & /@ Range[0.1, 0.9, 0.1] // AbsoluteTiming

(* {2.97234, {0.101521, 0.212579, 0.345012, 0.51682, 0.760525, 1.14449, 1.83958, 
  3.37093, 7.15482}} *)

MomentGeneratingFunction[dist[0.5], t] // AbsoluteTiming

(* {0.12198, 
 E^t (0.606531 + 0.264003 E^t + 0.0910935 E^(2 t) + 0.0274169 E^(3 t) + 
    0.00784581 E^(4 t) + 0.0022263 E^(5 t) + 0.000632653 E^(6 t) + 
    0.000180012 E^(7 t) + 0.0000512389 E^(8 t) + 0.0000145858 E^(9 t) + 
    4.15209*10^-6 E^(10 t) + 1.18197*10^-6 E^(11 t) + 
    3.36467*10^-7 E^(12 t) + 9.57826*10^-8 E^(13 t) + 
    2.72672*10^-8 E^(14 t) + 7.75901*10^-9 E^(15 t) + 
    2.19825*10^-9 E^(16 t) + 6.3983*10^-10 E^(17 t) + 
    1.56433*10^-10 E^(18 t) - 6.36646*10^-12 E^(19 t))} *)
$\endgroup$
1
  • $\begingroup$ Thanks so much! $\endgroup$
    – PiE
    Feb 19, 2022 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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