Możliwe jest wykonanie dokładnych obliczeń.Na przykład w R

  rzuca <-100
odwraca <- 600
dett < - rep (1/6, 6)
for (n in 2: rolls) {
  drazy <- (c (0, densional, 0,0,0,0,0) + c (0,0, densional, 0,0,0,0) + c (0,0,0, densional, 0,0,0) +
            c (0,0,0,0, densional, 0,0) + c (0,0,0,0,0, densional, 0) + c (0,0,0,0,0,0, densional)) / 6}
suma (drazy * (1-pbinom (1: przewroty, przewroty, 1/2))) # prawdopodobieństwo monety więcej
# 0.00809003
suma (drazy * dbinom (1: przerzuty, przewroty, 1/2)) # równość prawdopodobieństwa
# 0,00111972
suma (dic * pbinom (0: (flips-1), flipy, 1/2)) # kości prawdopodobieństwa więcej
# 0.99079025
 

z ostatnią liczbą pasującą do symulacji BruceET

Interesujące części funkcji masy prawdopodobieństwa wyglądają następująco (moneta rzuca się na czerwono, suma kości na niebiesko)

Innym sposobem jest symulacja miliona dopasowań między X USD a Y USD do przybliżenia $ P (X > Y) = 0.9907 \ pm 0,0002. $ [Symulacja w R.]

  set.seed (825)
d = replicate (10 ^ 6, sum (sample (1: 6,100, rep = T)) - rbinom (1,600, .5))
średnia (d > 0)
[1] 0,990736
2 * sd (d > 0) / 1000
[1] 0.0001916057 # aprx 95% margines błędu symulacji
 

Uwagi za komentarz @ AntoniParellada:

W R funkcja sample (1: 6, 100, rep = T) symuluje 100 rzutów uczciwą kostką; suma tego symuluje $ X $ . Również rbinom to kod R do symulacji dwumianowa zmienna losowa; tutaj jest $ Y. $ Różnica to $ D = X - Y. $ Procedura replicate tworzy wektor z milionem różnic d . Wtedy (d > 0) jest wektorem logicznym miliona TRUE si FALSE , czyli mean co jest jego udziałem w TRUE s - naszej odpowiedzi. Wreszcie ostatnia wypowiedź daje margines błędu 95% przedziału ufności proporcji z TRUE s (używając 2 zamiast 1,96), jako rzeczywistego sprawdzenia dokładności symulowanej odpowiedzi. [Z milionem iteracji, których zwykle się oczekuje Dokładność 2 lub 3 dziesiętne dla prawdopodobieństw - czasami większa dla prawdopodobieństwa do tej pory od 1/2.]

Trochę dokładniej:

Wariancja sumy lub różnicy dwóch niezależnych zmiennych losowych jest sumą ich wariancji. Mamy więc rozkład o średniej równej 50 $ i odchyleniu standardowym $ \ sqrt {292 + 150} \ ok. 21 $ . Jeśli chcemy wiedzieć, jak często spodziewamy się, że ta zmienna będzie poniżej 0, możemy spróbować przybliżyć naszą różnicę przez rozkład normalny i musimy sprawdzić $ z $ span> -score for $ z = \ frac {50} {21} \ ok. 2,38 $ . Oczywiście nasza rzeczywista dystrybucja będzie nieco szersza (ponieważ splatamy dwumianowy plik PDF z jednolitym rozkładem PDF), ale miejmy nadzieję, że nie będzie to zbyt niedokładne. Prawdopodobieństwo, że nasza suma będzie dodatnia, zgodnie z tabelą $ z $ -score, wynosi około 0,992 $ span >.

Przeprowadziłem szybki eksperyment w Pythonie, wykonując 10000 iteracji i otrzymałem $ \ frac {9923} {10000} $ pozytywów. Niezbyt daleko.

Mój kod:

  import numpy jako np
c = np.random.randint (0, 2, size = (10000, 100, 6)). sum (oś = -1)
d = np.random.randint (1, 7, size = (10000, 100))
(d.sum (axis = -1) > c.sum (axis = -1)). sum ()
--> 9923
 

Dokładna odpowiedź jest łatwa do obliczenia numerycznego - nie jest wymagana żadna symulacja. Dla celów edukacyjnych, oto podstawowy skrypt Python 3, który to robi, bez użycia gotowych bibliotek statystycznych.

  z kolekcji import defaultdict

