Dies ist eine Verallgemeinerung des berühmten Geburtstagsproblems : Wenn Personen mit zufälligen, gleichmäßig verteilten "Geburtstagen" unter einer Reihe von d = 6 Möglichkeiten sind, wie groß ist die Wahrscheinlichkeit, dass kein Geburtstag von mehr als m = 20 Personen geteilt wird?n=100d=6m=20
Eine genaue Berechnung ergibt die Antwort (mit doppelter Genauigkeit). Ich werde die Theorie skizzieren und den Code für allgemeines n , m , d bereitstellen . Das asymptotische Timing des Codes ist O ( n 2 log ( d ) ) , was ihn für eine sehr große Anzahl von Geburtstagen d geeignet machtund eine angemessene Leistung bietet, bis n in den Tausenden liegt. Zu diesem Zeitpunkt sollte die Poisson-Näherung, die unterAusdehnung des Geburtstagsparadoxons auf mehr als 2 Personenerörtert wurde,in den meisten Fällen gut funktionieren.0.267747907805267n,m,d.O ( n2Log( d) )dn
Erklärung der Lösung
Die Wahrscheinlichkeitserzeugungsfunktion (pgf) für die Ergebnisse von unabhängigen Würfeln eines d- seitigen Würfels istnd
d- nfn( x1, x2, … , X.d) = d- n( x1+ x2+ ⋯ + xd)n.
Der Koeffizient von bei der Erweiterung dieses Multinomials gibt die Anzahl der Möglichkeiten an, wie das Gesicht i genau e i mal erscheinen kann , i = 1 , 2 , … , d .xe11xe22⋯ xeddicheichi=1,2,…,d.
Die Beschränkung unseres Interesses auf nicht mehr als Erscheinungen eines Gesichts ist gleichbedeutend mit der Bewertung von f n modulo, dem Ideal I, das durch x m + 1 1 , x m + 1 2 , … , x m + 1 d erzeugt wird . Um diese Auswertung durchzuführen, verwenden Sie den Binomialsatz rekursiv, um zu erhaltenmfnIxm+11,xm+12,…,xm+1d.
fn(x1,…,xd)=((x1+⋯+xr)+(xr+1+xr+2+⋯+x2r))n=∑k=0n(nk)(x1+⋯+xr)k(xr+1+⋯+x2r)n−k=∑k=0n(nk)fk(x1,…,xr)fn−k(xr+1,…,x2r)
wenn ist gerade. Wenn wir f ( d ) n = f n ( 1 , 1 , … , 1 ) ( d Terme) schreiben , haben wird=2rf(d)n=fn(1,1,…,1)d
f( 2 r )n= ∑k = 0n( nk) f(r)kf(r)n−k.(a)
Wenn ungerade ist, verwenden Sie eine analoge Zerlegungd=2r+1
fn(x1,…,xd)=((x1+⋯+x2r)+x2r+1)n=∑k=0n(nk)fk(x1,…,x2r)fn−k(x2r+1),
geben
f(2r+1)n=∑k=0n(nk)f(2r)kf(1)n−k.(b)
In beiden Fällen können wir auch alles Modulo reduzieren , was von Anfang an leicht durchzuführen istI
fn(xj)≅{xn0n≤mn>mmodI,
Bereitstellung der Startwerte für die Rekursion,
f(1)n={10n≤mn>m
Was diese effizient macht , ist , dass die durch die Spaltung von Variablen in zwei gleich große Gruppen von r Variablen je und Einstellung aller Variablenwerte 1 , wir müssen nur alles beurteilen einmal für eine Gruppe und dann die Ergebnisse kombinieren. Dies erfordert die Berechnung von bis zu n + 1 Termen, von denen jeder eine O ( n ) -Berechnung für die Kombination benötigt. Wir brauchen nicht einmal ein 2D-Array, um f ( r ) n zu speichern , denn bei der Berechnung von f ( d ) n wird nur fdr1,n+1O(n)f(r)nf(d)n, undf ( 1 ) n sind erforderlich.f(r)nf(1)n
Die Gesamtzahl der Schritte ist eins weniger als die Anzahl der Stellen in der binären Erweiterung von (die die Teilungen in Formel ( a ) in gleiche Gruppen zählt ) plus die Anzahl der Stellen in der Erweiterung (die alle Male ungerade zählt Wert angetroffen wird, der die Anwendung der Formel ( b ) erfordert ). Das sind immer noch nur O ( log ( d ) ) Schritte.d(a)(b)O(log(d))
Auf R
einer zehn Jahre alten Workstation war die Arbeit in 0,007 Sekunden erledigt. Der Code ist am Ende dieses Beitrags aufgeführt. Es werden Logarithmen der Wahrscheinlichkeiten anstelle der Wahrscheinlichkeiten selbst verwendet, um mögliche Überläufe oder zu viele Unterläufe zu vermeiden. Dies ermöglicht es, den - Faktor in der Lösung zu entfernen, damit wir die Zählungen berechnen können, die den Wahrscheinlichkeiten zugrunde liegen.d−n
Beachten Sie, dass dieses Verfahren zur Berechnung der gesamten Folge von Wahrscheinlichkeiten auf einmal führt, wodurch wir leicht untersuchen können, wie sich die Chancen mit n ändern .f0,f1,…,fnn
Anwendungen
Die Verteilung im verallgemeinerten Geburtstagsproblem wird von der Funktion berechnet tmultinom.full
. Die einzige Herausforderung besteht darin, eine Obergrenze für die Anzahl der Personen zu finden, die anwesend sein müssen, bevor die Wahrscheinlichkeit einer Kollision zu groß wird. Der folgende Code tut dies mit roher Gewalt, beginnend mit kleinem n und verdoppelt es, bis es groß genug ist. Die gesamte Berechnung benötigt daher O ( n 2 log ( n ) log ( d ) ) Zeit, wobei n die Lösung ist. Die gesamte Wahrscheinlichkeitsverteilung für die Anzahl der Personen bis n wird berechnet.m+1nO(n2log(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)
}
798birthday(7)
365
Code
# 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)
}
Die Antwort erhalten Sie mit
print(tmultinom(100,20,6), digits=15)
0,267747907805267
Brute-Force-Berechnung
Dieser Code dauert auf meinem Laptop einige Sekunden
Ausgabe: 0,2677479
Dennoch könnte es interessant sein, eine direktere Methode zu finden, wenn Sie viele dieser Berechnungen durchführen oder höhere Werte verwenden möchten oder nur um eine elegantere Methode zu erhalten.
Zumindest ergibt diese Berechnung eine einfach berechnete, aber gültige Zahl, um andere (kompliziertere) Methoden zu überprüfen.
quelle