Calculate sum of probabilities in multinomial model

Question

The question at hand:

On average in every 7th chocoloate egg there is a figure to be drawn from a known list of special figures (e.g 15 distinct pieces) the other draws are on average fails.

With what probability in a given sample size (n) can you at least get 1/2/3/4/5 distinct specific figures out of these 15 figures (duplicates do not matter, though I could ask for it as well)?

The case with at least one distinct specific figure is: $$\small 1-\left(\frac{104}{105}\right)^n$$

For two I get:

$$\small -2^{3 n+1} \left(\frac{13}{105}\right)^n+\left(\frac{103}{105}\right)^n+1$$

For three I get: $$\small 105^{-n} \left(-102^n+3\dot\ 103^n-3\dot\ 104^n+105^n\right)$$

Which looks kind of nice.

For four I used again the modified code:

p1[a_, b_, n_] := 
 n - Sum[Subscript[k, a - i + 1], {i, 1, b}] - ((a - b))

bubu[a_, n_] := 
 Flatten[Append[
   Reverse[Table[{Subscript[k, i - 1], 1, p1[a - 1, a - i, n]}, {i, 1,
       a}]], {l, 0, p1[a - 1, a - 1, n] - Subscript[k, 0]}], 0]

(*Get at least 4 distinct special figures in your sample n, add up \
probability of all applicable possible series*)

calc22[n_] := 
 Sum[(n!/((k0)!*(k1)!*(k2)!*(k3)!*
       l!*(n - l - k0 - k1 - k2 - 
          k3)!))*(1/(105))^(k0)*(1/(105))^(k1)*(1/(105))^(k2)*(1/(105)\
)^(k3)*((11/105)^l)*((90/105)^(n - l - k0 - k1 - k2 - k3)), {k3, 1, 
   n - 3}, {k2, 1, n - k3 - 2}, {k1, 1, n - k2 - k3 - 1}, {k0, 1, 
   n - k1 - k2 - k3}, {l, 0, n - k0 - k1 - k2 - k3}] 

For more I have written the general code:

p1[a_, b_, n_] := 
 n - Sum[Subscript[k, a - i + 1], {i, 1, b}] - ((a - b))

bubu[a_, n_] := 
 Flatten[Append[
   Reverse[Table[{Subscript[k, i - 1], 1, p1[a - 1, a - i, n]}, {i, 1,
       a}]], {l, 0, p1[a - 1, a - 1, n] - Subscript[k, 0]}], 0]


(*Get at least #(A+1) distinct special figures in your sample n, add \
up probability of all applicable possible series*)
(*Choose low sample size <30*)

calc24[n_, A_] := 
 Sum[(n!/(Product[(Subscript[k, A - i])!, {i, 0, A}]*
       l!*((n - l - 
           Sum[Subscript[k, A - i], {i, 0, A}])!)))*(1/(105))^(Sum[
      Subscript[k, A - i], {i, 0, A}])*(((15 - (A + 1))/105)^
     l)*((90/105)^(n - l - (Sum[Subscript[k, A - i], {i, 0, A}]))), 
  Evaluate[Sequence @@ bubu[A + 1, n]]]

Using the Multinomial function you can substitute the last block beneath the comment with these two:

buba[A_, n_] := 
 Flatten[Append[
   Table[(Subscript[k, A - i]), {i, 0, A}], {l, 
    n - l - Sum[Subscript[k, A - i], {i, 0, A}]}], 1]

calc24[n_, A_] := 
 Sum[(Multinomial[
     Evaluate[Sequence @@ buba[A, n]]])*(1/(105))^(Sum[Subscript[k, 
      A - i], {i, 0, A}])*(((15 - (A + 1))/105)^
     l)*((90/105)^(n - l - (Sum[Subscript[k, A - i], {i, 0, A}]))), 
  Evaluate[Sequence @@ bubu[A + 1, n]]]

Somehow it takes too long to solve for larger values of A (number of distinct figures -1) or n (Sample size).

On this site [crossreference] the probability can be calculated but without stating the mechanics.

My number is slightly off after about 10 digits, I am still astonished that the code actually works to some precision.

Am I doing it overly complicated or is there a way to work with probability in complements for this kind of question?

Thank you.

Answer 1

This will give exact answer:

setp[eggs_, pPerEgg_, setSize_, wanted_] := 
 With[{r = Range[1, wanted]}, 
  1 - Tr[Binomial[wanted, r]*(1 - r*pPerEgg/setSize)^eggs*(-1)^(r - 1)]];

Example: 500 eggs collected, 1/10 chance of egg having item, 50 items available, you want 20 particular:

setp[500, 1/10, 50, 20] 

For huge arguments, you can use inexact arguments (machine precision, e.g.) to speed further if you don't require exact result.

As an aside,

setMeanTo[pPerEgg_, setSize_, wanted_] := (setSize*HarmonicNumber[wanted])/pPerEgg

will give you the mean eggs to collect to reach your goal.

Do note the web site ref'd is not accurate past about the tenth digit shown there...