# określ dystrybucje pojedynczej monety i kostki
coin = tuple ((i, 1/2) for i in (0, 1))
die = tuple ((i, 1/6) for i in (1, 2, 3, 4, 5, 6))

# prosta funkcja obliczająca sumę dwóch zmiennych losowych
def add_rv (a, b):
  sum = defaultdict (float)
  dla i, p in a:
    dla j, q in b:
      suma [i + j] + = p * q
  return tuple (sum.items ())

# oblicz sumy 600 monet i 100 kostek
coin_sum = dice_sum = ((0, 1),)
dla _ in range (600): coin_sum = add_rv (coin_sum, coin)
for _ in range (100): dice_sum = add_rv (dice_sum, die)

# obliczyć prawdopodobieństwo, że suma kostek będzie wyższa
prob = 0
dla i, p w dice_sum:
  dla j, q w coin_sum:
    jeśli i > j: prob + = p * q

print ("prawdopodobieństwo, że 100 kostek w sumie daje więcej niż 600 monet =% .10f"% prawdopodobieństwa)
 

Wypróbuj online!

Powyższy skrypt przedstawia dyskretny rozkład prawdopodobieństwa jako listę par (wartość, prawdopodobieństwo) i używa prostej pary zagnieżdżonych pętli do obliczenia rozkładu sumy dwóch zmiennych losowych (iterując po wszystkich możliwych wartościach każdej z szczyty). Niekoniecznie jest to najskuteczniejsza możliwa reprezentacja, ale jest łatwa w obsłudze i wystarczająco szybka do tego celu.

(FWIW, ta reprezentacja rozkładów prawdopodobieństwa jest również kompatybilna z zbiorem funkcji użytkowych do modelowania bardziej złożonych rzutów kostką, który napisałem dla posta w naszej siostrzanej witrynie chwilę temu.)

Oczywiście istnieją również biblioteki specyficzne dla domeny, a nawet całe języki programowania do takich obliczeń. Używając jednego z takich narzędzi online, zwanego AnyDice, te same obliczenia można zapisać bardziej zwięźle:

  X: 100d6
R: 600 dni {0,1}
wyjście X > Y o nazwie „1, jeśli X > Y, w przeciwnym razie 0”
 

Pod maską uważam, że AnyDice oblicza wynik prawie tak, jak robi to mój skrypt w Pythonie, może z wyjątkiem nieco większej liczby optymalizacji.W każdym razie obie dają takie samo prawdopodobieństwo 0,9907902497, gdy suma kości jest większa niż liczba orłów.

Jeśli chcesz, AnyDice może również wykreślić rozkład obu sum za Ciebie.Aby uzyskać podobne wykresy z kodu Pythona, należałoby wprowadzić listy dice_sum i coin_sum do biblioteki wykresów, takiej jak pyplot.

Poniższa odpowiedź jest nieco nudna, ale wydaje się być jedyną jak dotąd, która zawiera autentycznie dokładną odpowiedź! Zwykłe przybliżenie lub symulacja, a nawet po prostu liczbowe obliczenie dokładnej odpowiedzi na rozsądny poziom dokładności, co nie zajmuje dużo czasu, jest prawdopodobnie lepszym sposobem - ale jeśli chcesz uzyskać „matematyczny” sposób uzyskania dokładnej odpowiedzi, :

Niech $ X $ oznacza sumę kropek, które widzimy w 100 $ rzutach kostką, z prawdopodobieństwem funkcja masy $ p_X (x) $ .

Niech $ Y $ oznacza liczbę reszek w 600 $ rzutach monetą z funkcją masy prawdopodobieństwa $ p_Y (y) $ .

Szukamy $ P (X > Y) = P (X - Y > 0) = P (D > 0) $ , gdzie $ D = X - Y $ to różnica między sumą kropek a liczbą głowic.

Niech $ Z = -Y $ , z funkcją masy prawdopodobieństwa $ p_Z (z) = p_Y (-z) $ . Wtedy różnica $ D = X - Y $ może zostać przepisana jako suma $ D = X + Z $ span > co oznacza, że ​​skoro $ X $ i $ Z $ są niezależne, możemy znaleźć funkcję masy prawdopodobieństwa $ D $ , biorąc dyskretny splot PMF z $ X $ span > i $ Z $ :

$$ p_D (d) = \ Pr (X + Z = d) = \ sum_ {k = - \ infty} ^ {\ infty} \ Pr (X = k \ cap Z = d - k) = \ sum_ {k = - \ infty} ^ {\ infty} p_X (k) p_Z (dk) $$

