Birthday problem: Difference between revisions
Rename Perl 6 -> Raku, alphabetize, minor clean-up
SqrtNegInf (talk | contribs) (Added Perl example) |
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
||
Line 466:
4 collision: 187 people, P = 0.503229 +/- 0.000645587
5 collision: 313 people, P = 0.501105 +/- 0.000221016</pre>
=={{header|Go}}==
{{trans|C}}
Line 986 ⟶ 987:
find(4)
find(5)</lang>
=={{header|Perl}}==
{{trans|Perl 6}}
<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</lang>
{{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}}
<lang 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; -- Currently 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</lang>
{{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}}==
Line 1,112 ⟶ 1,236:
1182 496957 503043 50.304% <-
1183 494414 505586 50.559% <-</pre>
=={{header|Python}}==
Line 1,470 ⟶ 1,448:
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.
<lang perl6>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;</lang>
{{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}}==
| |||