I would like to start another thread, but in essence it is a continuation of this thread

http://community.wolfram.com/groups/-/m/t/208255?p_p_auth=OYun4CSwbut would like it to be a main thread too and not get lost in an older post. Anyway, I would like to present my version or way to calculate probabilities for similar types of problems, e.g. birthday problem, balls in urns etc. It uses a M3 Multinomial method see

"The Handbook of Mathematical Functions" by Abramowitz and Stegun (chapter 24 p 831/832 on Combinatorial Analysis).

http://apps.nrbook.com/abramowitz_and_stegun/index.htmlThe program posted is set to m=9 (say people, or balls etc) and n=365 (say days, or urns etc). When x is set to 0 (default) you will get a full list of all partitions. You can then ask questions like how many of those 9 people share a birthday? the partition would be {1,1,1,1,1,1,1,2} read as 7 people have different birthdays and 2 share a birthday. The probability from the table is 0.0912. or what is probability that 3 people and also 2 other persons share a birthday, partition {1,1,1,1,2,3} etc. To ask the typical birthday problem what is the probability that no one shares a birthday would be the partition {1,1,1,1,1,1,1,1,1}. Change the m to 23 and re-run and look at the bottom probability to see it is the exact probability expected, Be careful going too much above m=23 unless you have lots of memory and plenty of time to wait till all results are printed. The value of n isn't too important to that end result, you just get bigger numbers and not more numbers.

Now changing x to 2 will filter out only those partitions that contain a 2, or a pair, or multiple pairs etc. The set of numbers printed under the table are the M3 values, you may wish to check them against the published ones in abramowitz_and_stegun above, though in the book they only show upto M10. Part of the program calculates this number. Even though I use the Multinomial in MMa it is used to calculate the M3 value (A built in M3 multinomil function would be nice). The 1.000.. printed at the beginning of the table is a check indicating the sum of all probabilities is = to 1. any different value would indicate the M3 has been miss-calculated. The value and order is very important.

n = 365; m = 9; x = 0; tt = 0; tot = 0; list = {};

u = IntegerPartitions[m];

Do[u[[q]] = Sort[u[[q]]], {q, 1, Length[u]}];

u = SortBy[u, Length];

runs[li : {__Integer}] := ((Length /@ Split[#])) &[Sort@li]; lt =

Table[Apply[Multinomial, u, {1}]/

Apply[Times, (runs /@ u)!, {1}], {1}]; lt = Flatten[lt];

list = Reap[Do[pr = (lt[[q]] *n!)/((n - Length[u[[q]]])! n^m);

tot = tot + pr;

If[x == 0, tt = tt + pr; Sow[{q, pr, N[pr, 10], u[[q]], N[tt]}]];

If[Count[u[[q]], x] > 0, tt = tt + pr;

Sow[{q, pr, N[pr, 10], u[[q]], N[tt]}]], {q, 1, Length[lt]}]];

list = Flatten[list[[2]], 1];

Print[N[tot, 20]];

PrependTo[list, {"Ptn No", "Prblty Fraction", "Prblty Decimal",

"Partition", "Accum Probability"}]; Print[

Grid[list, Background -> {None, {Yellow}}, Alignment -> Left,

Frame -> All, Spacings -> {1, 1}]];

Print[lt]