W praktyce sumę należy obliczyć tylko na wartościach $ k $ , dla których oczywiście prawdopodobieństwa są różne od zera. Pomysł jest dokładnie tym, co zrobił @IlmariKaronen, chciałem tylko napisać matematyczną podstawę.

Teraz nie powiedziałem, jak znaleźć PMF $ X $ , który jest pozostawiony jako ćwiczenie, ale zwróć uwagę, że jeśli $ X_1, X_2, \ dots, X_ {100} $ to liczba kropek na każdym ze 100 niezależnych rzutów kośćmi, każdy z dyskretnymi, jednolitymi PMF na $ \ {1, 2, 3, 4, 5, 6 \} $ , a następnie $ X = X_1 + X_2 + \ dots + X_ {100} $ span> i tak ...

  # Przechowuj pliki PMF zmiennych jako ramki danych z kolumnami "value" i "prob".
# Ważne wartości są następujące po sobie i rosną dla spójności podczas konwoli,
# więc w razie potrzeby uwzględnij wartości pośrednie z prawdopodobieństwem 0!

# Funkcja sprawdzająca, czy dataframe jest zgodna z powyższą definicją PMF
# Użyj komunikatu message_intro, aby wyjaśnić, na czym polega błąd sprawdzania
is.pmf <- function (x, message_intro = "") {
  if (! is.data.frame (x)) {stop (paste0 (message_intro, "Not a dataframe"))}
  if (! nrow (x) > 0) {stop (paste0 (message_intro, "Dataframe nie ma wierszy"))}
  if (! "value"% in% colnames (x)) {stop (paste0 (message_intro, "Brak kolumny 'wartość'"))}
  if (! "prob"% in% colnames (x)) {stop (paste0 (message_intro, "Brak kolumny 'prob'"))}
  if (! is.numeric (x  $ value)) {stop (paste0 (message_intro, "'wartość' kolumna nie numeryczna"))}
  if (! all (is.finite (x $  value))) {stop (paste0 (message_intro, "Czy 'value' zawiera NA, Inf, NaN itd.?"))}
  if (! all (diff (x  $ value) == 1)) {stop (paste0 (message_intro, "'value' not consecutive and ascending"))}
  if (! is.numeric (x $  prob)) {stop (paste0 (komunikat_intro, "'prob' kolumna nie numeryczna"))}
if (! all (is.finite (x  $ prob))) {stop (paste0 (message_intro, "Czy 'prob' zawiera NA, Inf, NaN itd.?"))}
  if (! all.equal (sum (x $  prob), 1)) {stop (paste0 (message_intro, kolumna "'prob' nie sumuje się do 1"))}
  return (TRUE)
}

# Funkcja do splatania PMF x i y
# Zauważ, że aby zawinąć w R, musimy odwrócić drugi wektor
# name1 i name2 są używane w raportowaniu błędów dla dwóch danych wejściowych
convolve.pmf <- function (x, y, name1 = "x", name2 = "y") {
  is.pmf (x, message_intro = paste0 ("Sprawdzanie", nazwa1, "jest poprawnym PMF:"))
  is.pmf (y, message_intro = paste0 ("Sprawdzanie", nazwa2, "jest poprawnym PMF:"))
  x_plus_y <- data.frame (
    wartość = seq (from = min (x  $ wartość) + min (y $  wartość),
                to = max (x  $ wartość) + max (y $  wartość),
                by = 1),
    prob = convolve (x  $ prob, rev (y $  prob), type = "open")
  )
  powrót (x_plus_y)
}

# Niech x_i będzie wynikiem pojedynczego rzutu kostką i
# Uwaga PMF x_i jest taki sam dla każdego i = 1 do i = 100)
x_i <- data.frame (
  wartość = 1: 6,
  prob = rep (1/6, 6)
)

# Niech t_i będzie sumą x_1, x_2, ..., x_i
# Będziemy przechowywać pliki PMF z t_1, t_2 ... na liście
t_i <- lista ()
t_i [[1]] <- x_i # t_1 to po prostu x_1, więc ma ten sam PMF
# PMF t_i jest splotem PMF t_ (i-1) i x_i
for (i in 2: 100) {
  t_i [[i]] <- convolve.pmf (t_i [[i-1]], x_i,
        nazwa1 = wklej0 ("t_i [[", i-1, "]]"), nazwa2 = "x_i")
}

