Birthday problem: Difference between revisions
→{{header|Lasso}}: adding lasso version of the birthday problem |
Add Scala implementation |
||
(63 intermediate revisions by 30 users not shown) | |||
Line 1: | Line 1: | ||
[[Category:Probability and statistics]] |
|||
[[Category:Discrete math]] |
|||
{{Wikipedia pre 15 June 2009|pagename=Birthday Problem|lang=en|oldid=296054030|timedate=21:44, 12 June 2009}} |
{{Wikipedia pre 15 June 2009|pagename=Birthday Problem|lang=en|oldid=296054030|timedate=21:44, 12 June 2009}} |
||
{{draft task}} |
|||
{{draft task}}[[Category:Probability and statistics]][[Category:Discrete math]] |
|||
In [[wp:probability theory|probability theory]], the '''birthday problem''', or '''birthday [[wp:paradox|paradox]]''' This is not a paradox in the sense of leading to a [[wp:logic|logic]]al contradiction, but is called a paradox because the mathematical truth contradicts naïve [[wp:intuition (knowledge)|intuition]]: most people estimate that the chance is much lower than 50%. pertains to the [[wp:probability|probability]] that in a set of [[wp:random|random]]ly chosen people some pair of them will have the same [[wp:birthday|birthday]]. In a group of at least 23 randomly chosen people, there is more than 50% probability that some pair of them will both have been born on the same day. For 57 or more people, the probability is more than 99%, and it reaches 100% when the number of people reaches 366 (by the [[wp:pigeon hole principle|pigeon hole principle]], ignoring leap years). The mathematics behind this problem leads to a well-known cryptographic attack called the [[wp:birthday attack|birthday attack]]. |
In [[wp:probability theory|probability theory]], the '''birthday problem''', or '''birthday [[wp:paradox|paradox]]''' This is not a paradox in the sense of leading to a [[wp:logic|logic]]al contradiction, but is called a paradox because the mathematical truth contradicts naïve [[wp:intuition (knowledge)|intuition]]: most people estimate that the chance is much lower than 50%. pertains to the [[wp:probability|probability]] that in a set of [[wp:random|random]]ly chosen people some pair of them will have the same [[wp:birthday|birthday]]. In a group of at least 23 randomly chosen people, there is more than 50% probability that some pair of them will both have been born on the same day. For 57 or more people, the probability is more than 99%, and it reaches 100% when the number of people reaches 366 (by the [[wp:pigeon hole principle|pigeon hole principle]], ignoring leap years). The mathematics behind this problem leads to a well-known cryptographic attack called the [[wp:birthday attack|birthday attack]]. |
||
;Task |
;Task |
||
Using simulation, estimate the number of independent people required in a |
Using simulation, estimate the number of independent people required in a group before we can expect a ''better than even chance'' that at least 2 independent people in a group share a common birthday. Furthermore: Simulate and thus estimate when we can expect a ''better than even chance'' that at least 3, 4 & 5 independent people of the group share a common birthday. For simplicity assume that all of the people are alive... |
||
;Suggestions for improvement |
;Suggestions for improvement |
||
Line 10: | Line 16: | ||
* Converging to the <math>n</math><sup>th</sup> solution using a root finding method, as opposed to using an extensive search. |
* Converging to the <math>n</math><sup>th</sup> solution using a root finding method, as opposed to using an extensive search. |
||
* Kudos (κῦδος) for finding the solution by proof (in a programming language) rather than by construction and simulation. |
* Kudos (κῦδος) for finding the solution by proof (in a programming language) rather than by construction and simulation. |
||
;See also |
;See also |
||
* {{Wolfram|Birthday|Problem}} |
* Wolfram entry: {{Wolfram|Birthday|Problem}} |
||
<br><br> |
|||
=={{header|11l}}== |
|||
{{trans|D}} |
|||
<syntaxhighlight lang="11l">F equal_birthdays(sharers, groupsize, rep) |
|||
V eq = 0 |
|||
L 0 .< rep |
|||
V group = [0] * 365 |
|||
L 0 .< groupsize |
|||
group[random:(group.len)]++ |
|||
I any(group.map(c -> c >= @sharers)) |
|||
eq++ |
|||
R (eq * 100.) / rep |
|||
V group_est = 2 |
|||
L(sharers) 2..5 |
|||
V groupsize = group_est + 1 |
|||
L equal_birthdays(sharers, groupsize, 100) < 50. |
|||
groupsize++ |
|||
L(gs) Int(groupsize - (groupsize - group_est) / 4.) .< groupsize + 999 |
|||
V eq = equal_birthdays(sharers, gs, 250) |
|||
I eq > 50. |
|||
groupsize = gs |
|||
L.break |
|||
L(gs) groupsize - 1 .< groupsize + 999 |
|||
V eq = equal_birthdays(sharers, gs, 50'000) |
|||
I eq > 50. |
|||
group_est = gs |
|||
print(‘#. independent people in a group of #. share a common birthday. (#3.1)’.format(sharers, gs, eq)) |
|||
L.break</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 independent people in a group of 23 share a common birthday. ( 50.6) |
|||
3 independent people in a group of 87 share a common birthday. ( 50.1) |
|||
4 independent people in a group of 187 share a common birthday. ( 50.5) |
|||
5 independent people in a group of 313 share a common birthday. ( 50.4) |
|||
</pre> |
|||
=={{header|Ada}}== |
|||
This solution assumes a 4-year cycle, with three 365-day years and one leap year. |
|||
<syntaxhighlight lang="ada">with Ada.Command_Line, Ada.Text_IO, Ada.Numerics.Discrete_random; |
|||
procedure Birthday_Test is |
|||
Samples: constant Positive := Integer'Value(Ada.Command_Line.Argument(1)); |
|||
-- our experiment: Generate a X (birth-)days and check for Y-collisions |
|||
-- the constant "Samples" is the number of repetitions of this experiment |
|||
subtype Day is integer range 0 .. 365; -- this includes leap_days |
|||
subtype Extended_Day is Integer range 0 .. 365*4; -- a four-year cycle |
|||
package ANDR is new Ada.Numerics.Discrete_Random(Extended_Day); |
|||
Random_Generator: ANDR.Generator; |
|||
function Random_Day return Day is (ANDR.Random(Random_Generator) / 4); |
|||
-- days 0 .. 364 are equally probable, leap-day 365 is 4* less probable |
|||
type Checkpoint is record |
|||
Multiplicity: Positive; |
|||
Person_Count: Positive; |
|||
end record; |
|||
Checkpoints: constant array(Positive range <>) of Checkpoint |
|||
:= ( (2, 22), (2, 23), (3, 86), (3, 87), (3, 88), |
|||
(4, 186), (4, 187), (5, 312), (5, 313), (5, 314) ); |
|||
type Result_Type is array(Checkpoints'Range) of Natural; |
|||
Result: Result_Type := (others => 0); |
|||
-- how often is a 2-collision in a group of 22 or 23, ..., a 5-collision |
|||
-- in a group of 312 .. 314 |
|||
procedure Experiment(Result: in out Result_Type) is |
|||
-- run the experiment once! |
|||
A_Year: array(Day) of Natural := (others => 0); |
|||
A_Day: Day; |
|||
Multiplicity: Natural := 0; |
|||
People: Positive := 1; |
|||
begin |
|||
for I in Checkpoints'Range loop |
|||
while People <= Checkpoints(I).Person_Count loop |
|||
A_Day := Random_Day; |
|||
A_Year(A_Day) := A_Year(A_Day)+1; |
|||
if A_Year(A_Day) > Multiplicity then |
|||
Multiplicity := Multiplicity + 1; |
|||
end if; |
|||
People := People + 1; |
|||
end loop; |
|||
if Multiplicity >= Checkpoints(I).Multiplicity then |
|||
Result(I) := Result(I) + 1; |
|||
-- found a Multipl.-collision in a group of Person_Cnt. |
|||
end if; |
|||
end loop; |
|||
end Experiment; |
|||
package TIO renames Ada.Text_IO; |
|||
package FIO is new TIO.Float_IO(Float); |
|||
begin |
|||
-- initialize the random generator |
|||
ANDR.Reset(Random_Generator); |
|||
-- repeat the experiment Samples times |
|||
for I in 1 .. Samples loop |
|||
Experiment(Result); |
|||
end loop; |
|||
-- print the results |
|||
TIO.Put_Line("Birthday-Test with" & Integer'Image(Samples) & " samples:"); |
|||
for I in Result'Range loop |
|||
FIO.Put(Float(Result(I))/Float(Samples), Fore => 3, Aft => 6, Exp => 0); |
|||
TIO.Put_Line |
|||
("% of groups with" & Integer'Image(Checkpoints(I).Person_Count) & |
|||
" have" & Integer'Image(Checkpoints(I).Multiplicity) & |
|||
" persons sharing a common birthday."); |
|||
end loop; |
|||
end Birthday_Test;</syntaxhighlight> |
|||
{{out}} |
|||
Running the program with a sample size 500_000_000 took about 25 minutes on a slow pc. |
|||
<pre>./birthday_test 500_000_000 |
|||
Birthday-Test with 500000000 samples: |
|||
0.475292% of groups with 22 have 2 persons sharing a common birthday. |
|||
0.506882% of groups with 23 have 2 persons sharing a common birthday. |
|||
0.487155% of groups with 86 have 3 persons sharing a common birthday. |
|||
0.498788% of groups with 87 have 3 persons sharing a common birthday. |
|||
0.510391% of groups with 88 have 3 persons sharing a common birthday. |
|||
0.494970% of groups with 186 have 4 persons sharing a common birthday. |
|||
0.501825% of groups with 187 have 4 persons sharing a common birthday. |
|||
0.495137% of groups with 312 have 5 persons sharing a common birthday. |
|||
0.500010% of groups with 313 have 5 persons sharing a common birthday. |
|||
0.504888% of groups with 314 have 5 persons sharing a common birthday.</pre> |
|||
An interesting observation: |
|||
The probability for groups of 313 persons having 5 persons sharing a common birthday is almost exactly 0.5. Note that a solution based on 365-day years, i.e., a solution ignoring leap days, would generate slightly but significantly larger probabilities. |
|||
=={{header|ALGOL 68}}== |
=={{header|ALGOL 68}}== |
||
{{works with|ALGOL 68|Revision 1}} |
{{works with|ALGOL 68|Revision 1}} |
||
{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-2.6 algol68g-2.6].}} |
{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-2.6 algol68g-2.6].}} |
||
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''.}} |
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''.}} |
||
'''File: Birthday_problem.a68'''< |
'''File: Birthday_problem.a68'''<syntaxhighlight lang="algol68">#!/usr/bin/a68g --script # |
||
# -*- coding: utf-8 -*- # |
# -*- coding: utf-8 -*- # |
||
Line 82: | Line 225: | ||
required common +:= 1 |
required common +:= 1 |
||
FI |
FI |
||
OD # group size #</ |
OD # group size #</syntaxhighlight>'''Output:''' |
||
<pre> |
<pre> |
||
upb year: 365.2500; upb common: 5; upb sample size: 100000; |
upb year: 365.2500; upb common: 5; upb sample size: 100000; |
||
Line 96: | Line 239: | ||
required common: 5; group size: 314; %age of years with required common birthdays: 50.66%; |
required common: 5; group size: 314; %age of years with required common birthdays: 50.66%; |
||
</pre> |
</pre> |
||
=={{header|C}}== |
|||
Computing probabilities to 5 sigmas of confidence. It's very slow, chiefly because to make sure a probability like 0.5006 is indeed above .5 instead of just statistical fluctuation, you have to run the simulation millions of times. |
|||
<syntaxhighlight lang="c">#include <stdio.h> |
|||
#include <stdlib.h> |
|||
#include <string.h> |
|||
#include <math.h> |
|||
#define DEBUG 0 // set this to 2 for a lot of numbers on output |
|||
#define DAYS 365 |
|||
#define EXCESS (RAND_MAX / DAYS * DAYS) |
|||
int days[DAYS]; |
|||
inline int rand_day(void) |
|||
{ |
|||
int n; |
|||
while ((n = rand()) >= EXCESS); |
|||
return n / (EXCESS / DAYS); |
|||
} |
|||
// given p people, if n of them have same birthday in one run |
|||
int simulate1(int p, int n) |
|||
{ |
|||
memset(days, 0, sizeof(days)); |
|||
while (p--) |
|||
if (++days[rand_day()] == n) return 1; |
|||
return 0; |
|||
} |
|||
// decide if the probability of n out of np people sharing a birthday |
|||
// is above or below p_thresh, with n_sigmas sigmas confidence |
|||
// note that if p_thresh is very low or hi, minimum runs need to be much higher |
|||
double prob(int np, int n, double n_sigmas, double p_thresh, double *std_dev) |
|||
{ |
|||
double p, d; // prob and std dev |
|||
int runs = 0, yes = 0; |
|||
do { |
|||
yes += simulate1(np, n); |
|||
p = (double) yes / ++runs; |
|||
d = sqrt(p * (1 - p) / runs); |
|||
if (DEBUG > 1) |
|||
printf("\t\t%d: %d %d %g %g \r", np, yes, runs, p, d); |
|||
} while (runs < 10 || fabs(p - p_thresh) < n_sigmas * d); |
|||
if (DEBUG > 1) putchar('\n'); |
|||
*std_dev = d; |
|||
return p; |
|||
} |
|||
// bisect for truth |
|||
int find_half_chance(int n, double *p, double *dev) |
|||
{ |
|||
int lo, hi, mid; |
|||
reset: |
|||
lo = 0; |
|||
hi = DAYS * (n - 1) + 1; |
|||
do { |
|||
mid = (hi + lo) / 2; |
|||
// 5 sigma confidence. Conventionally people think 3 sigmas are good |
|||
// enough, but for case of 5 people sharing birthday, 3 sigmas actually |
|||
// sometimes give a slightly wrong answer |
|||
*p = prob(mid, n, 5, .5, dev); |
|||
if (DEBUG) |
|||
printf("\t%d %d %d %g %g\n", lo, mid, hi, *p, *dev); |
|||
if (*p < .5) lo = mid + 1; |
|||
else hi = mid; |
|||
if (hi < lo) { |
|||
// this happens when previous precisions were too low; |
|||
// easiest fix: reset |
|||
if (DEBUG) puts("\tMade a mess, will redo."); |
|||
goto reset; |
|||
} |
|||
} while (lo < mid || *p < .5); |
|||
return mid; |
|||
} |
|||
int main(void) |
|||
{ |
|||
int n, np; |
|||
double p, d; |
|||
srand(time(0)); |
|||
for (n = 2; n <= 5; n++) { |
|||
np = find_half_chance(n, &p, &d); |
|||
printf("%d collision: %d people, P = %g +/- %g\n", |
|||
n, np, p, d); |
|||
} |
|||
return 0; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 collision: 23 people, P = 0.508741 +/- 0.00174794 |
|||
3 collision: 88 people, P = 0.509034 +/- 0.00180628 |
|||
4 collision: 187 people, P = 0.501812 +/- 0.000362394 |
|||
5 collision: 313 people, P = 0.500641 +/- 0.000128174 |
|||
</pre> |
|||
=={{header|C++}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="cpp">#include <iostream> |
|||
#include <random> |
|||
#include <vector> |
|||
double equalBirthdays(int nSharers, int groupSize, int nRepetitions) { |
|||
std::default_random_engine generator; |
|||
std::uniform_int_distribution<int> distribution(0, 364); |
|||
std::vector<int> group(365); |
|||
int eq = 0; |
|||
for (int i = 0; i < nRepetitions; i++) { |
|||
std::fill(group.begin(), group.end(), 0); |
|||
for (int j = 0; j < groupSize; j++) { |
|||
int day = distribution(generator); |
|||
group[day]++; |
|||
} |
|||
if (std::any_of(group.cbegin(), group.cend(), [nSharers](int c) { return c >= nSharers; })) { |
|||
eq++; |
|||
} |
|||
} |
|||
return (100.0 * eq) / nRepetitions; |
|||
} |
|||
int main() { |
|||
int groupEst = 2; |
|||
for (int sharers = 2; sharers < 6; sharers++) { |
|||
// Coarse |
|||
int groupSize = groupEst + 1; |
|||
while (equalBirthdays(sharers, groupSize, 100) < 50.0) { |
|||
groupSize++; |
|||
} |
|||
// Finer |
|||
int inf = (int)(groupSize - (groupSize - groupEst) / 4.0f); |
|||
for (int gs = inf; gs < groupSize + 999; gs++) { |
|||
double eq = equalBirthdays(sharers, groupSize, 250); |
|||
if (eq > 50.0) { |
|||
groupSize = gs; |
|||
break; |
|||
} |
|||
} |
|||
// Finest |
|||
for (int gs = groupSize - 1; gs < groupSize + 999; gs++) { |
|||
double eq = equalBirthdays(sharers, gs, 50000); |
|||
if (eq > 50.0) { |
|||
groupEst = gs; |
|||
printf("%d independant people in a group of %d share a common birthday. (%5.1f)\n", sharers, gs, eq); |
|||
break; |
|||
} |
|||
} |
|||
} |
|||
return 0; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 independant people in a group of 23 share a common birthday. ( 50.6) |
|||
3 independant people in a group of 87 share a common birthday. ( 50.2) |
|||
4 independant people in a group of 186 share a common birthday. ( 50.0) |
|||
5 independant people in a group of 313 share a common birthday. ( 50.5)</pre> |
|||
=={{header|D}}== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="d">import std.stdio, std.random, std.algorithm, std.conv; |
|||
/// For sharing common birthday must all share same common day. |
|||
double equalBirthdays(in uint nSharers, in uint groupSize, |
|||
in uint nRepetitions, ref Xorshift rng) { |
|||
uint eq = 0; |
|||
foreach (immutable _; 0 .. nRepetitions) { |
|||
uint[365] group; |
|||
foreach (immutable __; 0 .. groupSize) |
|||
group[uniform(0, $, rng)]++; |
|||
eq += group[].any!(c => c >= nSharers); |
|||
} |
|||
return (eq * 100.0) / nRepetitions; |
|||
} |
|||
void main() { |
|||
auto rng = 1.Xorshift; // Fixed seed. |
|||
auto groupEst = 2; |
|||
foreach (immutable sharers; 2 .. 6) { |
|||
// Coarse. |
|||
auto groupSize = groupEst + 1; |
|||
while (equalBirthdays(sharers, groupSize, 100, rng) < 50.0) |
|||
groupSize++; |
|||
// Finer. |
|||
immutable inf = to!int(groupSize - (groupSize - groupEst) / 4.0); |
|||
foreach (immutable gs; inf .. groupSize + 999) { |
|||
immutable eq = equalBirthdays(sharers, groupSize, 250, rng); |
|||
if (eq > 50.0) { |
|||
groupSize = gs; |
|||
break; |
|||
} |
|||
} |
|||
// Finest. |
|||
foreach (immutable gs; groupSize - 1 .. groupSize + 999) { |
|||
immutable eq = equalBirthdays(sharers, gs, 50_000, rng); |
|||
if (eq > 50.0) { |
|||
groupEst = gs; |
|||
writefln("%d independent people in a group of %s share a common birthday. (%5.1f)", |
|||
sharers, gs, eq); |
|||
break; |
|||
} |
|||
} |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 independent people in a group of 23 share a common birthday. ( 50.5) |
|||
3 independent people in a group of 87 share a common birthday. ( 50.1) |
|||
4 independent people in a group of 187 share a common birthday. ( 50.2) |
|||
5 independent people in a group of 313 share a common birthday. ( 50.3)</pre> |
|||
Run-time about 10.4 seconds with ldc2 compiler. |
|||
Alternative version: |
|||
{{trans|C}} |
|||
<syntaxhighlight lang="d">import std.stdio, std.random, std.math; |
|||
enum nDays = 365; |
|||
// 5 sigma confidence. Conventionally people think 3 sigmas are good |
|||
// enough, but for case of 5 people sharing birthday, 3 sigmas |
|||
// actually sometimes give a slightly wrong answer. |
|||
enum double nSigmas = 3.0; // Currently 3 for smaller run time. |
|||
/// Given n people, if m of them have same birthday in one run. |
|||
bool simulate1(in uint nPeople, in uint nCollisions, ref Xorshift rng) |
|||
/*nothrow*/ @safe /*@nogc*/ { |
|||
static uint[nDays] days; |
|||
days[] = 0; |
|||
foreach (immutable _; 0 .. nPeople) { |
|||
immutable day = uniform(0, days.length, rng); |
|||
days[day]++; |
|||
if (days[day] == nCollisions) |
|||
return true; |
|||
} |
|||
return false; |
|||
} |
|||
/** Decide if the probablity of n out of np people sharing a birthday |
|||
is above or below pThresh, with nSigmas sigmas confidence. |
|||
If pThresh is very low or hi, minimum runs need to be much higher. */ |
|||
double prob(in uint np, in uint nCollisions, in double pThresh, |
|||
out double stdDev, ref Xorshift rng) { |
|||
double p, d; // Probablity and standard deviation. |
|||
uint nRuns = 0, yes = 0; |
|||
do { |
|||
yes += simulate1(np, nCollisions, rng); |
|||
nRuns++; |
|||
p = double(yes) / nRuns; |
|||
d = sqrt(p * (1 - p) / nRuns); |
|||
debug if (yes % 50_000 == 0) |
|||
printf("\t\t%d: %d %d %g %g \r", np, yes, nRuns, p, d); |
|||
} while (nRuns < 10 || abs(p - pThresh) < (nSigmas * d)); |
|||
debug '\n'.putchar; |
|||
stdDev = d; |
|||
return p; |
|||
} |
|||
/// Bisect for truth. |
|||
uint findHalfChance(in uint nCollisions, out double p, out double dev, ref Xorshift rng) { |
|||
uint mid; |
|||
RESET: |
|||
uint lo = 0; |
|||
uint hi = nDays * (nCollisions - 1) + 1; |
|||
do { |
|||
mid = (hi + lo) / 2; |
|||
p = prob(mid, nCollisions, 0.5, dev, rng); |
|||
debug printf("\t%d %d %d %g %g\n", lo, mid, hi, p, dev); |
|||
if (p < 0.5) |
|||
lo = mid + 1; |
|||
else |
|||
hi = mid; |
|||
if (hi < lo) { |
|||
// This happens when previous precisions were too low; |
|||
// easiest fix: reset. |
|||
debug "\tMade a mess, will redo.".puts; |
|||
goto RESET; |
|||
} |
|||
} while (lo < mid || p < 0.5); |
|||
return mid; |
|||
} |
|||
void main() { |
|||
auto rng = Xorshift(unpredictableSeed); |
|||
foreach (immutable uint nCollisions; 2 .. 6) { |
|||
double p, d; |
|||
immutable np = findHalfChance(nCollisions, p, d, rng); |
|||
writefln("%d collision: %d people, P = %g +/- %g", nCollisions, np, p, d); |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 collision: 23 people, P = 0.521934 +/- 0.00728933 |
|||
3 collision: 88 people, P = 0.512367 +/- 0.00411469 |
|||
4 collision: 187 people, P = 0.506974 +/- 0.00232306 |
|||
5 collision: 313 people, P = 0.501588 +/- 0.000529277</pre> |
|||
Output with nSigmas = 5.0: |
|||
<pre>2 collision: 23 people, P = 0.508607 +/- 0.00172133 |
|||
3 collision: 88 people, P = 0.511945 +/- 0.00238885 |
|||
4 collision: 187 people, P = 0.503229 +/- 0.000645587 |
|||
5 collision: 313 people, P = 0.501105 +/- 0.000221016</pre> |
|||
=={{header|Delphi}}== |
|||
{{libheader| System.SysUtils}} |
|||
{{Trans|C}} |
|||
<syntaxhighlight lang="delphi"> |
|||
program Birthday_problem; |
|||
{$APPTYPE CONSOLE} |
|||
uses |
|||
System.SysUtils; |
|||
const |
|||
_DAYS = 365; |
|||
var |
|||
days: array[0..(_DAYS) - 1] of Integer; |
|||
runs: Integer; |
|||
function rand_day: Integer; inline; |
|||
begin |
|||
Result := random(_DAYS); |
|||
end; |
|||
/// <summary> |
|||
/// given p people, if n of them have same birthday in one run |
|||
/// </summary> |
|||
function simulate1(p, n: Integer): Integer; |
|||
var |
|||
index, i: Integer; |
|||
begin |
|||
for i := 0 to High(days) do |
|||
begin |
|||
days[i] := 0; |
|||
end; |
|||
for i := 0 to p - 1 do |
|||
begin |
|||
index := rand_day(); |
|||
inc(days[index]); |
|||
if days[index] = n then |
|||
Exit(1); |
|||
end; |
|||
Exit(0); |
|||
end; |
|||
/// <summary> |
|||
/// decide if the probability of n out of np people sharing a birthday |
|||
/// is above or below p_thresh, with n_sigmas sigmas confidence |
|||
/// note that if p_thresh is very low or hi, minimum runs need to be much higher |
|||
/// </summary> |
|||
function prob(np, n: Integer; n_sigmas, p_thresh: Double; var d: Double): Double; |
|||
var |
|||
p: Double; |
|||
runs, yes: Integer; |
|||
begin |
|||
runs := 0; |
|||
yes := 0; |
|||
repeat |
|||
yes := yes + (simulate1(np, n)); |
|||
inc(runs); |
|||
p := yes / runs; |
|||
d := sqrt(p * (1 - p) / runs); |
|||
until ((runs >= 10) and (abs(p - p_thresh) >= (n_sigmas * d))); |
|||
Exit(p); |
|||
end; |
|||
/// <summary> |
|||
/// bisect for truth |
|||
/// </summary> |
|||
function find_half_chance(n: Integer; var p: Double; var d: Double): Integer; |
|||
var |
|||
lo, hi, mid: Integer; |
|||
label |
|||
reset; |
|||
begin |
|||
reset: |
|||
lo := 0; |
|||
hi := _DAYS * (n - 1) + 1; |
|||
repeat |
|||
mid := (hi + lo) div 2; |
|||
p := prob(mid, n, 3, 0.5, d); |
|||
if p < 0.5 then |
|||
lo := mid + 1 |
|||
else |
|||
hi := mid; |
|||
if hi < lo then |
|||
goto reset; |
|||
until ((lo >= mid) and (p >= 0.5)); |
|||
Exit(mid); |
|||
end; |
|||
var |
|||
n, np: Integer; |
|||
p, d: Double; |
|||
begin |
|||
Randomize; |
|||
writeln('Wait for calculate'); |
|||
for n := 2 to 5 do |
|||
begin |
|||
np := find_half_chance(n, p, d); |
|||
writeln(format('%d collision: %d people, P = %.8f +/- %.8f', [n, np, p, d])); |
|||
end; |
|||
writeln('Press enter to exit'); |
|||
readln; |
|||
end.</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Wait for calculate |
|||
2 collision: 23 people, P = 0,50411600 +/- 0,00137151 |
|||
3 collision: 88 people, P = 0,51060651 +/- 0,00352751 |
|||
4 collision: 188 people, P = 0,50856493 +/- 0,00285022 |
|||
5 collision: 313 people, P = 0,50093354 +/- 0,00031116 |
|||
Press enter to exit</pre> |
|||
=={{header|FreeBASIC}}== |
|||
{{trans|XPL0}} |
|||
<syntaxhighlight lang="freebasic">Function Simulacion(N As Integer) As Integer |
|||
Dim As Integer i, dias(365) |
|||
For i = 0 To 365-1 |
|||
dias(i) = 0 |
|||
Next i |
|||
Dim As Integer R, personas = 0 |
|||
Do |
|||
R = Rnd * 365 |
|||
dias(R) += 1 |
|||
personas += 1 |
|||
If dias(R) = N Then Return personas |
|||
Loop |
|||
End Function |
|||
Dim As Integer N, grupo, t |
|||
For N = 2 To 5 |
|||
grupo = 0 |
|||
For t = 1 To 10000 |
|||
grupo += Simulacion(N) |
|||
Next t |
|||
Print Using "Average of # people in a population of ### share birthdays"; N; Int(grupo/10000) |
|||
Next N |
|||
Sleep</syntaxhighlight> |
|||
=={{header|Go}}== |
|||
{{trans|C}} |
|||
<syntaxhighlight lang="go">package main |
|||
import ( |
|||
"fmt" |
|||
"math" |
|||
"math/rand" |
|||
"time" |
|||
) |
|||
const ( |
|||
DEBUG = 0 |
|||
DAYS = 365 |
|||
n_sigmas = 5. |
|||
WORKERS = 16 // concurrent worker processes |
|||
RUNS = 1000 // runs per flight |
|||
) |
|||
func simulate1(p, n int, r *rand.Rand) int { |
|||
var days [DAYS]int |
|||
for i := 0; i < p; i++ { |
|||
days[r.Intn(DAYS)]++ |
|||
} |
|||
for _, d := range days { |
|||
if d >= n { |
|||
return 1 |
|||
} |
|||
} |
|||
return 0 |
|||
} |
|||
// send yes's per fixed number of simulate1 runs until canceled |
|||
func work(p, n int, ych chan int, cancel chan bool) { |
|||
r := rand.New(rand.NewSource(time.Now().Unix() + rand.Int63())) |
|||
for { |
|||
select { |
|||
case <-cancel: |
|||
return |
|||
default: |
|||
} |
|||
y := 0 |
|||
for i := 0; i < RUNS; i++ { |
|||
y += simulate1(p, n, r) |
|||
} |
|||
ych <- y |
|||
} |
|||
} |
|||
func prob(np, n int) (p, d float64) { |
|||
ych := make(chan int, WORKERS) |
|||
cancel := make(chan bool) |
|||
for i := 0; i < WORKERS; i++ { |
|||
go work(np, n, ych, cancel) |
|||
} |
|||
var runs, yes int |
|||
for { |
|||
yes += <-ych |
|||
runs += RUNS |
|||
fr := float64(runs) |
|||
p = float64(yes) / fr |
|||
d = math.Sqrt(p * (1 - p) / fr) |
|||
if DEBUG > 1 { |
|||
fmt.Println("\t\t", np, yes, runs, p, d) |
|||
} |
|||
// .5 here is the "even chance" threshold |
|||
if !(math.Abs(p-.5) < n_sigmas*d) { |
|||
close(cancel) |
|||
break |
|||
} |
|||
} |
|||
if DEBUG > 1 { |
|||
fmt.Println() |
|||
} |
|||
return |
|||
} |
|||
func find_half_chance(n int) (mid int, p, dev float64) { |
|||
reset: |
|||
lo := 0 |
|||
hi := DAYS*(n-1) + 1 |
|||
for { |
|||
mid = (hi + lo) / 2 |
|||
p, dev = prob(mid, n) |
|||
if DEBUG > 0 { |
|||
fmt.Println("\t", lo, mid, hi, p, dev) |
|||
} |
|||
if p < .5 { |
|||
lo = mid + 1 |
|||
} else { |
|||
hi = mid |
|||
} |
|||
if hi < lo { |
|||
if DEBUG > 0 { |
|||
fmt.Println("\tMade a mess, will redo.") |
|||
} |
|||
goto reset |
|||
} |
|||
if !(lo < mid || p < .5) { |
|||
break |
|||
} |
|||
} |
|||
return |
|||
} |
|||
func main() { |
|||
for n := 2; n <= 5; n++ { |
|||
np, p, d := find_half_chance(n) |
|||
fmt.Printf("%d collision: %d people, P = %.4f ± %.4f\n", |
|||
n, np, p, d) |
|||
} |
|||
}</syntaxhighlight> |
|||
<pre> |
|||
2 collision: 23 people, P = 0.5081 ± 0.0016 |
|||
3 collision: 88 people, P = 0.5155 ± 0.0029 |
|||
4 collision: 187 people, P = 0.5041 ± 0.0008 |
|||
5 collision: 313 people, P = 0.5015 ± 0.0003 |
|||
</pre> |
|||
'''Also based on the C version:''' |
|||
<syntaxhighlight lang="go">package main |
|||
import ( |
|||
"fmt" |
|||
"math" |
|||
"math/rand" |
|||
"runtime" |
|||
"time" |
|||
) |
|||
type ProbeRes struct { |
|||
np int |
|||
p, d float64 |
|||
} |
|||
type Frac struct { |
|||
n int |
|||
d int |
|||
} |
|||
var DaysInYear int = 365 |
|||
func main() { |
|||
sigma := 5.0 |
|||
for i := 2; i <= 5; i++ { |
|||
res := GetNP(i, sigma, 0.5) |
|||
fmt.Printf("%d collision: %d people, P = %.4f ± %.4f\n", |
|||
i, res.np, res.p, res.d) |
|||
} |
|||
} |
|||
func GetNP(n int, n_sigmas, p_thresh float64) (res ProbeRes) { |
|||
res.np = DaysInYear * (n - 1) |
|||
for i := 0; i < DaysInYear*(n-1); i++ { |
|||
tmp := probe(i, n, n_sigmas, p_thresh) |
|||
if tmp.p > p_thresh && tmp.np < res.np { |
|||
res = tmp |
|||
} |
|||
} |
|||
return |
|||
} |
|||
var numCPU = runtime.NumCPU() |
|||
func probe(np, n int, n_sigmas, p_thresh float64) ProbeRes { |
|||
var p, d float64 |
|||
var runs, yes int |
|||
cRes := make(chan Frac, numCPU) |
|||
for i := 0; i < numCPU; i++ { |
|||
go SimN(np, n, 25, cRes) |
|||
} |
|||
for math.Abs(p-p_thresh) < n_sigmas*d || runs < 100 { |
|||
f := <-cRes |
|||
yes += f.n |
|||
runs += f.d |
|||
p = float64(yes) / float64(runs) |
|||
d = math.Sqrt(p * (1 - p) / float64(runs)) |
|||
go SimN(np, n, runs/3, cRes) |
|||
} |
|||
return ProbeRes{np, p, d} |
|||
} |
|||
func SimN(np, n, ssize int, c chan Frac) { |
|||
r := rand.New(rand.NewSource(time.Now().UnixNano() + rand.Int63())) |
|||
yes := 0 |
|||
for i := 0; i < ssize; i++ { |
|||
if Sim(np, n, r) { |
|||
yes++ |
|||
} |
|||
} |
|||
c <- Frac{yes, ssize} |
|||
} |
|||
func Sim(p, n int, r *rand.Rand) (res bool) { |
|||
Cal := make([]int, DaysInYear) |
|||
for i := 0; i < p; i++ { |
|||
Cal[r.Intn(DaysInYear)]++ |
|||
} |
|||
for _, v := range Cal { |
|||
if v >= n { |
|||
res = true |
|||
} |
|||
} |
|||
return |
|||
}</syntaxhighlight> |
|||
{{Out}} |
|||
<pre> |
|||
2 collision: 23 people, P = 0.5068 ± 0.0013 |
|||
3 collision: 88 people, P = 0.5148 ± 0.0028 |
|||
4 collision: 187 people, P = 0.5020 ± 0.0004 |
|||
5 collision: 313 people, P = 0.5011 ± 0.0002 |
|||
</pre> |
|||
=={{header|Hy}}== |
|||
We use a simple but not very accurate simulation method. |
|||
<syntaxhighlight lang="lisp">(import |
|||
[numpy :as np] |
|||
[random [randint]]) |
|||
(defmacro incf (place) |
|||
`(+= ~place 1)) |
|||
(defn birthday [required &optional [reps 20000] [ndays 365]] |
|||
(setv days (np.zeros (, reps ndays) np.int_)) |
|||
(setv qualifying-reps (np.zeros reps np.bool_)) |
|||
(setv group-size 1) |
|||
(setv count 0) |
|||
(while True |
|||
;(print group-size) |
|||
(for [r (range reps)] |
|||
(unless (get qualifying-reps r) |
|||
(setv day (randint 0 (dec ndays))) |
|||
(incf (get days (, r day))) |
|||
(when (= (get days (, r day)) required) |
|||
(setv (get qualifying-reps r) True) |
|||
(incf count)))) |
|||
(when (> (/ (float count) reps) .5) |
|||
(break)) |
|||
(incf group-size)) |
|||
group-size) |
|||
(print (birthday 2)) |
|||
(print (birthday 3)) |
|||
(print (birthday 4)) |
|||
(print (birthday 5))</syntaxhighlight> |
|||
=={{header|J}}== |
|||
Quicky approach (use a population of 1e5 people to get a quick estimate and then refine against a population of 1e8 people): |
|||
<syntaxhighlight lang="j">PopSmall=: 1e5 ?@# 365 |
|||
PopBig=: 1e8 ?@# 365 |
|||
countShared=: [: >./ #/.~ |
|||
avg=: +/ % # |
|||
probShared=: (1 :0)("0) |
|||
: |
|||
NB. y: shared birthday count |
|||
NB. m: population |
|||
NB. x: sample size |
|||
avg ,y <: (-x) countShared\ m |
|||
) |
|||
estGroupSz=: 3 :0 |
|||
approx=. (PopSmall probShared&y i.365) I. 0.5 |
|||
n=. approx-(2+y) |
|||
refine=. n+(PopBig probShared&y approx+i:2+y) I. 0.5 |
|||
assert. (2+y) > |approx-refine |
|||
refine, refine PopBig probShared y |
|||
)</syntaxhighlight> |
|||
Task cases: |
|||
<syntaxhighlight lang="j"> estGroupSz 2 |
|||
23 0.507254 |
|||
estGroupSz 3 |
|||
88 0.510737 |
|||
estGroupSz 4 |
|||
187 0.502878 |
|||
estGroupSz 5 |
|||
313 0.500903</syntaxhighlight> |
|||
So, for example, we need a group of 88 to have at least a 50% chance of 3 people in the group having the same birthday in a year of 365 days. And, in that case, the simulated probability was 51.0737% |
|||
=={{header|Java}}== |
|||
Translation of [[Birthday_problem#Python|Python]] via [[Birthday_problem#D|D]] |
|||
{{works with|Java|8}} |
|||
<syntaxhighlight lang="java">import static java.util.Arrays.stream; |
|||
import java.util.Random; |
|||
public class Test { |
|||
static double equalBirthdays(int nSharers, int groupSize, int nRepetitions) { |
|||
Random rand = new Random(1); |
|||
int eq = 0; |
|||
for (int i = 0; i < nRepetitions; i++) { |
|||
int[] group = new int[365]; |
|||
for (int j = 0; j < groupSize; j++) |
|||
group[rand.nextInt(group.length)]++; |
|||
eq += stream(group).anyMatch(c -> c >= nSharers) ? 1 : 0; |
|||
} |
|||
return (eq * 100.0) / nRepetitions; |
|||
} |
|||
public static void main(String[] a) { |
|||
int groupEst = 2; |
|||
for (int sharers = 2; sharers < 6; sharers++) { |
|||
// Coarse. |
|||
int groupSize = groupEst + 1; |
|||
while (equalBirthdays(sharers, groupSize, 100) < 50.0) |
|||
groupSize++; |
|||
// Finer. |
|||
int inf = (int) (groupSize - (groupSize - groupEst) / 4.0); |
|||
for (int gs = inf; gs < groupSize + 999; gs++) { |
|||
double eq = equalBirthdays(sharers, groupSize, 250); |
|||
if (eq > 50.0) { |
|||
groupSize = gs; |
|||
break; |
|||
} |
|||
} |
|||
// Finest. |
|||
for (int gs = groupSize - 1; gs < groupSize + 999; gs++) { |
|||
double eq = equalBirthdays(sharers, gs, 50_000); |
|||
if (eq > 50.0) { |
|||
groupEst = gs; |
|||
System.out.printf("%d independent people in a group of " |
|||
+ "%s share a common birthday. (%5.1f)%n", |
|||
sharers, gs, eq); |
|||
break; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
}</syntaxhighlight> |
|||
<pre>2 independent people in a group of 23 share a common birthday. ( 50,6) |
|||
3 independent people in a group of 87 share a common birthday. ( 50,4) |
|||
4 independent people in a group of 187 share a common birthday. ( 50,1) |
|||
5 independent people in a group of 314 share a common birthday. ( 50,2)</pre> |
|||
=={{header|Julia}}== |
|||
{{works with|Julia|0.6}} |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="julia">function equalbirthdays(sharers::Int, groupsize::Int; nrep::Int = 10000) |
|||
eq = 0 |
|||
for _ in 1:nrep |
|||
group = rand(1:365, groupsize) |
|||
grset = Set(group) |
|||
if groupsize - length(grset) ≥ sharers - 1 && |
|||
any(count(x -> x == d, group) ≥ sharers for d in grset) |
|||
eq += 1 |
|||
end |
|||
end |
|||
return eq / nrep |
|||
end |
|||
gsizes = [2] |
|||
for sh in (2, 3, 4, 5) |
|||
local gsize = gsizes[end] |
|||
local freq |
|||
# Coarse |
|||
while equalbirthdays(sh, gsize; nrep = 100) < .5 |
|||
gsize += 1 |
|||
end |
|||
# Finer |
|||
for gsize in trunc(Int, gsize - (gsize - gsizes[end]) / 4):(gsize + 999) |
|||
if equalbirthdays(sh, gsize; nrep = 250) > 0.5 |
|||
break |
|||
end |
|||
end |
|||
# Finest |
|||
for gsize in (gsize - 1):(gsize + 999) |
|||
freq = equalbirthdays(sh, gsize; nrep = 50000) |
|||
if freq > 0.5 |
|||
break |
|||
end |
|||
end |
|||
push!(gsizes, gsize) |
|||
@printf("%i independent people in a group of %s share a common birthday. (%5.3f)\n", sh, gsize, freq) |
|||
end</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 independent people in a group of 23 share a common birthday. (0.506) |
|||
3 independent people in a group of 88 share a common birthday. (0.510) |
|||
4 independent people in a group of 187 share a common birthday. (0.500) |
|||
5 independent people in a group of 314 share a common birthday. (0.507)</pre> |
|||
=={{header|Kotlin}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="scala">// version 1.1.3 |
|||
import java.util.Random |
|||
fun equalBirthdays(nSharers: Int, groupSize: Int, nRepetitions: Int): Double { |
|||
val rand = Random(1L) |
|||
var eq = 0 |
|||
for (i in 0 until nRepetitions) { |
|||
val group = IntArray(365) |
|||
for (j in 0 until groupSize) { |
|||
group[rand.nextInt(group.size)]++ |
|||
} |
|||
eq += if (group.any { it >= nSharers}) 1 else 0 |
|||
} |
|||
return eq * 100.0 / nRepetitions |
|||
} |
|||
fun main(args: Array<String>) { |
|||
var groupEst = 2 |
|||
for (sharers in 2..5) { |
|||
// Coarse |
|||
var groupSize = groupEst + 1 |
|||
while (equalBirthdays(sharers, groupSize, 100) < 50.0) groupSize++ |
|||
// Finer |
|||
val inf = (groupSize - (groupSize - groupEst) / 4.0).toInt() |
|||
for (gs in inf until groupSize + 999) { |
|||
val eq = equalBirthdays(sharers, groupSize, 250) |
|||
if (eq > 50.0) { |
|||
groupSize = gs |
|||
break |
|||
} |
|||
} |
|||
// Finest |
|||
for (gs in groupSize - 1 until groupSize + 999) { |
|||
val eq = equalBirthdays(sharers, gs, 50_000) |
|||
if (eq > 50.0) { |
|||
groupEst = gs |
|||
print("$sharers independent people in a group of ${"%3d".format(gs)} ") |
|||
println("share a common birthday (${"%2.1f%%".format(eq)})") |
|||
break |
|||
} |
|||
} |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
Expect runtime of about 15 seconds on a modest laptop: |
|||
<pre> |
|||
2 independent people in a group of 23 share a common birthday (50.6%) |
|||
3 independent people in a group of 87 share a common birthday (50.4%) |
|||
4 independent people in a group of 187 share a common birthday (50.1%) |
|||
5 independent people in a group of 314 share a common birthday (50.2%) |
|||
</pre> |
|||
=={{header|Lasso}}== |
=={{header|Lasso}}== |
||
< |
<syntaxhighlight lang="lasso">if(sys_listunboundmethods !>> 'randomgen') => { |
||
define randomgen(len::integer,max::integer)::array => { |
define randomgen(len::integer,max::integer)::array => { |
||
#len <= 0 ? return |
#len <= 0 ? return |
||
Line 130: | Line 1,189: | ||
'Threshold: '+#threshold+', qty: '+#qty+' - probability: '+#probability+'\r' |
'Threshold: '+#threshold+', qty: '+#qty+' - probability: '+#probability+'\r' |
||
#qty += 1 |
#qty += 1 |
||
^}</ |
^}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 154: | Line 1,213: | ||
Threshold: 5, qty: 314 - probability: 50.000000 |
Threshold: 5, qty: 314 - probability: 50.000000 |
||
</pre> |
</pre> |
||
=={{header|Nim}}== |
|||
{{trans|Kotlin}} |
|||
<syntaxhighlight lang="nim">import random, sequtils, strformat |
|||
proc equalBirthdays(nSharers, groupSize, nRepetitions: int): float = |
|||
randomize(1) |
|||
var eq = 0 |
|||
for _ in 1..nRepetitions: |
|||
var group: array[1..365, int] |
|||
for _ in 1..groupSize: |
|||
inc group[rand(1..group.len)] |
|||
eq += ord(group.anyIt(it >= nSharers)) |
|||
result = eq * 100 / nRepetitions |
|||
proc main() = |
|||
var groupEst = 2 |
|||
for sharers in 2..5: |
|||
# Coarse. |
|||
var groupSize = groupEst + 1 |
|||
while equalBirthdays(sharers, groupSize, 100) < 50: |
|||
inc groupSize |
|||
# Finer. |
|||
let inf = (groupSize.toFloat - (groupSize - groupEst) / 4).toInt() |
|||
for gs in inf..(groupSize+998): |
|||
let eq = equalBirthdays(sharers, groupSize, 250) |
|||
if eq > 50: |
|||
groupSize = gs |
|||
break |
|||
# Finest. |
|||
for gs in (groupSize-1)..(groupSize+998): |
|||
let eq = equalBirthdays(sharers, gs, 50_000) |
|||
if eq > 50: |
|||
groupEst = gs |
|||
echo &"{sharers} independent people in a group of {gs:3} ", |
|||
&"share a common birthday ({eq:4.1f}%)" |
|||
break |
|||
main()</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 independent people in a group of 23 share a common birthday (50.7%) |
|||
3 independent people in a group of 87 share a common birthday (50.0%) |
|||
4 independent people in a group of 187 share a common birthday (50.2%) |
|||
5 independent people in a group of 313 share a common birthday (50.2%)</pre> |
|||
=={{header|PARI/GP}}== |
|||
<syntaxhighlight lang="parigp">simulate(n)=my(v=vecsort(vector(n,i,random(365))),t,c=1); for(i=2,n,if(v[i]>v[i-1],t=max(t,c);c=1,c++)); t |
|||
find(n)=my(guess=365*n-342,t);while(1, t=sum(i=1,1e3,simulate(guess)>=n)/1e3; if(t>550, guess--); if(t<450, guess++); if(450<=t && t<=550, return(guess))) |
|||
find(2) |
|||
find(3) |
|||
find(4) |
|||
find(5)</syntaxhighlight> |
|||
=={{header|Perl}}== |
|||
{{trans|Raku}} |
|||
<syntaxhighlight lang="perl">use strict; |
|||
use warnings; |
|||
use List::AllUtils qw(max min uniqnum count_by any); |
|||
use Math::Random qw(random_uniform_integer); |
|||
sub simulation { |
|||
my($c) = shift; |
|||
my $max_trials = 1_000_000; |
|||
my $min_trials = 10_000; |
|||
my $n = int 47 * ($c-1.5)**1.5; # OEIS/A050256: 16 86 185 307 |
|||
my $N = min $max_trials, max $min_trials, 1000 * sqrt $n; |
|||
while (1) { |
|||
my $yes = 0; |
|||
for (1..$N) { |
|||
my %birthday_freq = count_by { $_ } random_uniform_integer($n, 1, 365); |
|||
$yes++ if any { $birthday_freq{$_} >= $c } keys %birthday_freq; |
|||
} |
|||
my $p = $yes/$N; |
|||
return($n, $p) if $p > 0.5; |
|||
$N = min $max_trials, max $min_trials, int 1000/(0.5-$p)**1.75; |
|||
$n++; |
|||
} |
|||
} |
|||
printf "$_ people in a group of %s share a common birthday. (%.4f)\n", simulation($_) for 2..5</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 people in a group of 23 share a common birthday. (0.5083) |
|||
3 people in a group of 88 share a common birthday. (0.5120) |
|||
4 people in a group of 187 share a common birthday. (0.5034) |
|||
5 people in a group of 313 share a common birthday. (0.5008)</pre> |
|||
=={{header|Phix}}== |
|||
{{trans|D}} |
|||
<!--<syntaxhighlight lang="phix">(phixonline)[very/too slow]--> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">nDays</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">365</span> |
|||
<span style="color: #000080;font-style:italic;">-- 5 sigma confidence. Conventionally people think 3 sigmas are |
|||
-- good enough, but for the case of 5 people sharing a birthday, |
|||
-- 3 sigmas actually sometimes gives a slightly wrong answer.</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">nSigmas</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">5.0</span><span style="color: #0000FF;">;</span> <span style="color: #000080;font-style:italic;">-- Change to 3 for smaller run time.</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">simulate1</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">nPeople</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000080;font-style:italic;">-- |
|||
-- Given n people, if m of them have same birthday in one run. |
|||
--</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">days</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nDays</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">nPeople</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">day</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">rand</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nDays</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">days</span><span style="color: #0000FF;">[</span><span style="color: #000000;">day</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">days</span><span style="color: #0000FF;">[</span><span style="color: #000000;">day</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">==</span> <span style="color: #000000;">nCollisions</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #004600;">true</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #004600;">false</span><span style="color: #0000FF;">;</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">prob</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">np</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">pThresh</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000080;font-style:italic;">-- |
|||
-- Decide if the probablity of n out of np people sharing a birthday |
|||
-- is above or below pThresh, with nSigmas sigmas confidence. |
|||
-- If pThresh is very low or hi, minimum runs need to be much higher. |
|||
--</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">;</span> <span style="color: #000080;font-style:italic;">-- Probablity and standard deviation.</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">nRuns</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">yes</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">;</span> |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">nRuns</span><span style="color: #0000FF;"><</span><span style="color: #000000;">10</span> <span style="color: #008080;">or</span> <span style="color: #0000FF;">(</span><span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">pThresh</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;"><</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">nSigmas</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">yes</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">simulate1</span><span style="color: #0000FF;">(</span><span style="color: #000000;">np</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">nRuns</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">yes</span><span style="color: #0000FF;">/</span><span style="color: #000000;">nRuns</span> |
|||
<span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span> <span style="color: #0000FF;">*</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">1</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">nRuns</span><span style="color: #0000FF;">);</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">findHalfChance</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000080;font-style:italic;">-- Bisect for truth.</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dev</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">mid</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">lo</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">hi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nDays</span> <span style="color: #0000FF;">*</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">nCollisions</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">;</span> |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">lo</span> <span style="color: #0000FF;"><</span> <span style="color: #000000;">mid</span> <span style="color: #008080;">or</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;"><</span> <span style="color: #000000;">0.5</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">mid</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">((</span><span style="color: #000000;">hi</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">lo</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dev</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">prob</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mid</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">p</span> <span style="color: #0000FF;"><</span> <span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000000;">lo</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mid</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">;</span> |
|||
<span style="color: #008080;">else</span> |
|||
<span style="color: #000000;">hi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mid</span><span style="color: #0000FF;">;</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">hi</span> <span style="color: #0000FF;"><</span> <span style="color: #000000;">lo</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">findHalfChance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- reset</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dev</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mid</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">6</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,</span><span style="color: #000000;">np</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">findHalfChance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%d collision: %d people, P = %g +/- %g\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">nCollisions</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">np</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
|||
<pre> |
|||
2 collision: 23 people, P = 0.520699 +/- 0.00688426 |
|||
3 collision: 88 people, P = 0.507159 +/- 0.00238534 |
|||
4 collision: 187 people, P = 0.504129 +/- 0.00137625 |
|||
5 collision: 313 people, P = 0.501219 +/- 0.000406284 |
|||
6 collision: 460 people, P = 0.502131 +/- 0.000710091 |
|||
</pre> |
|||
Output with nSigmas = 5.0: |
|||
<pre> |
|||
2 collision: 23 people, P = 0.507817 +/- 0.00156278 |
|||
3 collision: 88 people, P = 0.512042 +/- 0.00240772 |
|||
4 collision: 187 people, P = 0.502546 +/- 0.000509275 |
|||
5 collision: 313 people, P = 0.501218 +/- 0.000243516 |
|||
6 collision: 460 people, P = 0.502901 +/- 0.000580137 |
|||
</pre> |
|||
=={{header|PL/I}}== |
=={{header|PL/I}}== |
||
< |
<syntaxhighlight lang="pl/i">*process source attributes xref; |
||
bd: Proc Options(main); |
bd: Proc Options(main); |
||
/*-------------------------------------------------------------------- |
/*-------------------------------------------------------------------- |
||
Line 206: | Line 1,442: | ||
End; |
End; |
||
End; |
End; |
||
End;</ |
End;</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 248: | Line 1,484: | ||
461 493948 506052 50.605% <- 50.622% <- 50.626% <- |
461 493948 506052 50.605% <- 50.622% <- 50.626% <- |
||
462 489416 510584 51.058% <- 51.029% <- 51.055% <-</pre> |
462 489416 510584 51.058% <- 51.029% <- 51.055% <-</pre> |
||
extended to verify REXX results: |
|||
<pre> 1000000 samples. Percentage of groups with at least 7 matches |
|||
Group size |
|||
621 503758 496242 49.624% |
|||
622 500320 499680 49.968% |
|||
623 497047 502953 50.295% <- |
|||
624 493679 506321 50.632% <- |
|||
625 491240 508760 50.876% <- |
|||
1000000 samples. Percentage of groups with at least 8 matches |
|||
=={{header|REXX}}== |
|||
Group size |
|||
===Version 1=== |
|||
796 504764 495236 49.524% |
|||
<lang rexx>/*REXX programs solves "birthday problem" via random number simulations.*/ |
|||
797 502537 497463 49.746% |
|||
parse arg grps samp seed . /*get optional arguments from CL.*/ |
|||
798 499488 500512 50.051% <- |
|||
if grps=='' | grps==',' then grps=5 /*Not specified? Use the default*/ |
|||
799 496658 503342 50.334% <- |
|||
if samp=='' | samp==',' then samp=100000 /*" " " " " */ |
|||
800 494773 505227 50.523% <- |
|||
if seed\==',' & seed \=='' then call random ,,seed /*repeatability? */ |
|||
diy =365 /*or: diy=365.25*/ /*the number of Days In a Year. */ |
|||
diyM=diy*100 /*this expands the RANDOM range. */ |
|||
/* [↓] get a rough estimate for %*/ |
|||
do g=2 to grps; s=0 /*perform through 2──►group size.*/ |
|||
do samp; @.=0 /*perform some number of trials. */ |
|||
do j=1 /*do until G dup birthdays found.*/ |
|||
day=random(1,diyM) % 100 /*expand random number generation*/ |
|||
@.day=@.day+1 /*record the # common birthdays.*/ |
|||
if @.day==g then leave /*'nuff G hits have occurred ··· */ |
|||
end /*j*/ |
|||
s=s+j /*add # of birthday hits to sum. */ |
|||
end /*samp*/ |
|||
1000000 samples. Percentage of groups with at least 9 matches |
|||
start.g=s/samp%1-g /*define where the try-outs start*/ |
|||
Group size |
|||
end /*g*/ |
|||
983 502613 497387 49.739% |
|||
___=' ' /*padding for easier eyeballing. */ |
|||
984 501665 498335 49.834% |
|||
say ___ ___ ___ 'sample size is ' samp /*show sample size of this run. */ |
|||
985 498606 501394 50.139% <- |
|||
say |
|||
986 497453 502547 50.255% <- |
|||
say ___ 'required' ___ 'group' ___ '% with required' |
|||
987 493816 506184 50.618% <- |
|||
say ___ ' common ' ___ ' size' ___ 'common birthdays' |
|||
say ___ '────────' ___ '─────' ___ '────────────────' |
|||
/* [↓] where the try-outs happen.*/ |
|||
do g=2 to grps /*perform through 2──►group size.*/ |
|||
do try=start.g; s=0 /*perform try-outs until avg>50%.*/ |
|||
do samp; @.=0 /*perform some number of trials. */ |
|||
do j=1 for try /*do until K dup birthdays found.*/ |
|||
day=random(1,diyM) % 100 /*expand random number generation*/ |
|||
@.day=@.day+1 /*record the # common birthdays.*/ |
|||
if @.day\==g then iterate /*not 'nuff G hits occurred? ··· */ |
|||
s=s+1 /*another common birthday found. */ |
|||
leave /* ··· and stop looking for more.*/ |
|||
end /*j*/ |
|||
end /*samp*/ |
|||
1000000 samples. Percentage of groups with at least10 matches |
|||
if s/samp>.5 then leave /*if the average is > 50%, stop. */ |
|||
Group size |
|||
end /*try*/ |
|||
1179 502910 497090 49.709% |
|||
1180 500906 499094 49.909% |
|||
1181 499079 500921 50.092% <- |
|||
1182 496957 503043 50.304% <- |
|||
1183 494414 505586 50.559% <-</pre> |
|||
=={{header|Python}}== |
|||
Note: the first (unused), version of function equal_birthdays() uses a different but equally valid interpretation of the phrase "common birthday". |
|||
<syntaxhighlight lang="python"> |
|||
from random import randint |
|||
def equal_birthdays(sharers=2, groupsize=23, rep=100000): |
|||
say ___ center(g,8) ___ right(try,5) ___ center((s/samp*100)'%',16) |
|||
'Note: 4 sharing common birthday may have 2 dates shared between two people each' |
|||
end /*g*/ |
|||
g = range(groupsize) |
|||
/*stick a fork in it, we're done.*/</lang> |
|||
sh = sharers - 1 |
|||
'''output''' |
|||
eq = sum((groupsize - len(set(randint(1,365) for i in g)) >= sh) |
|||
<pre style="overflow:scroll"> |
|||
for j in range(rep)) |
|||
return (eq * 100.) / rep |
|||
def equal_birthdays(sharers=2, groupsize=23, rep=100000): |
|||
required group % with required |
|||
common |
'Note: 4 sharing common birthday must all share same common day' |
||
g = range(groupsize) |
|||
──────── ───── ──────────────── |
|||
sh = sharers - 1 |
|||
2 23 50.75100% |
|||
eq = 0 |
|||
3 88 51.14700% |
|||
for j in range(rep): |
|||
4 187 50.30400% |
|||
group = [randint(1,365) for i in g] |
|||
5 313 50.14100% |
|||
if (groupsize - len(set(group)) >= sh and |
|||
any( group.count(member) >= sharers for member in set(group))): |
|||
eq += 1 |
|||
return (eq * 100.) / rep |
|||
group_est = [2] |
|||
for sharers in (2, 3, 4, 5): |
|||
groupsize = group_est[-1]+1 |
|||
while equal_birthdays(sharers, groupsize, 100) < 50.: |
|||
# Coarse |
|||
groupsize += 1 |
|||
for groupsize in range(int(groupsize - (groupsize - group_est[-1])/4.), groupsize + 999): |
|||
# Finer |
|||
eq = equal_birthdays(sharers, groupsize, 250) |
|||
if eq > 50.: |
|||
break |
|||
for groupsize in range(groupsize - 1, groupsize +999): |
|||
# Finest |
|||
eq = equal_birthdays(sharers, groupsize, 50000) |
|||
if eq > 50.: |
|||
break |
|||
group_est.append(groupsize) |
|||
print("%i independent people in a group of %s share a common birthday. (%5.1f)" % (sharers, groupsize, eq))</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 independent people in a group of 23 share a common birthday. ( 50.9) |
|||
3 independent people in a group of 87 share a common birthday. ( 50.0) |
|||
4 independent people in a group of 188 share a common birthday. ( 50.9) |
|||
5 independent people in a group of 314 share a common birthday. ( 50.6)</pre> |
|||
===Enumeration method=== |
|||
The following enumerates all birthday distributation of n people in a year. It's patentedly unscalable. |
|||
<syntaxhighlight lang="python">from collections import defaultdict |
|||
days = 365 |
|||
def find_half(c): |
|||
# inc_people takes birthday combinations of n people and generates the |
|||
# new set for n+1 |
|||
def inc_people(din, over): |
|||
# 'over' is the number of combinations that have at least c people |
|||
# sharing a birthday. These are not contained in the set. |
|||
dout,over = defaultdict(int), over * days |
|||
for k,s in din.items(): |
|||
for i,v in enumerate(k): |
|||
if v + 1 >= c: |
|||
over += s |
|||
else: |
|||
dout[tuple(sorted(k[0:i] + (v + 1,) + k[i+1:]))] += s |
|||
dout[(1,) + k] += s * (days - len(k)) |
|||
return dout, over |
|||
d, combos, good, n = {():1}, 1, 0, 0 |
|||
# increase number of people until at least half of the cases have at |
|||
# at least c people sharing a birthday |
|||
while True: |
|||
n += 1 |
|||
combos *= days # or, combos = sum(d.values()) + good |
|||
d,good = inc_people(d, good) |
|||
#!!! print d.items() |
|||
if good * 2 >= combos: |
|||
return n, good, combos |
|||
# In all fairness, I don't know if the code works for x >= 4: I probably don't |
|||
# have enough RAM for it, and certainly not enough patience. But it should. |
|||
# In theory. |
|||
for x in range(2, 5): |
|||
n, good, combos = find_half(x) |
|||
print "%d of %d people sharing birthday: %d out of %d combos"% (x, n, good, combos) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 of 23 people sharing birthday: 43450860051057961364418604769486195435604861663267741453125 out of 85651679353150321236814267844395152689354622364044189453125 combos |
|||
3 of 88 people sharing birthday: 1549702400401473425983277424737696914087385196361193892581987189461901608374448849589919219974092878625057027641693544686424625999709818279964664633586995549680467629183956971001416481439048256933422687688148710727691650390625 out of 3032299345394764867793392128292779133654078653518318790345269064871742118915665927782934165016667902517875712171754287171746462419635313222013443107339730598579399174951673950890087953259632858049599235528148710727691650390625 combos |
|||
...? |
|||
</pre> |
</pre> |
||
=== |
===Enumeration method #2=== |
||
<syntaxhighlight lang="python"># ought to use a memoize class for all this |
|||
<lang rexx> /*-------------------------------------------------------------------- |
|||
# factorial |
|||
def fact(n, cache={0:1}): |
|||
if not n in cache: |
|||
cache[n] = n * fact(n - 1) |
|||
return cache[n] |
|||
# permutations |
|||
def perm(n, k, cache={}): |
|||
if not (n,k) in cache: |
|||
cache[(n,k)] = fact(n) / fact(n - k) |
|||
return cache[(n,k)] |
|||
def choose(n, k, cache={}): |
|||
if not (n,k) in cache: |
|||
cache[(n,k)] = perm(n, k) / fact(k) |
|||
return cache[(n, k)] |
|||
# ways of distribute p people's birthdays into d days, with |
|||
# no more than m sharing any one day |
|||
def combos(d, p, m, cache={}): |
|||
if not p: return 1 |
|||
if not m: return 0 |
|||
if p <= m: return d**p # any combo would satisfy |
|||
k = (d, p, m) |
|||
if not k in cache: |
|||
result = 0 |
|||
for x in range(0, p//m + 1): |
|||
c = combos(d - x, p - x * m, m - 1) |
|||
# ways to occupy x days with m people each |
|||
if c: result += c * choose(d, x) * perm(p, x * m) / fact(m)**x |
|||
cache[k] = result |
|||
return cache[k] |
|||
def find_half(m): |
|||
n = 0 |
|||
while True: |
|||
n += 1 |
|||
total = 365 ** n |
|||
c = total - combos(365, n, m - 1) |
|||
if c * 2 >= total: |
|||
print "%d of %d people: %d/%d combos" % (n, m, c, total) |
|||
return |
|||
for x in range(2, 6): find_half(x)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
23 of 2 people: 43450860....3125/85651679....3125 combos |
|||
88 of 3 people: 15497...50390625/30322...50390625 combos |
|||
187 of 4 people: 708046698...0703125/1408528546...0703125 combos |
|||
313 of 5 people: 498385488882289...2578125/99464149835930...2578125 combos |
|||
</pre> |
|||
=={{header|Racket}}== |
|||
{{trans|Python}}Based on the Python task. For three digits precision use 250000 repetitions. For four digits precision use 25000000 repetitions, but it’s very slow. See discussion page. |
|||
<syntaxhighlight lang="racket">#lang racket |
|||
#;(define repetitions 25000000) ; for \sigma=1/10000 |
|||
(define repetitions 250000) ; for \sigma=1/1000 |
|||
(define coarse-repetitions 2500) |
|||
(define (vector-inc! v pos) |
|||
(vector-set! v pos (add1 (vector-ref v pos)))) |
|||
(define (equal-birthdays sharers group-size repetitions) |
|||
(/ (for/sum ([j (in-range repetitions)]) |
|||
(let ([days (make-vector 365 0)]) |
|||
(for ([person (in-range group-size)]) |
|||
(vector-inc! days (random 365))) |
|||
(if (>= (apply max (vector->list days)) sharers) |
|||
1 0))) |
|||
repetitions)) |
|||
(define (search-coarse-group-size sharers) |
|||
(let loop ([coarse-group-size 2]) |
|||
(let ([coarse-probability |
|||
(equal-birthdays sharers coarse-group-size coarse-repetitions)]) |
|||
(if (> coarse-probability .5) |
|||
coarse-group-size |
|||
(loop (add1 coarse-group-size)))))) |
|||
(define (search-upwards sharers group-size) |
|||
(let ([probability (equal-birthdays sharers group-size repetitions)]) |
|||
(if (> probability .5) |
|||
(values group-size probability) |
|||
(search-upwards sharers (add1 group-size))))) |
|||
(define (search-downwards sharers group-size last-probability) |
|||
(let ([probability (equal-birthdays sharers group-size repetitions)]) |
|||
(if (> probability .5) |
|||
(search-downwards sharers (sub1 group-size) probability) |
|||
(values (add1 group-size) last-probability)))) |
|||
(define (search-from sharers group-size) |
|||
(let ([probability (equal-birthdays sharers group-size repetitions)]) |
|||
(if (> probability .5) |
|||
(search-downwards sharers (sub1 group-size) probability) |
|||
(search-upwards sharers (add1 group-size))))) |
|||
(for ([sharers (in-range 2 6)]) |
|||
(let-values ([(group-size probability) |
|||
(search-from sharers (search-coarse-group-size sharers))]) |
|||
(printf "~a independent people in a group of ~a share a common birthday. (~a%)\n" |
|||
sharers group-size (~r (* probability 100) #:precision '(= 2)))))</syntaxhighlight> |
|||
'''Output''' |
|||
<pre>2 independent people in a group of 23 share a common birthday. (50.80%) |
|||
3 independent people in a group of 88 share a common birthday. (51.19%) |
|||
4 independent people in a group of 187 share a common birthday. (50.18%) |
|||
5 independent people in a group of 313 share a common birthday. (50.17%) |
|||
</pre> |
|||
=={{header|Raku}}== |
|||
(formerly Perl 6) |
|||
Gives correct answers, but more of a proof-of-concept at this point, even with max-trials at 250K it is too slow to be practical. |
|||
<syntaxhighlight lang="raku" line>sub simulation ($c) { |
|||
my $max-trials = 250_000; |
|||
my $min-trials = 5_000; |
|||
my $n = floor 47 * ($c-1.5)**1.5; # OEIS/A050256: 16 86 185 307 |
|||
my $N = min $max-trials, max $min-trials, 1000 * sqrt $n; |
|||
loop { |
|||
my $p = $N R/ elems grep { .elems > 0 }, ((grep { $_>=$c }, values bag (^365).roll($n)) xx $N); |
|||
return($n, $p) if $p > 0.5; |
|||
$N = min $max-trials, max $min-trials, floor 1000/(0.5-$p); |
|||
$n++; |
|||
} |
|||
} |
|||
printf "$_ people in a group of %s share a common birthday. (%.3f)\n", simulation($_) for 2..5;</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 people in a group of 23 share a common birthday. (0.506) |
|||
3 people in a group of 88 share a common birthday. (0.511) |
|||
4 people in a group of 187 share a common birthday. (0.500) |
|||
5 people in a group of 313 share a common birthday. (0.507)</pre> |
|||
=={{header|REXX}}== |
|||
===version 1=== |
|||
The root finding method used is to find the average number of people to share a birthday, and then use the '''floor''' |
|||
<br>of that value (less the group size) as a starting point to find a new group size with an expected size that exceeds |
|||
<br>50% duplicate birthdays of the required size. |
|||
This REXX version doesn't need a precalculated group size to find the percentage required to exceed 50%. |
|||
<syntaxhighlight lang="rexx">/*REXX pgm examines the birthday problem via random# simulation (with specifiable parms)*/ |
|||
parse arg dups samp seed . /*get optional arguments from the CL. */ |
|||
if dups=='' | dups=="," then dups= 10 /*Not specified? Then use the default.*/ |
|||
if samp=='' | samp=="," then samp= 10000 /* " " " " " " */ |
|||
if datatype(seed, 'W') then call random ,,seed /*RANDOM seed given for repeatability ?*/ |
|||
diy = 365 /*alternative: diy=365.25*/ /*the number of Days In a Year. */ |
|||
diyM= diy*100 /*this expands the RANDOM (BIF) range.*/ |
|||
do g=2 to dups; s= 0 /*perform through 2 ──► duplicate size*/ |
|||
do samp; @.= 0 /*perform some number of trials. */ |
|||
do j=0 until @.day==g /*perform until G dup. birthdays found.*/ |
|||
day= random(1, diyM) % 100 /*expand range RANDOM number generation*/ |
|||
@.day= @.day + 1 /*record the number of common birthdays*/ |
|||
end /*j*/ /* [↓] adjust for the DO loop index.*/ |
|||
s= s+j /*add number of birthday hits to sum. */ |
|||
end /*samp*/ /* [↓] % 1 rounds down the division.*/ |
|||
start.g= s/samp % 1 - g /*define where the try─outs start. */ |
|||
end /*g*/ /* [↑] get a rough estimate for %. */ |
|||
say right('sample size is ' samp, 40); say /*display this run's sample size. */ |
|||
say ' required trial % with required' |
|||
say ' duplicates size common birthdays' |
|||
say ' ──────────── ─────── ──────────────────' |
|||
do g=2 to dups /*perform through 2 ──► duplicate size*/ |
|||
do try=start.g until s/samp>=.5; s= 0 /* " try─outs until average ≥ 50%.*/ |
|||
do samp; @.= 0 /* " some number of trials. */ |
|||
do try; day= random(1, diyM) % 100 /* " until G dup. birthdays found.*/ |
|||
@.day= @.day + 1 /*record the number of common birthdays*/ |
|||
if @.day==g then do; s=s+1; leave; end /*found enough G (birthday) hits ? */ |
|||
end /*try;*/ |
|||
end /*samp*/ |
|||
end /*try=start.g*/ /* [↑] where the try─outs happen. */ |
|||
say right(g, 15) right(try, 15) center( format( s / samp * 100, , 4)'%', 30) |
|||
end /*g*/ /*stick a fork in it, we're all done. */</syntaxhighlight> |
|||
{{out|output|text= when using the default inputs:}} |
|||
<pre> |
|||
sample size is 10000 |
|||
required trial % with required |
|||
duplicates size common birthdays |
|||
──────────── ─────── ────────────────── |
|||
2 23 50.2300% |
|||
3 87 50.2400% |
|||
4 187 50.3800% |
|||
5 312 50.0100% |
|||
6 458 50.5200% |
|||
7 622 50.3900% |
|||
8 798 50.1700% |
|||
9 984 50.5700% |
|||
10 1182 51.4000% |
|||
</pre> |
|||
===version 2=== |
|||
<syntaxhighlight lang="rexx"> /*-------------------------------------------------------------------- |
|||
* 04.11.2013 Walter Pachl translated from PL/I |
* 04.11.2013 Walter Pachl translated from PL/I |
||
* Take samp samples of groups with gs persons and check |
* Take samp samples of groups with gs persons and check |
||
Line 340: | Line 1,838: | ||
Say format(gs,10) cnt.0 cnt.1 100*hits||'%'||arrow |
Say format(gs,10) cnt.0 cnt.1 100*hits||'%'||arrow |
||
End |
End |
||
End</ |
End</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 382: | Line 1,880: | ||
461 49316 50684 50.68400% <- |
461 49316 50684 50.68400% <- |
||
462 49121 50879 50.87900% <-</pre> |
462 49121 50879 50.87900% <-</pre> |
||
=={{header|Ruby}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="ruby">def equalBirthdays(nSharers, groupSize, nRepetitions) |
|||
eq = 0 |
|||
for i in 1 .. nRepetitions |
|||
group = [0] * 365 |
|||
for j in 1 .. groupSize |
|||
group[rand(group.length)] += 1 |
|||
end |
|||
eq += group.any? { |n| n >= nSharers } ? 1 : 0 |
|||
end |
|||
return (eq * 100.0) / nRepetitions |
|||
end |
|||
def main |
|||
groupEst = 2 |
|||
for sharers in 2 .. 5 |
|||
# Coarse |
|||
groupSize = groupEst + 1 |
|||
while equalBirthdays(sharers, groupSize, 100) < 50.0 |
|||
groupSize += 1 |
|||
end |
|||
# Finer |
|||
inf = (groupSize - (groupSize - groupEst) / 4.0).floor |
|||
for gs in inf .. groupSize + 999 |
|||
eq = equalBirthdays(sharers, groupSize, 250) |
|||
if eq > 50.0 then |
|||
groupSize = gs |
|||
break |
|||
end |
|||
end |
|||
# Finest |
|||
for gs in groupSize - 1 .. groupSize + 999 |
|||
eq = equalBirthdays(sharers, gs, 50000) |
|||
if eq > 50.0 then |
|||
groupEst = gs |
|||
print "%d independant people in a group of %s share a common birthday. (%5.1f)\n" % [sharers, gs, eq] |
|||
break |
|||
end |
|||
end |
|||
end |
|||
end |
|||
main()</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 independant people in a group of 23 share a common birthday. ( 50.5) |
|||
3 independant people in a group of 90 share a common birthday. ( 53.3) |
|||
4 independant people in a group of 187 share a common birthday. ( 50.3) |
|||
5 independant people in a group of 313 share a common birthday. ( 50.3)</pre> |
|||
=={{header|Scala}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="Scala"> |
|||
import scala.util.Random |
|||
import scala.util.control.Breaks._ |
|||
object Test { |
|||
def equalBirthdays( |
|||
nSharers: Int, |
|||
groupSize: Int, |
|||
nRepetitions: Int |
|||
): Double = { |
|||
val rand = new Random(1) |
|||
var eq = 0 |
|||
for (_ <- 1 to nRepetitions) { |
|||
val group = new Array[Int](365) |
|||
for (_ <- 1 to groupSize) { |
|||
val birthday = rand.nextInt(group.length) |
|||
group(birthday) += 1 |
|||
} |
|||
if (group.exists(_ >= nSharers)) eq += 1 |
|||
} |
|||
(eq * 100.0) / nRepetitions |
|||
} |
|||
def main(args: Array[String]): Unit = { |
|||
var groupEst = 2 |
|||
for (sharers <- 2 until 6) { |
|||
// Coarse. |
|||
var groupSize = groupEst + 1 |
|||
while (equalBirthdays(sharers, groupSize, 100) < 50.0) |
|||
groupSize += 1 |
|||
// Finer. |
|||
val inf = (groupSize - (groupSize - groupEst) / 4.0).toInt |
|||
breakable({ |
|||
for (gs <- inf until groupSize + 999) { |
|||
val eq = equalBirthdays(sharers, groupSize, 250) |
|||
if (eq > 50.0) { |
|||
groupSize = gs |
|||
break |
|||
} |
|||
} |
|||
}) |
|||
// Finest. |
|||
breakable({ |
|||
for (gs <- (groupSize - 1) until groupSize + 999) { |
|||
val eq = equalBirthdays(sharers, gs, 50000) |
|||
if (eq > 50.0) { |
|||
groupEst = gs |
|||
println( |
|||
f"$sharers independent people in a group of $gs share a common birthday. ($eq%5.1f)" |
|||
) |
|||
break |
|||
} |
|||
} |
|||
}) |
|||
} |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 independent people in a group of 23 share a common birthday. ( 50.6) |
|||
3 independent people in a group of 87 share a common birthday. ( 50.4) |
|||
4 independent people in a group of 187 share a common birthday. ( 50.1) |
|||
5 independent people in a group of 314 share a common birthday. ( 50.2) |
|||
</pre> |
|||
=={{header|SQL}}== |
|||
birthday.sql |
|||
<syntaxhighlight lang="sql"> |
|||
with |
|||
c as |
|||
( |
|||
select |
|||
500 nrep, |
|||
50 maxgsiz |
|||
from dual |
|||
), |
|||
reps as |
|||
( |
|||
select level rep |
|||
from dual |
|||
cross join c |
|||
connect by level <= c.nrep |
|||
), |
|||
pers as |
|||
( |
|||
select |
|||
round(sqrt(2*level)) npers |
|||
from dual |
|||
cross join c |
|||
connect by level <= c.maxgsiz*(c.maxgsiz+1)/2 |
|||
), |
|||
bds as |
|||
( |
|||
select |
|||
reps.rep, |
|||
pers.npers, |
|||
floor(dbms_random.value(1,366)) bd |
|||
from |
|||
reps |
|||
cross join pers |
|||
), |
|||
mtch as |
|||
( |
|||
select |
|||
bds.npers, |
|||
case count(distinct bds.bd ) when bds.npers then 0 else 1 end match |
|||
from bds |
|||
group by |
|||
bds.rep, |
|||
bds.npers, |
|||
null |
|||
order by |
|||
bds.npers |
|||
), |
|||
nm as |
|||
( |
|||
select mtch.npers, sum (mtch.match) nmatch |
|||
from mtch |
|||
group by mtch.npers |
|||
), |
|||
sol as |
|||
( |
|||
select first_value ( nm.npers ) over ( order by abs ( nm.nmatch - c.nrep / 2 ) ) npers |
|||
from |
|||
nm |
|||
cross join c |
|||
) |
|||
select npers |
|||
from sol where rownum = 1 |
|||
; |
|||
</syntaxhighlight> |
|||
SQL> @ birthday.sql |
|||
Connected. |
|||
NPERS |
|||
---------- |
|||
23 |
|||
=={{header|Tcl}}== |
|||
<syntaxhighlight lang="tcl">proc birthdays {num {same 2}} { |
|||
for {set i 0} {$i < $num} {incr i} { |
|||
set b [expr {int(rand() * 365)}] |
|||
if {[incr bs($b)] >= $same} { |
|||
return 1 |
|||
} |
|||
} |
|||
return 0 |
|||
} |
|||
proc estimateBirthdayChance {num same} { |
|||
# Gives a reasonably close estimate with minimal execution time; the idea |
|||
# is to keep the amount that one random value may influence the result |
|||
# fairly constant. |
|||
set count [expr {$num * 100 / $same}] |
|||
set x 0 |
|||
for {set i 0} {$i < $count} {incr i} { |
|||
incr x [birthdays $num $same] |
|||
} |
|||
return [expr {double($x) / $count}] |
|||
} |
|||
foreach {count from to} {2 20 25 3 85 90 4 183 190 5 310 315} { |
|||
puts "identifying level for $count people with same birthday" |
|||
for {set i $from} {$i <= $to} {incr i} { |
|||
set chance [estimateBirthdayChance $i $count] |
|||
puts [format "%d people => %%%.2f chance of %d people with same birthday" \ |
|||
$i [expr {$chance * 100}] $count] |
|||
if {$chance >= 0.5} { |
|||
puts "level found: $i people" |
|||
break |
|||
} |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
identifying level for 2 people with same birthday |
|||
20 people => %43.40 chance of 2 people with same birthday |
|||
21 people => %44.00 chance of 2 people with same birthday |
|||
22 people => %46.91 chance of 2 people with same birthday |
|||
23 people => %53.48 chance of 2 people with same birthday |
|||
level found: 23 people |
|||
identifying level for 3 people with same birthday |
|||
85 people => %47.97 chance of 3 people with same birthday |
|||
86 people => %48.46 chance of 3 people with same birthday |
|||
87 people => %49.55 chance of 3 people with same birthday |
|||
88 people => %50.66 chance of 3 people with same birthday |
|||
level found: 88 people |
|||
identifying level for 4 people with same birthday |
|||
183 people => %48.02 chance of 4 people with same birthday |
|||
184 people => %47.67 chance of 4 people with same birthday |
|||
185 people => %48.89 chance of 4 people with same birthday |
|||
186 people => %49.98 chance of 4 people with same birthday |
|||
187 people => %50.99 chance of 4 people with same birthday |
|||
level found: 187 people |
|||
identifying level for 5 people with same birthday |
|||
310 people => %48.52 chance of 5 people with same birthday |
|||
311 people => %48.14 chance of 5 people with same birthday |
|||
312 people => %49.07 chance of 5 people with same birthday |
|||
313 people => %49.63 chance of 5 people with same birthday |
|||
314 people => %49.59 chance of 5 people with same birthday |
|||
315 people => %51.79 chance of 5 people with same birthday |
|||
level found: 315 people |
|||
</pre> |
|||
=={{header|Wren}}== |
|||
{{trans|Kotlin}} |
|||
{{libheader|Wren-fmt}} |
|||
<syntaxhighlight lang="wren">import "random" for Random |
|||
import "./fmt" for Fmt |
|||
var equalBirthdays = Fn.new { |nSharers, groupSize, nRepetitions| |
|||
var rand = Random.new(12345) |
|||
var eq = 0 |
|||
for (i in 0...nRepetitions) { |
|||
var group = List.filled(365, 0) |
|||
for (j in 0...groupSize) { |
|||
var r = rand.int(group.count) |
|||
group[r] = group[r] + 1 |
|||
} |
|||
eq = eq + (group.any { |i| i >= nSharers } ? 1 : 0) |
|||
} |
|||
return eq * 100 / nRepetitions |
|||
} |
|||
var groupEst = 2 |
|||
for (sharers in 2..5) { |
|||
// Coarse |
|||
var groupSize = groupEst + 1 |
|||
while (equalBirthdays.call(sharers, groupSize, 100) < 50) groupSize = groupSize + 1 |
|||
// Finer |
|||
var inf = (groupSize - (groupSize - groupEst) / 4).floor |
|||
for (gs in inf...groupSize + 999) { |
|||
var eq = equalBirthdays.call(sharers, groupSize, 250) |
|||
if (eq > 50) { |
|||
groupSize = gs |
|||
break |
|||
} |
|||
} |
|||
// Finest |
|||
for (gs in groupSize - 1...groupSize + 999) { |
|||
var eq = equalBirthdays.call(sharers, gs, 50000) |
|||
if (eq > 50) { |
|||
groupEst = gs |
|||
Fmt.write("$d independent people in a group of $3d ", sharers, gs) |
|||
Fmt.print("share a common birthday $2.1f\%", eq) |
|||
break |
|||
} |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 independent people in a group of 23 share a common birthday 51.0% |
|||
3 independent people in a group of 88 share a common birthday 51.0% |
|||
4 independent people in a group of 187 share a common birthday 50.1% |
|||
5 independent people in a group of 314 share a common birthday 50.7% |
|||
</pre> |
|||
=={{header|XPL0}}== |
|||
<syntaxhighlight lang="xpl0">func Sim(N); |
|||
\Simulate birthdays and return number of people to have N same days |
|||
int N, I, People, R; |
|||
char Days(365); |
|||
[for I:= 0 to 365-1 do Days(I):= 0; |
|||
People:= 0; |
|||
loop [R:= Ran(365); |
|||
Days(R):= Days(R)+1; |
|||
People:= People+1; |
|||
if Days(R) = N then return People; |
|||
]; |
|||
]; |
|||
int N, Sum, Trial; |
|||
[for N:= 2 to 5 do |
|||
[Sum:= 0; |
|||
for Trial:= 1 to 10000 do |
|||
Sum:= Sum + Sim(N); |
|||
IntOut(0, N); Text(0, ": "); IntOut(0, Sum/10000); CrLf(0); |
|||
]; |
|||
]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2: 24 |
|||
3: 88 |
|||
4: 187 |
|||
5: 311 |
|||
</pre> |
|||
=={{header|zkl}}== |
|||
Pure simulation; adding a person to a population until there are the required number of collisions, then repeating that a bunch of times to get an average. |
|||
<syntaxhighlight lang="zkl">fcn bdays(N){ // N is shared birthdays in a population |
|||
year:=(0).pump(365,List.createLong(365).write,0); // 365 days == one year |
|||
shared:=people:=0; do{ // add a new person to population |
|||
bday:=(0).random(365); // with this birthday [0..364] |
|||
shared=shared.max(year[bday]+=1); people+=1; |
|||
}while(shared<N); |
|||
people // size of simulated population that contains N shared birthdays |
|||
} |
|||
fcn simulate(N,T){ avg:=0.0; do(T){ avg+=bdays(N) } avg/=T; } // N shared, T trials |
|||
foreach n in ([1..5]){ |
|||
println("Average of %d people in a populatation of %s share birthdays" |
|||
.fmt(n,simulate(n,0d10_000))); |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Average of 1 people in a populatation of 1 share birthdays |
|||
Average of 2 people in a populatation of 24.7199 share birthdays |
|||
Average of 3 people in a populatation of 88.6416 share birthdays |
|||
Average of 4 people in a populatation of 186.849 share birthdays |
|||
Average of 5 people in a populatation of 312.399 share birthdays |
|||
</pre> |
Latest revision as of 11:58, 20 January 2024
This page uses content from Wikipedia. The current wikipedia article is at Birthday Problem. The original RosettaCode article was extracted from the wikipedia article № 296054030 of 21:44, 12 June 2009 . The list of authors can be seen in the page history. As with Rosetta Code, the pre 5 June 2009 text of Wikipedia is available under the GNU FDL. (See links for details on variance) |
In probability theory, the birthday problem, or birthday paradox This is not a paradox in the sense of leading to a logical contradiction, but is called a paradox because the mathematical truth contradicts naïve intuition: most people estimate that the chance is much lower than 50%. pertains to the probability that in a set of randomly chosen people some pair of them will have the same birthday. In a group of at least 23 randomly chosen people, there is more than 50% probability that some pair of them will both have been born on the same day. For 57 or more people, the probability is more than 99%, and it reaches 100% when the number of people reaches 366 (by the pigeon hole principle, ignoring leap years). The mathematics behind this problem leads to a well-known cryptographic attack called the birthday attack.
- Task
Using simulation, estimate the number of independent people required in a group before we can expect a better than even chance that at least 2 independent people in a group share a common birthday. Furthermore: Simulate and thus estimate when we can expect a better than even chance that at least 3, 4 & 5 independent people of the group share a common birthday. For simplicity assume that all of the people are alive...
- Suggestions for improvement
- Estimating the error in the estimate to help ensure the estimate is accurate to 4 decimal places.
- Converging to the th solution using a root finding method, as opposed to using an extensive search.
- Kudos (κῦδος) for finding the solution by proof (in a programming language) rather than by construction and simulation.
- See also
- Wolfram entry: Birthday Problem
11l
F equal_birthdays(sharers, groupsize, rep)
V eq = 0
L 0 .< rep
V group = [0] * 365
L 0 .< groupsize
group[random:(group.len)]++
I any(group.map(c -> c >= @sharers))
eq++
R (eq * 100.) / rep
V group_est = 2
L(sharers) 2..5
V groupsize = group_est + 1
L equal_birthdays(sharers, groupsize, 100) < 50.
groupsize++
L(gs) Int(groupsize - (groupsize - group_est) / 4.) .< groupsize + 999
V eq = equal_birthdays(sharers, gs, 250)
I eq > 50.
groupsize = gs
L.break
L(gs) groupsize - 1 .< groupsize + 999
V eq = equal_birthdays(sharers, gs, 50'000)
I eq > 50.
group_est = gs
print(‘#. independent people in a group of #. share a common birthday. (#3.1)’.format(sharers, gs, eq))
L.break
- Output:
2 independent people in a group of 23 share a common birthday. ( 50.6) 3 independent people in a group of 87 share a common birthday. ( 50.1) 4 independent people in a group of 187 share a common birthday. ( 50.5) 5 independent people in a group of 313 share a common birthday. ( 50.4)
Ada
This solution assumes a 4-year cycle, with three 365-day years and one leap year.
with Ada.Command_Line, Ada.Text_IO, Ada.Numerics.Discrete_random;
procedure Birthday_Test is
Samples: constant Positive := Integer'Value(Ada.Command_Line.Argument(1));
-- our experiment: Generate a X (birth-)days and check for Y-collisions
-- the constant "Samples" is the number of repetitions of this experiment
subtype Day is integer range 0 .. 365; -- this includes leap_days
subtype Extended_Day is Integer range 0 .. 365*4; -- a four-year cycle
package ANDR is new Ada.Numerics.Discrete_Random(Extended_Day);
Random_Generator: ANDR.Generator;
function Random_Day return Day is (ANDR.Random(Random_Generator) / 4);
-- days 0 .. 364 are equally probable, leap-day 365 is 4* less probable
type Checkpoint is record
Multiplicity: Positive;
Person_Count: Positive;
end record;
Checkpoints: constant array(Positive range <>) of Checkpoint
:= ( (2, 22), (2, 23), (3, 86), (3, 87), (3, 88),
(4, 186), (4, 187), (5, 312), (5, 313), (5, 314) );
type Result_Type is array(Checkpoints'Range) of Natural;
Result: Result_Type := (others => 0);
-- how often is a 2-collision in a group of 22 or 23, ..., a 5-collision
-- in a group of 312 .. 314
procedure Experiment(Result: in out Result_Type) is
-- run the experiment once!
A_Year: array(Day) of Natural := (others => 0);
A_Day: Day;
Multiplicity: Natural := 0;
People: Positive := 1;
begin
for I in Checkpoints'Range loop
while People <= Checkpoints(I).Person_Count loop
A_Day := Random_Day;
A_Year(A_Day) := A_Year(A_Day)+1;
if A_Year(A_Day) > Multiplicity then
Multiplicity := Multiplicity + 1;
end if;
People := People + 1;
end loop;
if Multiplicity >= Checkpoints(I).Multiplicity then
Result(I) := Result(I) + 1;
-- found a Multipl.-collision in a group of Person_Cnt.
end if;
end loop;
end Experiment;
package TIO renames Ada.Text_IO;
package FIO is new TIO.Float_IO(Float);
begin
-- initialize the random generator
ANDR.Reset(Random_Generator);
-- repeat the experiment Samples times
for I in 1 .. Samples loop
Experiment(Result);
end loop;
-- print the results
TIO.Put_Line("Birthday-Test with" & Integer'Image(Samples) & " samples:");
for I in Result'Range loop
FIO.Put(Float(Result(I))/Float(Samples), Fore => 3, Aft => 6, Exp => 0);
TIO.Put_Line
("% of groups with" & Integer'Image(Checkpoints(I).Person_Count) &
" have" & Integer'Image(Checkpoints(I).Multiplicity) &
" persons sharing a common birthday.");
end loop;
end Birthday_Test;
- Output:
Running the program with a sample size 500_000_000 took about 25 minutes on a slow pc.
./birthday_test 500_000_000 Birthday-Test with 500000000 samples: 0.475292% of groups with 22 have 2 persons sharing a common birthday. 0.506882% of groups with 23 have 2 persons sharing a common birthday. 0.487155% of groups with 86 have 3 persons sharing a common birthday. 0.498788% of groups with 87 have 3 persons sharing a common birthday. 0.510391% of groups with 88 have 3 persons sharing a common birthday. 0.494970% of groups with 186 have 4 persons sharing a common birthday. 0.501825% of groups with 187 have 4 persons sharing a common birthday. 0.495137% of groups with 312 have 5 persons sharing a common birthday. 0.500010% of groups with 313 have 5 persons sharing a common birthday. 0.504888% of groups with 314 have 5 persons sharing a common birthday.
An interesting observation: The probability for groups of 313 persons having 5 persons sharing a common birthday is almost exactly 0.5. Note that a solution based on 365-day years, i.e., a solution ignoring leap days, would generate slightly but significantly larger probabilities.
ALGOL 68
File: Birthday_problem.a68
#!/usr/bin/a68g --script #
# -*- coding: utf-8 -*- #
REAL desired probability := 0.5; # 50% #
REAL upb year = 365 + 1/4 # - 3/400 but alive, ignore those born prior to 1901 #,
INT upb sample size = 100 000,
upb common = 5 ;
FORMAT name int fmt = $g": "g(-0)"; "$,
name real fmt = $g": "g(-0,4)"; "$,
name percent fmt = $g": "g(-0,2)"%; "$;
printf((
name real fmt,
"upb year",upb year,
name int fmt,
"upb common",upb common,
"upb sample size",upb sample size,
$l$
));
INT required common := 1; # initial value #
FOR group size FROM required common WHILE required common <= upb common DO
INT sample with no required common := 0;
TO upb sample size DO
# generate sample #
[group size]INT sample;
FOR i TO UPB sample DO sample[i] := ENTIER(random * upb year) + 1 OD;
FOR birthday i TO UPB sample DO
INT birthday = sample[birthday i];
INT number in common := 1;
# special case = 1 #
IF number in common >= required common THEN
found required common
FI;
FOR birthday j FROM birthday i + 1 TO UPB sample DO
IF birthday = sample[birthday j] THEN
number in common +:= 1;
IF number in common >= required common THEN
found required common
FI
FI
OD
OD # days in year #;
sample with no required common +:= 1;
found required common: SKIP
OD # sample size #;
REAL portion of years with required common birthdays =
(upb sample size - sample with no required common) / upb sample size;
print(".");
IF portion of years with required common birthdays > desired probability THEN
printf((
$l$,
name int fmt,
"required common",required common,
"group size",group size,
# "sample with no required common",sample with no required common, #
name percent fmt,
"%age of years with required common birthdays",portion of years with required common birthdays*100,
$l$
));
required common +:= 1
FI
OD # group size #
Output:
upb year: 365.2500; upb common: 5; upb sample size: 100000; . required common: 1; group size: 1; %age of years with required common birthdays: 100.00%; ...................... required common: 2; group size: 23; %age of years with required common birthdays: 50.71%; ................................................................. required common: 3; group size: 88; %age of years with required common birthdays: 50.90%; ................................................................................................... required common: 4; group size: 187; %age of years with required common birthdays: 50.25%; ............................................................................................................................... required common: 5; group size: 314; %age of years with required common birthdays: 50.66%;
C
Computing probabilities to 5 sigmas of confidence. It's very slow, chiefly because to make sure a probability like 0.5006 is indeed above .5 instead of just statistical fluctuation, you have to run the simulation millions of times.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define DEBUG 0 // set this to 2 for a lot of numbers on output
#define DAYS 365
#define EXCESS (RAND_MAX / DAYS * DAYS)
int days[DAYS];
inline int rand_day(void)
{
int n;
while ((n = rand()) >= EXCESS);
return n / (EXCESS / DAYS);
}
// given p people, if n of them have same birthday in one run
int simulate1(int p, int n)
{
memset(days, 0, sizeof(days));
while (p--)
if (++days[rand_day()] == n) return 1;
return 0;
}
// decide if the probability of n out of np people sharing a birthday
// is above or below p_thresh, with n_sigmas sigmas confidence
// note that if p_thresh is very low or hi, minimum runs need to be much higher
double prob(int np, int n, double n_sigmas, double p_thresh, double *std_dev)
{
double p, d; // prob and std dev
int runs = 0, yes = 0;
do {
yes += simulate1(np, n);
p = (double) yes / ++runs;
d = sqrt(p * (1 - p) / runs);
if (DEBUG > 1)
printf("\t\t%d: %d %d %g %g \r", np, yes, runs, p, d);
} while (runs < 10 || fabs(p - p_thresh) < n_sigmas * d);
if (DEBUG > 1) putchar('\n');
*std_dev = d;
return p;
}
// bisect for truth
int find_half_chance(int n, double *p, double *dev)
{
int lo, hi, mid;
reset:
lo = 0;
hi = DAYS * (n - 1) + 1;
do {
mid = (hi + lo) / 2;
// 5 sigma confidence. Conventionally people think 3 sigmas are good
// enough, but for case of 5 people sharing birthday, 3 sigmas actually
// sometimes give a slightly wrong answer
*p = prob(mid, n, 5, .5, dev);
if (DEBUG)
printf("\t%d %d %d %g %g\n", lo, mid, hi, *p, *dev);
if (*p < .5) lo = mid + 1;
else hi = mid;
if (hi < lo) {
// this happens when previous precisions were too low;
// easiest fix: reset
if (DEBUG) puts("\tMade a mess, will redo.");
goto reset;
}
} while (lo < mid || *p < .5);
return mid;
}
int main(void)
{
int n, np;
double p, d;
srand(time(0));
for (n = 2; n <= 5; n++) {
np = find_half_chance(n, &p, &d);
printf("%d collision: %d people, P = %g +/- %g\n",
n, np, p, d);
}
return 0;
}
- Output:
2 collision: 23 people, P = 0.508741 +/- 0.00174794 3 collision: 88 people, P = 0.509034 +/- 0.00180628 4 collision: 187 people, P = 0.501812 +/- 0.000362394 5 collision: 313 people, P = 0.500641 +/- 0.000128174
C++
#include <iostream>
#include <random>
#include <vector>
double equalBirthdays(int nSharers, int groupSize, int nRepetitions) {
std::default_random_engine generator;
std::uniform_int_distribution<int> distribution(0, 364);
std::vector<int> group(365);
int eq = 0;
for (int i = 0; i < nRepetitions; i++) {
std::fill(group.begin(), group.end(), 0);
for (int j = 0; j < groupSize; j++) {
int day = distribution(generator);
group[day]++;
}
if (std::any_of(group.cbegin(), group.cend(), [nSharers](int c) { return c >= nSharers; })) {
eq++;
}
}
return (100.0 * eq) / nRepetitions;
}
int main() {
int groupEst = 2;
for (int sharers = 2; sharers < 6; sharers++) {
// Coarse
int groupSize = groupEst + 1;
while (equalBirthdays(sharers, groupSize, 100) < 50.0) {
groupSize++;
}
// Finer
int inf = (int)(groupSize - (groupSize - groupEst) / 4.0f);
for (int gs = inf; gs < groupSize + 999; gs++) {
double eq = equalBirthdays(sharers, groupSize, 250);
if (eq > 50.0) {
groupSize = gs;
break;
}
}
// Finest
for (int gs = groupSize - 1; gs < groupSize + 999; gs++) {
double eq = equalBirthdays(sharers, gs, 50000);
if (eq > 50.0) {
groupEst = gs;
printf("%d independant people in a group of %d share a common birthday. (%5.1f)\n", sharers, gs, eq);
break;
}
}
}
return 0;
}
- Output:
2 independant people in a group of 23 share a common birthday. ( 50.6) 3 independant people in a group of 87 share a common birthday. ( 50.2) 4 independant people in a group of 186 share a common birthday. ( 50.0) 5 independant people in a group of 313 share a common birthday. ( 50.5)
D
import std.stdio, std.random, std.algorithm, std.conv;
/// For sharing common birthday must all share same common day.
double equalBirthdays(in uint nSharers, in uint groupSize,
in uint nRepetitions, ref Xorshift rng) {
uint eq = 0;
foreach (immutable _; 0 .. nRepetitions) {
uint[365] group;
foreach (immutable __; 0 .. groupSize)
group[uniform(0, $, rng)]++;
eq += group[].any!(c => c >= nSharers);
}
return (eq * 100.0) / nRepetitions;
}
void main() {
auto rng = 1.Xorshift; // Fixed seed.
auto groupEst = 2;
foreach (immutable sharers; 2 .. 6) {
// Coarse.
auto groupSize = groupEst + 1;
while (equalBirthdays(sharers, groupSize, 100, rng) < 50.0)
groupSize++;
// Finer.
immutable inf = to!int(groupSize - (groupSize - groupEst) / 4.0);
foreach (immutable gs; inf .. groupSize + 999) {
immutable eq = equalBirthdays(sharers, groupSize, 250, rng);
if (eq > 50.0) {
groupSize = gs;
break;
}
}
// Finest.
foreach (immutable gs; groupSize - 1 .. groupSize + 999) {
immutable eq = equalBirthdays(sharers, gs, 50_000, rng);
if (eq > 50.0) {
groupEst = gs;
writefln("%d independent people in a group of %s share a common birthday. (%5.1f)",
sharers, gs, eq);
break;
}
}
}
}
- Output:
2 independent people in a group of 23 share a common birthday. ( 50.5) 3 independent people in a group of 87 share a common birthday. ( 50.1) 4 independent people in a group of 187 share a common birthday. ( 50.2) 5 independent people in a group of 313 share a common birthday. ( 50.3)
Run-time about 10.4 seconds with ldc2 compiler.
Alternative version:
import std.stdio, std.random, std.math;
enum nDays = 365;
// 5 sigma confidence. Conventionally people think 3 sigmas are good
// enough, but for case of 5 people sharing birthday, 3 sigmas
// actually sometimes give a slightly wrong answer.
enum double nSigmas = 3.0; // Currently 3 for smaller run time.
/// Given n people, if m of them have same birthday in one run.
bool simulate1(in uint nPeople, in uint nCollisions, ref Xorshift rng)
/*nothrow*/ @safe /*@nogc*/ {
static uint[nDays] days;
days[] = 0;
foreach (immutable _; 0 .. nPeople) {
immutable day = uniform(0, days.length, rng);
days[day]++;
if (days[day] == nCollisions)
return true;
}
return false;
}
/** Decide if the probablity of n out of np people sharing a birthday
is above or below pThresh, with nSigmas sigmas confidence.
If pThresh is very low or hi, minimum runs need to be much higher. */
double prob(in uint np, in uint nCollisions, in double pThresh,
out double stdDev, ref Xorshift rng) {
double p, d; // Probablity and standard deviation.
uint nRuns = 0, yes = 0;
do {
yes += simulate1(np, nCollisions, rng);
nRuns++;
p = double(yes) / nRuns;
d = sqrt(p * (1 - p) / nRuns);
debug if (yes % 50_000 == 0)
printf("\t\t%d: %d %d %g %g \r", np, yes, nRuns, p, d);
} while (nRuns < 10 || abs(p - pThresh) < (nSigmas * d));
debug '\n'.putchar;
stdDev = d;
return p;
}
/// Bisect for truth.
uint findHalfChance(in uint nCollisions, out double p, out double dev, ref Xorshift rng) {
uint mid;
RESET:
uint lo = 0;
uint hi = nDays * (nCollisions - 1) + 1;
do {
mid = (hi + lo) / 2;
p = prob(mid, nCollisions, 0.5, dev, rng);
debug printf("\t%d %d %d %g %g\n", lo, mid, hi, p, dev);
if (p < 0.5)
lo = mid + 1;
else
hi = mid;
if (hi < lo) {
// This happens when previous precisions were too low;
// easiest fix: reset.
debug "\tMade a mess, will redo.".puts;
goto RESET;
}
} while (lo < mid || p < 0.5);
return mid;
}
void main() {
auto rng = Xorshift(unpredictableSeed);
foreach (immutable uint nCollisions; 2 .. 6) {
double p, d;
immutable np = findHalfChance(nCollisions, p, d, rng);
writefln("%d collision: %d people, P = %g +/- %g", nCollisions, np, p, d);
}
}
- Output:
2 collision: 23 people, P = 0.521934 +/- 0.00728933 3 collision: 88 people, P = 0.512367 +/- 0.00411469 4 collision: 187 people, P = 0.506974 +/- 0.00232306 5 collision: 313 people, P = 0.501588 +/- 0.000529277
Output with nSigmas = 5.0:
2 collision: 23 people, P = 0.508607 +/- 0.00172133 3 collision: 88 people, P = 0.511945 +/- 0.00238885 4 collision: 187 people, P = 0.503229 +/- 0.000645587 5 collision: 313 people, P = 0.501105 +/- 0.000221016
Delphi
program Birthday_problem;
{$APPTYPE CONSOLE}
uses
System.SysUtils;
const
_DAYS = 365;
var
days: array[0..(_DAYS) - 1] of Integer;
runs: Integer;
function rand_day: Integer; inline;
begin
Result := random(_DAYS);
end;
/// <summary>
/// given p people, if n of them have same birthday in one run
/// </summary>
function simulate1(p, n: Integer): Integer;
var
index, i: Integer;
begin
for i := 0 to High(days) do
begin
days[i] := 0;
end;
for i := 0 to p - 1 do
begin
index := rand_day();
inc(days[index]);
if days[index] = n then
Exit(1);
end;
Exit(0);
end;
/// <summary>
/// decide if the probability of n out of np people sharing a birthday
/// is above or below p_thresh, with n_sigmas sigmas confidence
/// note that if p_thresh is very low or hi, minimum runs need to be much higher
/// </summary>
function prob(np, n: Integer; n_sigmas, p_thresh: Double; var d: Double): Double;
var
p: Double;
runs, yes: Integer;
begin
runs := 0;
yes := 0;
repeat
yes := yes + (simulate1(np, n));
inc(runs);
p := yes / runs;
d := sqrt(p * (1 - p) / runs);
until ((runs >= 10) and (abs(p - p_thresh) >= (n_sigmas * d)));
Exit(p);
end;
/// <summary>
/// bisect for truth
/// </summary>
function find_half_chance(n: Integer; var p: Double; var d: Double): Integer;
var
lo, hi, mid: Integer;
label
reset;
begin
reset:
lo := 0;
hi := _DAYS * (n - 1) + 1;
repeat
mid := (hi + lo) div 2;
p := prob(mid, n, 3, 0.5, d);
if p < 0.5 then
lo := mid + 1
else
hi := mid;
if hi < lo then
goto reset;
until ((lo >= mid) and (p >= 0.5));
Exit(mid);
end;
var
n, np: Integer;
p, d: Double;
begin
Randomize;
writeln('Wait for calculate');
for n := 2 to 5 do
begin
np := find_half_chance(n, p, d);
writeln(format('%d collision: %d people, P = %.8f +/- %.8f', [n, np, p, d]));
end;
writeln('Press enter to exit');
readln;
end.
- Output:
Wait for calculate 2 collision: 23 people, P = 0,50411600 +/- 0,00137151 3 collision: 88 people, P = 0,51060651 +/- 0,00352751 4 collision: 188 people, P = 0,50856493 +/- 0,00285022 5 collision: 313 people, P = 0,50093354 +/- 0,00031116 Press enter to exit
FreeBASIC
Function Simulacion(N As Integer) As Integer
Dim As Integer i, dias(365)
For i = 0 To 365-1
dias(i) = 0
Next i
Dim As Integer R, personas = 0
Do
R = Rnd * 365
dias(R) += 1
personas += 1
If dias(R) = N Then Return personas
Loop
End Function
Dim As Integer N, grupo, t
For N = 2 To 5
grupo = 0
For t = 1 To 10000
grupo += Simulacion(N)
Next t
Print Using "Average of # people in a population of ### share birthdays"; N; Int(grupo/10000)
Next N
Sleep
Go
package main
import (
"fmt"
"math"
"math/rand"
"time"
)
const (
DEBUG = 0
DAYS = 365
n_sigmas = 5.
WORKERS = 16 // concurrent worker processes
RUNS = 1000 // runs per flight
)
func simulate1(p, n int, r *rand.Rand) int {
var days [DAYS]int
for i := 0; i < p; i++ {
days[r.Intn(DAYS)]++
}
for _, d := range days {
if d >= n {
return 1
}
}
return 0
}
// send yes's per fixed number of simulate1 runs until canceled
func work(p, n int, ych chan int, cancel chan bool) {
r := rand.New(rand.NewSource(time.Now().Unix() + rand.Int63()))
for {
select {
case <-cancel:
return
default:
}
y := 0
for i := 0; i < RUNS; i++ {
y += simulate1(p, n, r)
}
ych <- y
}
}
func prob(np, n int) (p, d float64) {
ych := make(chan int, WORKERS)
cancel := make(chan bool)
for i := 0; i < WORKERS; i++ {
go work(np, n, ych, cancel)
}
var runs, yes int
for {
yes += <-ych
runs += RUNS
fr := float64(runs)
p = float64(yes) / fr
d = math.Sqrt(p * (1 - p) / fr)
if DEBUG > 1 {
fmt.Println("\t\t", np, yes, runs, p, d)
}
// .5 here is the "even chance" threshold
if !(math.Abs(p-.5) < n_sigmas*d) {
close(cancel)
break
}
}
if DEBUG > 1 {
fmt.Println()
}
return
}
func find_half_chance(n int) (mid int, p, dev float64) {
reset:
lo := 0
hi := DAYS*(n-1) + 1
for {
mid = (hi + lo) / 2
p, dev = prob(mid, n)
if DEBUG > 0 {
fmt.Println("\t", lo, mid, hi, p, dev)
}
if p < .5 {
lo = mid + 1
} else {
hi = mid
}
if hi < lo {
if DEBUG > 0 {
fmt.Println("\tMade a mess, will redo.")
}
goto reset
}
if !(lo < mid || p < .5) {
break
}
}
return
}
func main() {
for n := 2; n <= 5; n++ {
np, p, d := find_half_chance(n)
fmt.Printf("%d collision: %d people, P = %.4f ± %.4f\n",
n, np, p, d)
}
}
2 collision: 23 people, P = 0.5081 ± 0.0016 3 collision: 88 people, P = 0.5155 ± 0.0029 4 collision: 187 people, P = 0.5041 ± 0.0008 5 collision: 313 people, P = 0.5015 ± 0.0003
Also based on the C version:
package main
import (
"fmt"
"math"
"math/rand"
"runtime"
"time"
)
type ProbeRes struct {
np int
p, d float64
}
type Frac struct {
n int
d int
}
var DaysInYear int = 365
func main() {
sigma := 5.0
for i := 2; i <= 5; i++ {
res := GetNP(i, sigma, 0.5)
fmt.Printf("%d collision: %d people, P = %.4f ± %.4f\n",
i, res.np, res.p, res.d)
}
}
func GetNP(n int, n_sigmas, p_thresh float64) (res ProbeRes) {
res.np = DaysInYear * (n - 1)
for i := 0; i < DaysInYear*(n-1); i++ {
tmp := probe(i, n, n_sigmas, p_thresh)
if tmp.p > p_thresh && tmp.np < res.np {
res = tmp
}
}
return
}
var numCPU = runtime.NumCPU()
func probe(np, n int, n_sigmas, p_thresh float64) ProbeRes {
var p, d float64
var runs, yes int
cRes := make(chan Frac, numCPU)
for i := 0; i < numCPU; i++ {
go SimN(np, n, 25, cRes)
}
for math.Abs(p-p_thresh) < n_sigmas*d || runs < 100 {
f := <-cRes
yes += f.n
runs += f.d
p = float64(yes) / float64(runs)
d = math.Sqrt(p * (1 - p) / float64(runs))
go SimN(np, n, runs/3, cRes)
}
return ProbeRes{np, p, d}
}
func SimN(np, n, ssize int, c chan Frac) {
r := rand.New(rand.NewSource(time.Now().UnixNano() + rand.Int63()))
yes := 0
for i := 0; i < ssize; i++ {
if Sim(np, n, r) {
yes++
}
}
c <- Frac{yes, ssize}
}
func Sim(p, n int, r *rand.Rand) (res bool) {
Cal := make([]int, DaysInYear)
for i := 0; i < p; i++ {
Cal[r.Intn(DaysInYear)]++
}
for _, v := range Cal {
if v >= n {
res = true
}
}
return
}
- Output:
2 collision: 23 people, P = 0.5068 ± 0.0013 3 collision: 88 people, P = 0.5148 ± 0.0028 4 collision: 187 people, P = 0.5020 ± 0.0004 5 collision: 313 people, P = 0.5011 ± 0.0002
Hy
We use a simple but not very accurate simulation method.
(import
[numpy :as np]
[random [randint]])
(defmacro incf (place)
`(+= ~place 1))
(defn birthday [required &optional [reps 20000] [ndays 365]]
(setv days (np.zeros (, reps ndays) np.int_))
(setv qualifying-reps (np.zeros reps np.bool_))
(setv group-size 1)
(setv count 0)
(while True
;(print group-size)
(for [r (range reps)]
(unless (get qualifying-reps r)
(setv day (randint 0 (dec ndays)))
(incf (get days (, r day)))
(when (= (get days (, r day)) required)
(setv (get qualifying-reps r) True)
(incf count))))
(when (> (/ (float count) reps) .5)
(break))
(incf group-size))
group-size)
(print (birthday 2))
(print (birthday 3))
(print (birthday 4))
(print (birthday 5))
J
Quicky approach (use a population of 1e5 people to get a quick estimate and then refine against a population of 1e8 people):
PopSmall=: 1e5 ?@# 365
PopBig=: 1e8 ?@# 365
countShared=: [: >./ #/.~
avg=: +/ % #
probShared=: (1 :0)("0)
:
NB. y: shared birthday count
NB. m: population
NB. x: sample size
avg ,y <: (-x) countShared\ m
)
estGroupSz=: 3 :0
approx=. (PopSmall probShared&y i.365) I. 0.5
n=. approx-(2+y)
refine=. n+(PopBig probShared&y approx+i:2+y) I. 0.5
assert. (2+y) > |approx-refine
refine, refine PopBig probShared y
)
Task cases:
estGroupSz 2
23 0.507254
estGroupSz 3
88 0.510737
estGroupSz 4
187 0.502878
estGroupSz 5
313 0.500903
So, for example, we need a group of 88 to have at least a 50% chance of 3 people in the group having the same birthday in a year of 365 days. And, in that case, the simulated probability was 51.0737%
Java
import static java.util.Arrays.stream;
import java.util.Random;
public class Test {
static double equalBirthdays(int nSharers, int groupSize, int nRepetitions) {
Random rand = new Random(1);
int eq = 0;
for (int i = 0; i < nRepetitions; i++) {
int[] group = new int[365];
for (int j = 0; j < groupSize; j++)
group[rand.nextInt(group.length)]++;
eq += stream(group).anyMatch(c -> c >= nSharers) ? 1 : 0;
}
return (eq * 100.0) / nRepetitions;
}
public static void main(String[] a) {
int groupEst = 2;
for (int sharers = 2; sharers < 6; sharers++) {
// Coarse.
int groupSize = groupEst + 1;
while (equalBirthdays(sharers, groupSize, 100) < 50.0)
groupSize++;
// Finer.
int inf = (int) (groupSize - (groupSize - groupEst) / 4.0);
for (int gs = inf; gs < groupSize + 999; gs++) {
double eq = equalBirthdays(sharers, groupSize, 250);
if (eq > 50.0) {
groupSize = gs;
break;
}
}
// Finest.
for (int gs = groupSize - 1; gs < groupSize + 999; gs++) {
double eq = equalBirthdays(sharers, gs, 50_000);
if (eq > 50.0) {
groupEst = gs;
System.out.printf("%d independent people in a group of "
+ "%s share a common birthday. (%5.1f)%n",
sharers, gs, eq);
break;
}
}
}
}
}
2 independent people in a group of 23 share a common birthday. ( 50,6) 3 independent people in a group of 87 share a common birthday. ( 50,4) 4 independent people in a group of 187 share a common birthday. ( 50,1) 5 independent people in a group of 314 share a common birthday. ( 50,2)
Julia
function equalbirthdays(sharers::Int, groupsize::Int; nrep::Int = 10000)
eq = 0
for _ in 1:nrep
group = rand(1:365, groupsize)
grset = Set(group)
if groupsize - length(grset) ≥ sharers - 1 &&
any(count(x -> x == d, group) ≥ sharers for d in grset)
eq += 1
end
end
return eq / nrep
end
gsizes = [2]
for sh in (2, 3, 4, 5)
local gsize = gsizes[end]
local freq
# Coarse
while equalbirthdays(sh, gsize; nrep = 100) < .5
gsize += 1
end
# Finer
for gsize in trunc(Int, gsize - (gsize - gsizes[end]) / 4):(gsize + 999)
if equalbirthdays(sh, gsize; nrep = 250) > 0.5
break
end
end
# Finest
for gsize in (gsize - 1):(gsize + 999)
freq = equalbirthdays(sh, gsize; nrep = 50000)
if freq > 0.5
break
end
end
push!(gsizes, gsize)
@printf("%i independent people in a group of %s share a common birthday. (%5.3f)\n", sh, gsize, freq)
end
- Output:
2 independent people in a group of 23 share a common birthday. (0.506) 3 independent people in a group of 88 share a common birthday. (0.510) 4 independent people in a group of 187 share a common birthday. (0.500) 5 independent people in a group of 314 share a common birthday. (0.507)
Kotlin
// version 1.1.3
import java.util.Random
fun equalBirthdays(nSharers: Int, groupSize: Int, nRepetitions: Int): Double {
val rand = Random(1L)
var eq = 0
for (i in 0 until nRepetitions) {
val group = IntArray(365)
for (j in 0 until groupSize) {
group[rand.nextInt(group.size)]++
}
eq += if (group.any { it >= nSharers}) 1 else 0
}
return eq * 100.0 / nRepetitions
}
fun main(args: Array<String>) {
var groupEst = 2
for (sharers in 2..5) {
// Coarse
var groupSize = groupEst + 1
while (equalBirthdays(sharers, groupSize, 100) < 50.0) groupSize++
// Finer
val inf = (groupSize - (groupSize - groupEst) / 4.0).toInt()
for (gs in inf until groupSize + 999) {
val eq = equalBirthdays(sharers, groupSize, 250)
if (eq > 50.0) {
groupSize = gs
break
}
}
// Finest
for (gs in groupSize - 1 until groupSize + 999) {
val eq = equalBirthdays(sharers, gs, 50_000)
if (eq > 50.0) {
groupEst = gs
print("$sharers independent people in a group of ${"%3d".format(gs)} ")
println("share a common birthday (${"%2.1f%%".format(eq)})")
break
}
}
}
}
- Output:
Expect runtime of about 15 seconds on a modest laptop:
2 independent people in a group of 23 share a common birthday (50.6%) 3 independent people in a group of 87 share a common birthday (50.4%) 4 independent people in a group of 187 share a common birthday (50.1%) 5 independent people in a group of 314 share a common birthday (50.2%)
Lasso
if(sys_listunboundmethods !>> 'randomgen') => {
define randomgen(len::integer,max::integer)::array => {
#len <= 0 ? return
local(out = array)
loop(#len) => { #out->insert(math_random(#max,1)) }
return #out
}
}
if(sys_listunboundmethods !>> 'hasdupe') => {
define hasdupe(a::array,threshold::integer) => {
with i in #a do => {
#a->find(#i)->size > #threshold-1 ? return true
}
return false
}
}
local(threshold = 2)
local(qty = 22, probability = 0.00, samplesize = 10000)
while(#probability < 50.00) => {^
local(dupeqty = 0)
loop(#samplesize) => {
local(x = randomgen(#qty,365))
hasdupe(#x,#threshold) ? #dupeqty++
}
#probability = (#dupeqty / decimal(#samplesize)) * 100
'Threshold: '+#threshold+', qty: '+#qty+' - probability: '+#probability+'\r'
#qty += 1
^}
- Output:
Threshold: 2, qty: 22 - probability: 47.810000 Threshold: 2, qty: 23 - probability: 51.070000 Threshold: 3, qty: 86 - probability: 48.400000 Threshold: 3, qty: 87 - probability: 49.200000 Threshold: 3, qty: 88 - probability: 52.900000 Threshold: 4, qty: 184 - probability: 48.000000 Threshold: 4, qty: 185 - probability: 49.800000 Threshold: 4, qty: 186 - probability: 49.600000 Threshold: 4, qty: 187 - probability: 48.900000 Threshold: 4, qty: 188 - probability: 50.700000 Threshold: 5, qty: 308 - probability: 48.130000 Threshold: 5, qty: 309 - probability: 48.430000 Threshold: 5, qty: 310 - probability: 48.640000 Threshold: 5, qty: 311 - probability: 49.370000 Threshold: 5, qty: 312 - probability: 49.180000 Threshold: 5, qty: 313 - probability: 49.540000 Threshold: 5, qty: 314 - probability: 50.000000
Nim
import random, sequtils, strformat
proc equalBirthdays(nSharers, groupSize, nRepetitions: int): float =
randomize(1)
var eq = 0
for _ in 1..nRepetitions:
var group: array[1..365, int]
for _ in 1..groupSize:
inc group[rand(1..group.len)]
eq += ord(group.anyIt(it >= nSharers))
result = eq * 100 / nRepetitions
proc main() =
var groupEst = 2
for sharers in 2..5:
# Coarse.
var groupSize = groupEst + 1
while equalBirthdays(sharers, groupSize, 100) < 50:
inc groupSize
# Finer.
let inf = (groupSize.toFloat - (groupSize - groupEst) / 4).toInt()
for gs in inf..(groupSize+998):
let eq = equalBirthdays(sharers, groupSize, 250)
if eq > 50:
groupSize = gs
break
# Finest.
for gs in (groupSize-1)..(groupSize+998):
let eq = equalBirthdays(sharers, gs, 50_000)
if eq > 50:
groupEst = gs
echo &"{sharers} independent people in a group of {gs:3} ",
&"share a common birthday ({eq:4.1f}%)"
break
main()
- Output:
2 independent people in a group of 23 share a common birthday (50.7%) 3 independent people in a group of 87 share a common birthday (50.0%) 4 independent people in a group of 187 share a common birthday (50.2%) 5 independent people in a group of 313 share a common birthday (50.2%)
PARI/GP
simulate(n)=my(v=vecsort(vector(n,i,random(365))),t,c=1); for(i=2,n,if(v[i]>v[i-1],t=max(t,c);c=1,c++)); t
find(n)=my(guess=365*n-342,t);while(1, t=sum(i=1,1e3,simulate(guess)>=n)/1e3; if(t>550, guess--); if(t<450, guess++); if(450<=t && t<=550, return(guess)))
find(2)
find(3)
find(4)
find(5)
Perl
use strict;
use warnings;
use List::AllUtils qw(max min uniqnum count_by any);
use Math::Random qw(random_uniform_integer);
sub simulation {
my($c) = shift;
my $max_trials = 1_000_000;
my $min_trials = 10_000;
my $n = int 47 * ($c-1.5)**1.5; # OEIS/A050256: 16 86 185 307
my $N = min $max_trials, max $min_trials, 1000 * sqrt $n;
while (1) {
my $yes = 0;
for (1..$N) {
my %birthday_freq = count_by { $_ } random_uniform_integer($n, 1, 365);
$yes++ if any { $birthday_freq{$_} >= $c } keys %birthday_freq;
}
my $p = $yes/$N;
return($n, $p) if $p > 0.5;
$N = min $max_trials, max $min_trials, int 1000/(0.5-$p)**1.75;
$n++;
}
}
printf "$_ people in a group of %s share a common birthday. (%.4f)\n", simulation($_) for 2..5
- Output:
2 people in a group of 23 share a common birthday. (0.5083) 3 people in a group of 88 share a common birthday. (0.5120) 4 people in a group of 187 share a common birthday. (0.5034) 5 people in a group of 313 share a common birthday. (0.5008)
Phix
constant nDays = 365 -- 5 sigma confidence. Conventionally people think 3 sigmas are -- good enough, but for the case of 5 people sharing a birthday, -- 3 sigmas actually sometimes gives a slightly wrong answer. constant nSigmas = 5.0; -- Change to 3 for smaller run time. function simulate1(integer nPeople, nCollisions) -- -- Given n people, if m of them have same birthday in one run. -- sequence days = repeat(0,nDays) for p=1 to nPeople do integer day = rand(nDays) days[day] += 1 if days[day] == nCollisions then return true end if end for return false; end function function prob(integer np, nCollisions, atom pThresh) -- -- Decide if the probablity of n out of np people sharing a birthday -- is above or below pThresh, with nSigmas sigmas confidence. -- If pThresh is very low or hi, minimum runs need to be much higher. -- atom p, d; -- Probablity and standard deviation. integer nRuns = 0, yes = 0; while nRuns<10 or (abs(p - pThresh) < (nSigmas * d)) do yes += simulate1(np, nCollisions) nRuns += 1 p = yes/nRuns d = sqrt(p * (1 - p) / nRuns); end while return {p,d} end function function findHalfChance(integer nCollisions) -- Bisect for truth. atom p, dev integer mid = 1, lo = 0, hi = nDays * (nCollisions - 1) + 1; while lo < mid or p < 0.5 do mid = floor((hi + lo) / 2) {p,dev} = prob(mid, nCollisions, 0.5) if (p < 0.5) then lo = mid + 1; else hi = mid; end if if (hi < lo) then return findHalfChance(nCollisions) -- reset end if end while return {p,dev,mid} end function for nCollisions=2 to 6 do atom {p,d,np} = findHalfChance(nCollisions) printf(1,"%d collision: %d people, P = %g +/- %g\n", {nCollisions, np, p, d}) end for
- Output:
2 collision: 23 people, P = 0.520699 +/- 0.00688426 3 collision: 88 people, P = 0.507159 +/- 0.00238534 4 collision: 187 people, P = 0.504129 +/- 0.00137625 5 collision: 313 people, P = 0.501219 +/- 0.000406284 6 collision: 460 people, P = 0.502131 +/- 0.000710091
Output with nSigmas = 5.0:
2 collision: 23 people, P = 0.507817 +/- 0.00156278 3 collision: 88 people, P = 0.512042 +/- 0.00240772 4 collision: 187 people, P = 0.502546 +/- 0.000509275 5 collision: 313 people, P = 0.501218 +/- 0.000243516 6 collision: 460 people, P = 0.502901 +/- 0.000580137
PL/I
*process source attributes xref;
bd: Proc Options(main);
/*--------------------------------------------------------------------
* 04.11.2013 Walter Pachl
* Take samp samples of groups with gs persons and check
*how many of the groups have at least match persons with same birthday
*-------------------------------------------------------------------*/
Dcl (float,random) Builtin;
Dcl samp Bin Fixed(31) Init(1000000);
Dcl arr(0:366) Bin Fixed(31);
Dcl r Bin fixed(31);
Dcl i Bin fixed(31);
Dcl ok Bin fixed(31);
Dcl g Bin fixed(31);
Dcl gs Bin fixed(31);
Dcl match Bin fixed(31);
Dcl cnt(0:1) Bin Fixed(31);
Dcl lo(6) Bin Fixed(31) Init(0,21,85,185,311,458);
Dcl hi(6) Bin Fixed(31) Init(0,25,89,189,315,462);
Dcl rf Bin Float(63);
Dcl hits Bin Float(63);
Dcl arrow Char(3);
Do match=2 To 6;
Put Edit(' ')(Skip,a);
Put Edit(samp,' samples. Percentage of groups with at least',
match,' matches')(Skip,f(8),a,f(2),a);
Put Edit('Group size')(Skip,a);
Do gs=lo(match) To hi(match);
cnt=0;
Do i=1 To samp;
ok=0;
arr=0;
Do g=1 To gs;
rf=random();
r=rf*365+1;
arr(r)+=1;
If arr(r)=match Then Do;
/* Put Edit(r)(Skip,f(4));*/
ok=1;
End;
End;
cnt(ok)+=1;
End;
hits=float(cnt(1))/samp;
If hits>=.5 Then arrow=' <-';
Else arrow='';
Put Edit(gs,cnt(0),cnt(1),100*hits,'%',arrow)
(Skip,f(10),2(f(7)),f(8,3),a,a);
End;
End;
End;
Output:
1000000 samples. Percentage of groups with at least 2 matches Group size 3000000 500000 samples 21 556903 443097 44.310% 44.343% 44.347% 22 524741 475259 47.526% 47.549% 47.521% 23 492034 507966 50.797% <- 50.735% <- 50.722% <- 24 462172 537828 53.783% <- 53.815% <- 53.838% <- 25 431507 568493 56.849% <- 56.849% <- 56.842% <- 1000000 samples. Percentage of groups with at least 3 matches Group size 85 523287 476713 47.671% 47.638% 47.631% 86 512219 487781 48.778% 48.776% 48.821% 87 499874 500126 50.013% <- 49.902% 49.903% 88 488197 511803 51.180% <- 51.127% <- 51.096% <- 89 478044 521956 52.196% <- 52.263% <- 52.290% <- 1000000 samples. Percentage of groups with at least 4 matches Group size 185 511352 488648 48.865% 48.868% 48.921% 186 503888 496112 49.611% 49.601% 49.568% 187 497844 502156 50.216% <- 50.258% <- 50.297% <- 188 490490 509510 50.951% <- 50.916% <- 50.946% <- 189 482893 517107 51.711% <- 51.645% <- 51.655% <- 1000000 samples. Percentage of groups with at least 5 matches Group size 311 508743 491257 49.126% 49.158% 49.164% 312 503524 496476 49.648% 49.631% 49.596% 313 498244 501756 50.176% <- 50.139% <- 50.095% <- 314 494032 505968 50.597% <- 50.636% <- 50.586% <- 315 489821 510179 51.018% <- 51.107% <- 51.114% <- 1000000 samples. Percentage of groups with at least 6 matches Group size 458 505225 494775 49.478% 49.498% 49.512% 459 501871 498129 49.813% 49.893% 49.885% 460 497719 502281 50.228% <- 50.278% <- 50.248% <- 461 493948 506052 50.605% <- 50.622% <- 50.626% <- 462 489416 510584 51.058% <- 51.029% <- 51.055% <-
extended to verify REXX results:
1000000 samples. Percentage of groups with at least 7 matches Group size 621 503758 496242 49.624% 622 500320 499680 49.968% 623 497047 502953 50.295% <- 624 493679 506321 50.632% <- 625 491240 508760 50.876% <- 1000000 samples. Percentage of groups with at least 8 matches Group size 796 504764 495236 49.524% 797 502537 497463 49.746% 798 499488 500512 50.051% <- 799 496658 503342 50.334% <- 800 494773 505227 50.523% <- 1000000 samples. Percentage of groups with at least 9 matches Group size 983 502613 497387 49.739% 984 501665 498335 49.834% 985 498606 501394 50.139% <- 986 497453 502547 50.255% <- 987 493816 506184 50.618% <- 1000000 samples. Percentage of groups with at least10 matches Group size 1179 502910 497090 49.709% 1180 500906 499094 49.909% 1181 499079 500921 50.092% <- 1182 496957 503043 50.304% <- 1183 494414 505586 50.559% <-
Python
Note: the first (unused), version of function equal_birthdays() uses a different but equally valid interpretation of the phrase "common birthday".
from random import randint
def equal_birthdays(sharers=2, groupsize=23, rep=100000):
'Note: 4 sharing common birthday may have 2 dates shared between two people each'
g = range(groupsize)
sh = sharers - 1
eq = sum((groupsize - len(set(randint(1,365) for i in g)) >= sh)
for j in range(rep))
return (eq * 100.) / rep
def equal_birthdays(sharers=2, groupsize=23, rep=100000):
'Note: 4 sharing common birthday must all share same common day'
g = range(groupsize)
sh = sharers - 1
eq = 0
for j in range(rep):
group = [randint(1,365) for i in g]
if (groupsize - len(set(group)) >= sh and
any( group.count(member) >= sharers for member in set(group))):
eq += 1
return (eq * 100.) / rep
group_est = [2]
for sharers in (2, 3, 4, 5):
groupsize = group_est[-1]+1
while equal_birthdays(sharers, groupsize, 100) < 50.:
# Coarse
groupsize += 1
for groupsize in range(int(groupsize - (groupsize - group_est[-1])/4.), groupsize + 999):
# Finer
eq = equal_birthdays(sharers, groupsize, 250)
if eq > 50.:
break
for groupsize in range(groupsize - 1, groupsize +999):
# Finest
eq = equal_birthdays(sharers, groupsize, 50000)
if eq > 50.:
break
group_est.append(groupsize)
print("%i independent people in a group of %s share a common birthday. (%5.1f)" % (sharers, groupsize, eq))
- Output:
2 independent people in a group of 23 share a common birthday. ( 50.9) 3 independent people in a group of 87 share a common birthday. ( 50.0) 4 independent people in a group of 188 share a common birthday. ( 50.9) 5 independent people in a group of 314 share a common birthday. ( 50.6)
Enumeration method
The following enumerates all birthday distributation of n people in a year. It's patentedly unscalable.
from collections import defaultdict
days = 365
def find_half(c):
# inc_people takes birthday combinations of n people and generates the
# new set for n+1
def inc_people(din, over):
# 'over' is the number of combinations that have at least c people
# sharing a birthday. These are not contained in the set.
dout,over = defaultdict(int), over * days
for k,s in din.items():
for i,v in enumerate(k):
if v + 1 >= c:
over += s
else:
dout[tuple(sorted(k[0:i] + (v + 1,) + k[i+1:]))] += s
dout[(1,) + k] += s * (days - len(k))
return dout, over
d, combos, good, n = {():1}, 1, 0, 0
# increase number of people until at least half of the cases have at
# at least c people sharing a birthday
while True:
n += 1
combos *= days # or, combos = sum(d.values()) + good
d,good = inc_people(d, good)
#!!! print d.items()
if good * 2 >= combos:
return n, good, combos
# In all fairness, I don't know if the code works for x >= 4: I probably don't
# have enough RAM for it, and certainly not enough patience. But it should.
# In theory.
for x in range(2, 5):
n, good, combos = find_half(x)
print "%d of %d people sharing birthday: %d out of %d combos"% (x, n, good, combos)
- Output:
2 of 23 people sharing birthday: 43450860051057961364418604769486195435604861663267741453125 out of 85651679353150321236814267844395152689354622364044189453125 combos 3 of 88 people sharing birthday: 1549702400401473425983277424737696914087385196361193892581987189461901608374448849589919219974092878625057027641693544686424625999709818279964664633586995549680467629183956971001416481439048256933422687688148710727691650390625 out of 3032299345394764867793392128292779133654078653518318790345269064871742118915665927782934165016667902517875712171754287171746462419635313222013443107339730598579399174951673950890087953259632858049599235528148710727691650390625 combos ...?
Enumeration method #2
# ought to use a memoize class for all this
# factorial
def fact(n, cache={0:1}):
if not n in cache:
cache[n] = n * fact(n - 1)
return cache[n]
# permutations
def perm(n, k, cache={}):
if not (n,k) in cache:
cache[(n,k)] = fact(n) / fact(n - k)
return cache[(n,k)]
def choose(n, k, cache={}):
if not (n,k) in cache:
cache[(n,k)] = perm(n, k) / fact(k)
return cache[(n, k)]
# ways of distribute p people's birthdays into d days, with
# no more than m sharing any one day
def combos(d, p, m, cache={}):
if not p: return 1
if not m: return 0
if p <= m: return d**p # any combo would satisfy
k = (d, p, m)
if not k in cache:
result = 0
for x in range(0, p//m + 1):
c = combos(d - x, p - x * m, m - 1)
# ways to occupy x days with m people each
if c: result += c * choose(d, x) * perm(p, x * m) / fact(m)**x
cache[k] = result
return cache[k]
def find_half(m):
n = 0
while True:
n += 1
total = 365 ** n
c = total - combos(365, n, m - 1)
if c * 2 >= total:
print "%d of %d people: %d/%d combos" % (n, m, c, total)
return
for x in range(2, 6): find_half(x)
- Output:
23 of 2 people: 43450860....3125/85651679....3125 combos 88 of 3 people: 15497...50390625/30322...50390625 combos 187 of 4 people: 708046698...0703125/1408528546...0703125 combos 313 of 5 people: 498385488882289...2578125/99464149835930...2578125 combos
Racket
Based on the Python task. For three digits precision use 250000 repetitions. For four digits precision use 25000000 repetitions, but it’s very slow. See discussion page.
#lang racket
#;(define repetitions 25000000) ; for \sigma=1/10000
(define repetitions 250000) ; for \sigma=1/1000
(define coarse-repetitions 2500)
(define (vector-inc! v pos)
(vector-set! v pos (add1 (vector-ref v pos))))
(define (equal-birthdays sharers group-size repetitions)
(/ (for/sum ([j (in-range repetitions)])
(let ([days (make-vector 365 0)])
(for ([person (in-range group-size)])
(vector-inc! days (random 365)))
(if (>= (apply max (vector->list days)) sharers)
1 0)))
repetitions))
(define (search-coarse-group-size sharers)
(let loop ([coarse-group-size 2])
(let ([coarse-probability
(equal-birthdays sharers coarse-group-size coarse-repetitions)])
(if (> coarse-probability .5)
coarse-group-size
(loop (add1 coarse-group-size))))))
(define (search-upwards sharers group-size)
(let ([probability (equal-birthdays sharers group-size repetitions)])
(if (> probability .5)
(values group-size probability)
(search-upwards sharers (add1 group-size)))))
(define (search-downwards sharers group-size last-probability)
(let ([probability (equal-birthdays sharers group-size repetitions)])
(if (> probability .5)
(search-downwards sharers (sub1 group-size) probability)
(values (add1 group-size) last-probability))))
(define (search-from sharers group-size)
(let ([probability (equal-birthdays sharers group-size repetitions)])
(if (> probability .5)
(search-downwards sharers (sub1 group-size) probability)
(search-upwards sharers (add1 group-size)))))
(for ([sharers (in-range 2 6)])
(let-values ([(group-size probability)
(search-from sharers (search-coarse-group-size sharers))])
(printf "~a independent people in a group of ~a share a common birthday. (~a%)\n"
sharers group-size (~r (* probability 100) #:precision '(= 2)))))
Output
2 independent people in a group of 23 share a common birthday. (50.80%) 3 independent people in a group of 88 share a common birthday. (51.19%) 4 independent people in a group of 187 share a common birthday. (50.18%) 5 independent people in a group of 313 share a common birthday. (50.17%)
Raku
(formerly Perl 6) Gives correct answers, but more of a proof-of-concept at this point, even with max-trials at 250K it is too slow to be practical.
sub simulation ($c) {
my $max-trials = 250_000;
my $min-trials = 5_000;
my $n = floor 47 * ($c-1.5)**1.5; # OEIS/A050256: 16 86 185 307
my $N = min $max-trials, max $min-trials, 1000 * sqrt $n;
loop {
my $p = $N R/ elems grep { .elems > 0 }, ((grep { $_>=$c }, values bag (^365).roll($n)) xx $N);
return($n, $p) if $p > 0.5;
$N = min $max-trials, max $min-trials, floor 1000/(0.5-$p);
$n++;
}
}
printf "$_ people in a group of %s share a common birthday. (%.3f)\n", simulation($_) for 2..5;
- Output:
2 people in a group of 23 share a common birthday. (0.506) 3 people in a group of 88 share a common birthday. (0.511) 4 people in a group of 187 share a common birthday. (0.500) 5 people in a group of 313 share a common birthday. (0.507)
REXX
version 1
The root finding method used is to find the average number of people to share a birthday, and then use the floor
of that value (less the group size) as a starting point to find a new group size with an expected size that exceeds
50% duplicate birthdays of the required size.
This REXX version doesn't need a precalculated group size to find the percentage required to exceed 50%.
/*REXX pgm examines the birthday problem via random# simulation (with specifiable parms)*/
parse arg dups samp seed . /*get optional arguments from the CL. */
if dups=='' | dups=="," then dups= 10 /*Not specified? Then use the default.*/
if samp=='' | samp=="," then samp= 10000 /* " " " " " " */
if datatype(seed, 'W') then call random ,,seed /*RANDOM seed given for repeatability ?*/
diy = 365 /*alternative: diy=365.25*/ /*the number of Days In a Year. */
diyM= diy*100 /*this expands the RANDOM (BIF) range.*/
do g=2 to dups; s= 0 /*perform through 2 ──► duplicate size*/
do samp; @.= 0 /*perform some number of trials. */
do j=0 until @.day==g /*perform until G dup. birthdays found.*/
day= random(1, diyM) % 100 /*expand range RANDOM number generation*/
@.day= @.day + 1 /*record the number of common birthdays*/
end /*j*/ /* [↓] adjust for the DO loop index.*/
s= s+j /*add number of birthday hits to sum. */
end /*samp*/ /* [↓] % 1 rounds down the division.*/
start.g= s/samp % 1 - g /*define where the try─outs start. */
end /*g*/ /* [↑] get a rough estimate for %. */
say right('sample size is ' samp, 40); say /*display this run's sample size. */
say ' required trial % with required'
say ' duplicates size common birthdays'
say ' ──────────── ─────── ──────────────────'
do g=2 to dups /*perform through 2 ──► duplicate size*/
do try=start.g until s/samp>=.5; s= 0 /* " try─outs until average ≥ 50%.*/
do samp; @.= 0 /* " some number of trials. */
do try; day= random(1, diyM) % 100 /* " until G dup. birthdays found.*/
@.day= @.day + 1 /*record the number of common birthdays*/
if @.day==g then do; s=s+1; leave; end /*found enough G (birthday) hits ? */
end /*try;*/
end /*samp*/
end /*try=start.g*/ /* [↑] where the try─outs happen. */
say right(g, 15) right(try, 15) center( format( s / samp * 100, , 4)'%', 30)
end /*g*/ /*stick a fork in it, we're all done. */
- output when using the default inputs:
sample size is 10000 required trial % with required duplicates size common birthdays ──────────── ─────── ────────────────── 2 23 50.2300% 3 87 50.2400% 4 187 50.3800% 5 312 50.0100% 6 458 50.5200% 7 622 50.3900% 8 798 50.1700% 9 984 50.5700% 10 1182 51.4000%
version 2
/*--------------------------------------------------------------------
* 04.11.2013 Walter Pachl translated from PL/I
* Take samp samples of groups with gs persons and check
*how many of the groups have at least match persons with same birthday
*-------------------------------------------------------------------*/
samp=100000
lo='0 21 85 185 311 458'
hi='0 25 89 189 315 462'
Do match=2 To 6
Say ' '
Say samp' samples . Percentage of groups with at least',
match ' matches'
Say 'Group size'
Do gs=word(lo,match) To word(hi,match)
cnt.=0
Do i=1 To samp
ok=0
arr.=0
Do g=1 To gs
r=random(1,365)
arr.r=arr.r+1
If arr.r=match Then
ok=1
End
cnt.ok=cnt.ok+1
End
hits=cnt.1/samp
If hits>=.5 Then arrow=' <-'
Else arrow=''
Say format(gs,10) cnt.0 cnt.1 100*hits||'%'||arrow
End
End
Output:
100000 samples . Percentage of groups with at least 2 matches Group size 21 55737 44263 44.26300% 22 52158 47842 47.84200% 23 49141 50859 50.85900% <- 24 46227 53773 53.77300% <- 25 43091 56909 56.90900% <- 100000 samples . Percentage of groups with at least 3 matches Group size 85 52193 47807 47.80700% 86 51489 48511 48.51100% 87 50146 49854 49.85400% 88 48790 51210 51.2100% <- 89 47771 52229 52.22900% <- 100000 samples . Percentage of groups with at least 4 matches Group size 185 50930 49070 49.0700% 186 50506 49494 49.49400% 187 49739 50261 50.26100% <- 188 49024 50976 50.97600% <- 189 48283 51717 51.71700% <- 100000 samples . Percentage of groups with at least 5 matches Group size 311 50909 49091 49.09100% 312 50441 49559 49.55900% 313 49912 50088 50.08800% <- 314 49425 50575 50.57500% <- 315 48930 51070 51.0700% <- 100000 samples . Percentage of groups with at least 6 matches Group size 458 50580 49420 49.4200% 459 49848 50152 50.15200% <- 460 49975 50025 50.02500% <- 461 49316 50684 50.68400% <- 462 49121 50879 50.87900% <-
Ruby
def equalBirthdays(nSharers, groupSize, nRepetitions)
eq = 0
for i in 1 .. nRepetitions
group = [0] * 365
for j in 1 .. groupSize
group[rand(group.length)] += 1
end
eq += group.any? { |n| n >= nSharers } ? 1 : 0
end
return (eq * 100.0) / nRepetitions
end
def main
groupEst = 2
for sharers in 2 .. 5
# Coarse
groupSize = groupEst + 1
while equalBirthdays(sharers, groupSize, 100) < 50.0
groupSize += 1
end
# Finer
inf = (groupSize - (groupSize - groupEst) / 4.0).floor
for gs in inf .. groupSize + 999
eq = equalBirthdays(sharers, groupSize, 250)
if eq > 50.0 then
groupSize = gs
break
end
end
# Finest
for gs in groupSize - 1 .. groupSize + 999
eq = equalBirthdays(sharers, gs, 50000)
if eq > 50.0 then
groupEst = gs
print "%d independant people in a group of %s share a common birthday. (%5.1f)\n" % [sharers, gs, eq]
break
end
end
end
end
main()
- Output:
2 independant people in a group of 23 share a common birthday. ( 50.5) 3 independant people in a group of 90 share a common birthday. ( 53.3) 4 independant people in a group of 187 share a common birthday. ( 50.3) 5 independant people in a group of 313 share a common birthday. ( 50.3)
Scala
import scala.util.Random
import scala.util.control.Breaks._
object Test {
def equalBirthdays(
nSharers: Int,
groupSize: Int,
nRepetitions: Int
): Double = {
val rand = new Random(1)
var eq = 0
for (_ <- 1 to nRepetitions) {
val group = new Array[Int](365)
for (_ <- 1 to groupSize) {
val birthday = rand.nextInt(group.length)
group(birthday) += 1
}
if (group.exists(_ >= nSharers)) eq += 1
}
(eq * 100.0) / nRepetitions
}
def main(args: Array[String]): Unit = {
var groupEst = 2
for (sharers <- 2 until 6) {
// Coarse.
var groupSize = groupEst + 1
while (equalBirthdays(sharers, groupSize, 100) < 50.0)
groupSize += 1
// Finer.
val inf = (groupSize - (groupSize - groupEst) / 4.0).toInt
breakable({
for (gs <- inf until groupSize + 999) {
val eq = equalBirthdays(sharers, groupSize, 250)
if (eq > 50.0) {
groupSize = gs
break
}
}
})
// Finest.
breakable({
for (gs <- (groupSize - 1) until groupSize + 999) {
val eq = equalBirthdays(sharers, gs, 50000)
if (eq > 50.0) {
groupEst = gs
println(
f"$sharers independent people in a group of $gs share a common birthday. ($eq%5.1f)"
)
break
}
}
})
}
}
}
- Output:
2 independent people in a group of 23 share a common birthday. ( 50.6) 3 independent people in a group of 87 share a common birthday. ( 50.4) 4 independent people in a group of 187 share a common birthday. ( 50.1) 5 independent people in a group of 314 share a common birthday. ( 50.2)
SQL
birthday.sql
with
c as
(
select
500 nrep,
50 maxgsiz
from dual
),
reps as
(
select level rep
from dual
cross join c
connect by level <= c.nrep
),
pers as
(
select
round(sqrt(2*level)) npers
from dual
cross join c
connect by level <= c.maxgsiz*(c.maxgsiz+1)/2
),
bds as
(
select
reps.rep,
pers.npers,
floor(dbms_random.value(1,366)) bd
from
reps
cross join pers
),
mtch as
(
select
bds.npers,
case count(distinct bds.bd ) when bds.npers then 0 else 1 end match
from bds
group by
bds.rep,
bds.npers,
null
order by
bds.npers
),
nm as
(
select mtch.npers, sum (mtch.match) nmatch
from mtch
group by mtch.npers
),
sol as
(
select first_value ( nm.npers ) over ( order by abs ( nm.nmatch - c.nrep / 2 ) ) npers
from
nm
cross join c
)
select npers
from sol where rownum = 1
;
SQL> @ birthday.sql Connected.
NPERS
23
Tcl
proc birthdays {num {same 2}} {
for {set i 0} {$i < $num} {incr i} {
set b [expr {int(rand() * 365)}]
if {[incr bs($b)] >= $same} {
return 1
}
}
return 0
}
proc estimateBirthdayChance {num same} {
# Gives a reasonably close estimate with minimal execution time; the idea
# is to keep the amount that one random value may influence the result
# fairly constant.
set count [expr {$num * 100 / $same}]
set x 0
for {set i 0} {$i < $count} {incr i} {
incr x [birthdays $num $same]
}
return [expr {double($x) / $count}]
}
foreach {count from to} {2 20 25 3 85 90 4 183 190 5 310 315} {
puts "identifying level for $count people with same birthday"
for {set i $from} {$i <= $to} {incr i} {
set chance [estimateBirthdayChance $i $count]
puts [format "%d people => %%%.2f chance of %d people with same birthday" \
$i [expr {$chance * 100}] $count]
if {$chance >= 0.5} {
puts "level found: $i people"
break
}
}
}
- Output:
identifying level for 2 people with same birthday 20 people => %43.40 chance of 2 people with same birthday 21 people => %44.00 chance of 2 people with same birthday 22 people => %46.91 chance of 2 people with same birthday 23 people => %53.48 chance of 2 people with same birthday level found: 23 people identifying level for 3 people with same birthday 85 people => %47.97 chance of 3 people with same birthday 86 people => %48.46 chance of 3 people with same birthday 87 people => %49.55 chance of 3 people with same birthday 88 people => %50.66 chance of 3 people with same birthday level found: 88 people identifying level for 4 people with same birthday 183 people => %48.02 chance of 4 people with same birthday 184 people => %47.67 chance of 4 people with same birthday 185 people => %48.89 chance of 4 people with same birthday 186 people => %49.98 chance of 4 people with same birthday 187 people => %50.99 chance of 4 people with same birthday level found: 187 people identifying level for 5 people with same birthday 310 people => %48.52 chance of 5 people with same birthday 311 people => %48.14 chance of 5 people with same birthday 312 people => %49.07 chance of 5 people with same birthday 313 people => %49.63 chance of 5 people with same birthday 314 people => %49.59 chance of 5 people with same birthday 315 people => %51.79 chance of 5 people with same birthday level found: 315 people
Wren
import "random" for Random
import "./fmt" for Fmt
var equalBirthdays = Fn.new { |nSharers, groupSize, nRepetitions|
var rand = Random.new(12345)
var eq = 0
for (i in 0...nRepetitions) {
var group = List.filled(365, 0)
for (j in 0...groupSize) {
var r = rand.int(group.count)
group[r] = group[r] + 1
}
eq = eq + (group.any { |i| i >= nSharers } ? 1 : 0)
}
return eq * 100 / nRepetitions
}
var groupEst = 2
for (sharers in 2..5) {
// Coarse
var groupSize = groupEst + 1
while (equalBirthdays.call(sharers, groupSize, 100) < 50) groupSize = groupSize + 1
// Finer
var inf = (groupSize - (groupSize - groupEst) / 4).floor
for (gs in inf...groupSize + 999) {
var eq = equalBirthdays.call(sharers, groupSize, 250)
if (eq > 50) {
groupSize = gs
break
}
}
// Finest
for (gs in groupSize - 1...groupSize + 999) {
var eq = equalBirthdays.call(sharers, gs, 50000)
if (eq > 50) {
groupEst = gs
Fmt.write("$d independent people in a group of $3d ", sharers, gs)
Fmt.print("share a common birthday $2.1f\%", eq)
break
}
}
}
- Output:
2 independent people in a group of 23 share a common birthday 51.0% 3 independent people in a group of 88 share a common birthday 51.0% 4 independent people in a group of 187 share a common birthday 50.1% 5 independent people in a group of 314 share a common birthday 50.7%
XPL0
func Sim(N);
\Simulate birthdays and return number of people to have N same days
int N, I, People, R;
char Days(365);
[for I:= 0 to 365-1 do Days(I):= 0;
People:= 0;
loop [R:= Ran(365);
Days(R):= Days(R)+1;
People:= People+1;
if Days(R) = N then return People;
];
];
int N, Sum, Trial;
[for N:= 2 to 5 do
[Sum:= 0;
for Trial:= 1 to 10000 do
Sum:= Sum + Sim(N);
IntOut(0, N); Text(0, ": "); IntOut(0, Sum/10000); CrLf(0);
];
]
- Output:
2: 24 3: 88 4: 187 5: 311
zkl
Pure simulation; adding a person to a population until there are the required number of collisions, then repeating that a bunch of times to get an average.
fcn bdays(N){ // N is shared birthdays in a population
year:=(0).pump(365,List.createLong(365).write,0); // 365 days == one year
shared:=people:=0; do{ // add a new person to population
bday:=(0).random(365); // with this birthday [0..364]
shared=shared.max(year[bday]+=1); people+=1;
}while(shared<N);
people // size of simulated population that contains N shared birthdays
}
fcn simulate(N,T){ avg:=0.0; do(T){ avg+=bdays(N) } avg/=T; } // N shared, T trials
foreach n in ([1..5]){
println("Average of %d people in a populatation of %s share birthdays"
.fmt(n,simulate(n,0d10_000)));
}
- Output:
Average of 1 people in a populatation of 1 share birthdays Average of 2 people in a populatation of 24.7199 share birthdays Average of 3 people in a populatation of 88.6416 share birthdays Average of 4 people in a populatation of 186.849 share birthdays Average of 5 people in a populatation of 312.399 share birthdays