Jest to uogólnienie znanego problemu urodzinowego : biorąc pod uwagę n=100 osób, które mają losowe, równomiernie rozmieszczone „urodziny” wśród zestawu d=6 możliwości, jaka jest szansa, że żadne urodziny nie są dzielone przez więcej niż m=20 osób?
Dokładne obliczenia dają odpowiedź (z podwójną precyzją). Naszkicuję teorię i podam kod dla ogólnych n , m , d . Czas asymptotyczny kodu wynosi O ( n 2 log ( d ) ), co czyni go odpowiednim dla bardzo dużej liczby urodzin d i zapewnia rozsądną wydajność, dopóki n nie będzie tysiąca. W tym momencie aproksymacja Poissona omawiana podczasrozszerzania paradoksu urodzinowego na więcej niż 2 osobypowinna w większości przypadków działać dobrze.0,267747907805267n , m , d.O ( n2)log( d) )ren
Objaśnienie rozwiązania
Funkcja generowania prawdopodobieństwa (PGF) w odniesieniu do wyników niezależne rolki o d -sided dyszowej mieścinre
d−nfn( x1, x2), … , Xre) = d- n( x1+ x2)+ ⋯ + xre)n.
Współczynnik w rozwinięciu tego wielomianu daje liczbę sposobów, w jakie twarz i może pojawić się dokładnie e i razy, i = 1 , 2 , … , d .xmi11xmi2)2)⋯ xmirerejamijai = 1 , 2 , … , d.
Ograniczenie naszego zainteresowania nie więcej niż występy dowolnej powierzchni jest równoznaczne oceniającej f n modulo idealne I generowany przez x m + 1 1 , x m + 1, 2 , ... , x m + 1 d . Aby wykonać tę ocenę, użyj rekurencyjnie twierdzenia dwumianowego w celu uzyskaniamfanjaxm + 11, xm + 12), … , Xm + 1re.
fan( x1, … , Xre)= ( ( x1+ ⋯ + xr) + ( xr + 1+ xr + 2+ ⋯ + x2 r) )n= ∑k = 0n( nk) (x1+ ⋯ + xr)k( xr + 1+ ⋯ + x2 r)n - k= ∑k = 0n( nk) fk( x1, … , Xr) fn - k( xr + 1, … , X2 r)
gdy jest parzyste. Pisząc ( terminy), mamyf ( d ) n = f n ( 1 , 1 , … , 1 ) dre= 2 rfa( d)n= fn( 1 , 1 , … , 1 )re
fa( 2 r )n= ∑k = 0n( nk) f( r )kfa( r )n - k.(za)
Gdy jest nieparzyste, zastosuj analogiczny rozkładre= 2 r + 1
fan( x1, … , Xre)= ( ( x1+ ⋯ + x2 r) + x2 r + 1)n= ∑k = 0n( nk) fk( x1, … , X2 r) fn - k( x2 r + 1) ,
dający
fa( 2 r + 1 )n= ∑k = 0n( nk) f(2r)kf(1)n−k.(b)
W obu przypadkach możemy również zredukować wszystko modulo , co można łatwo przeprowadzić od początkuI
fn( xjot) ≅{ xn0n ≤ mn > mmodja,
podając wartości początkowe rekurencji,
fa( 1 )n= { 10n ≤ mn > m
To, co sprawia, że jest to wydajne, polega na tym, że dzieląc zmiennych na dwie równe grupy zmiennych każda i ustawiając wszystkie wartości zmiennych na musimy ocenić wszystko tylko raz dla jednej grupy, a następnie połączyć wyniki. Wymaga to obliczenia do terminów, z których każdy wymaga obliczenia dla kombinacji. Nie potrzebujemy nawet tablicy 2D do przechowywania , ponieważ podczas obliczania wymagane są tylko i .r 1 , n + 1 O ( n ) f ( r ) n f ( d ) n , f ( r ) n f ( 1 ) nrer1 ,n + 1O ( n )fa( r )nfa( d)n,fa( r )nfa( 1 )n
Całkowita liczba kroków jest o jeden mniejsza niż liczba cyfr w binarnym rozwinięciu (który liczy podziały na równe grupy we wzorze ) plus liczba jednych w rozwinięciu (które liczą wszystkie razy nieparzyste napotkano wartość wymagającą zastosowania wzoru ). To wciąż tylko kroki .( a ) ( b ) O ( log ( d ) )re( )( b )O ( log( d) )
Na R
dziesięcioletniej stacji roboczej praca została wykonana w 0,007 sekundy. Kod znajduje się na końcu tego postu. Wykorzystuje logarytmy prawdopodobieństwa, a nie same prawdopodobieństwa, aby uniknąć ewentualnych przelewów lub nagromadzenia zbyt dużego niedomiaru. Umożliwia to usunięcie czynnika w rozwiązaniu, dzięki czemu możemy obliczyć liczby leżące u podstaw prawdopodobieństw.re- n
Zauważ, że ta procedura skutkuje obliczeniem całej sekwencji prawdopodobieństwa jednocześnie, co pozwala nam łatwo zbadać, jak szanse zmieniają się z . nfa0, f1, … , Fnn
Aplikacje
Rozkład w uogólnionym problemie urodzinowym jest obliczany przez funkcję tmultinom.full
. Jedyne wyzwanie polega na znalezieniu górnej granicy dla liczby osób, które muszą być obecne, zanim szansa kolizji stanie się zbyt duża. Poniższy kod robi to brutalną siłą, zaczynając od małego i podwajając go, aż będzie wystarczająco duży. Całe obliczenie zajmuje zatem czas gdzie jest rozwiązaniem. Obliczany jest cały rozkład prawdopodobieństwa dla liczby osób w górę przez .n O ( n 2 log ( n ) log ( d ) ) n nm + 1nO ( n2)log( n ) log( d) )nn
#
# The birthday problem: find the number of people where the chance of
# a collision of `m+1` birthdays first exceeds `alpha`.
#
birthday <- function(m=1, d=365, alpha=0.50) {
n <- 8
while((p <- tmultinom.full(n, m, d))[n] > alpha) n <- n * 2
return(p)
}
Przykładowo, jak wynika z obliczeń , minimalna liczba osób potrzebnych w tłumie, aby zwiększyć prawdopodobieństwo, że co najmniej osiem z nich ma wspólne urodziny, wynosi . To zajmuje tylko kilka sekund. Oto wykres części wyniku:798birthday(7)
Specjalna wersja tego problemu dotyczy rozszerzenia paradoksu urodzinowego na więcej niż 2 osoby , który dotyczy przypadku stronnej kostki, która jest rzucana bardzo wiele razy.365
Kod
# Compute the chance that in `n` independent rolls of a `d`-sided die,
# no side appears more than `m` times.
#
tmultinom <- function(n, m, d, count=FALSE) tmultinom.full(n, m, d, count)[n+1]
#
# Compute the chances that in 0, 1, 2, ..., `n` independent rolls of a
# `d`-sided die, no side appears more than `m` times.
#
tmultinom.full <- function(n, m, d, count=FALSE) {
if (n < 0) return(numeric(0))
one <- rep(1.0, n+1); names(one) <- 0:n
if (d <= 0 || m >= n) return(one)
if(count) log.p <- 0 else log.p <- -log(d)
f <- function(n, m, d) { # The recursive solution
if (d==1) return(one) # Base case
r <- floor(d/2)
x <- double(f(n, m, r), m) # Combine two equal values
if (2*r < d) x <- combine(x, one, m) # Treat odd `d`
return(x)
}
one <- c(log.p*(0:m), rep(-Inf, n-m)) # Reduction modulo x^(m+1)
double <- function(x, m) combine(x, x, m)
combine <- function(x, y, m) { # The Binomial Theorem
z <- sapply(1:length(x), function(n) { # Need all powers 0..n
z <- x[1:n] + lchoose(n-1, 1:n-1) + y[n:1]
z.max <- max(z)
log(sum(exp(z - z.max), na.rm=TRUE)) + z.max
})
return(z)
}
x <- exp(f(n, m, d)); names(x) <- 0:n
return(x)
}
Odpowiedź uzyskuje się za pomocą
print(tmultinom(100,20,6), digits=15)
0,267747907805267
Obliczanie siły brutalnej
Ten kod zajmuje kilka sekund na moim laptopie
wyjście: 0,2677479
Ale nadal może być interesujące znalezienie bardziej bezpośredniej metody na wypadek, gdybyś chciał wykonać wiele z tych obliczeń lub użyć wyższych wartości, lub po prostu w celu uzyskania bardziej eleganckiej metody.
Przynajmniej to obliczenie daje uproszczoną, ale prawidłową liczbę do sprawdzenia innych (bardziej skomplikowanych) metod.
źródło