# Niech x będzie sumą wyników wszystkich 100 niezależnych rzutów kośćmi
x <- t_i [[100]]
is.pmf (x, message_intro = "Sprawdzanie x jest poprawne PMF:")

# Niech y będzie liczbą orłów w 600 rzutach monetą, więc ma rozkład dwumianowy (600, 0,5):
y <- data.frame (wartość = 0: 600)
y  $ prob <- dbinom (wartość y $ , rozmiar = 600, prob = 0,5)
is.pmf (y, message_intro = "Sprawdzanie y jest poprawne PMF:")
# Niech z będzie ujemną wartością y (zauważ, że odwracamy kolejność, aby wartości rosły)
z <- data.frame (value = -rev (y  $ value), prob = rev (y $  prob))
is.pmf (z, message_intro = "Sprawdzanie z jest poprawne PMF:")

# Niech d będzie różnicą, d = x - y = x + z
d <- convolve.pmf (x, z, name1 = "x", name2 = "z")
is.pmf (d, message_intro = "Sprawdzanie, czy d jest poprawne PMF:")

# Prob (X > Y) = Prob (D > 0)
suma (d [d $ wartość > 0, "prob"])
# [1] 0,9907902
 

Wypróbuj online!

Nie ma to znaczenia w praktyce, jeśli zależy Ci na rozsądnej dokładności, ponieważ powyższy kod i tak działa w ułamku sekundy, ale istnieje skrót do zrobienia zwojów dla sumy 100 niezależnych zmiennych o identycznym rozkładzie: ponieważ 100 = 64 + 32 + 4 wyrażone jako suma potęg 2, możesz w miarę możliwości dalej łączyć ze sobą swoje pośrednie odpowiedzi. Zapisywanie sum pośrednich dla pierwszych $ i $ rzutów kośćmi jako $ T_i = \ sum_ {k = 1} ^ {k = i} X_k $ możemy uzyskać PMF z $ T_2 = X_1 + X_2 $ , $ T_4 = T_2 + T_2 '$ (gdzie $ T_2' $ jest niezależne od $ T_2 $ , ale ma ten sam PMF) i podobnie $ T_8 = T_4 + T_4 '$ , $ T_ {16} = T_8 + T_8' $ , $ T_ {32} = T_ {16} + T_ {16} '$ i $ T_ {64} = T_ {32} + T_ {32} '$ . Potrzebujemy jeszcze dwóch zwojów, aby znaleźć łączny wynik wszystkich 100 kości jako sumę trzech niezależnych zmiennych, $ X = T_ {100} = (T_ {64} + T_ {32} '') + T_4 '' $ i końcowe splot dla $ D = X + Z $ . Więc myślę, że w sumie potrzebujesz tylko dziewięciu zwojów - a na koniec możesz ograniczyć się do tych części splotu, które dają dodatnią wartość dla $ D $ . Lub, jeśli jest to mniej kłopotliwe, części, które podają niedodatnie wartości dla $ D $ , a następnie przyjmują dopełnienie. Pod warunkiem, że wybierzesz najbardziej efektywny sposób, myślę, że oznacza to, że Twój najgorszy przypadek to efektywnie osiem i pół zwojów. EDYCJA: i jak sugeruje @whuber, to też niekoniecznie jest optymalne!

Używając zidentyfikowanej przeze mnie metody dziewięciu splotów, z pakietem gmp, mogłem pracować z obiektami bigq i pisać niezoptymalizowaną pętlę do wykonania zwoje (ponieważ wbudowana metoda R nie obsługuje danych wejściowych bigq ), obliczenie dokładnego, uproszczonego ułamka zajęło tylko kilka sekund:

1342994286789364913259466589226414913145071640552263974478047652925028002001448330257335942966819418087658458889485712017471984746983053946540181650207455490497876104509955761041797420425037042000821811370562452822223052224332163891926447848261758144860052289/1355477899826721990460331878897812400287035152117007099242967137806414779868504848322476153909567683818236244909105993544861767898849017476783551366983047536680132501682168520276732248143444078295080865383592365060506205489222306287318639217916612944423026688

co rzeczywiście zaokrągla do 0,9907902. Teraz, aby uzyskać dokładną odpowiedź, nie chciałbym tego robić przy zbyt wielu zwojach, mogłem poczuć, jak koła zębate mojego laptopa zaczynają skrzypieć!