Prawdopodobieństwo serii k sukcesów w sekwencji n prób Bernoulliego

13

Próbuję znaleźć prawdopodobieństwo prawidłowego wykonania 8 prób z rzędu w bloku 25 prób, masz 8 wszystkich bloków (z 25 prób), aby uzyskać 8 prób z rzędu. Prawdopodobieństwo, że jakakolwiek próba będzie poprawna w oparciu o zgadywanie, wynosi 1/3, po uzyskaniu poprawności 8 z rzędu bloki się zakończą (więc uzyskanie więcej niż 8 z rzędu poprawności nie jest technicznie możliwe). Jak mógłbym zająć się znalezieniem prawdopodobieństwa tego wystąpienia? Zastanawiałem się nad zastosowaniem (1/3) ^ 8 jako prawdopodobieństwa prawidłowego otrzymania 8 z rzędu, istnieje 17 możliwych szans na uzyskanie 8 z rzędu w bloku 25 prób, jeśli pomnożę 17 możliwości * 8 bloków Dostaję 136, czy 1- (1- (1/3) ^ 8) ^ 136 dałoby mi prawdopodobieństwo uzyskania 8 w rzędzie poprawnej w tej sytuacji, czy też brakuje mi tutaj czegoś fundamentalnego?

probability binomial AcidNynex
źródło

1

Uważam, że problem z podanym argumentem polega na tym, że rozpatrywane zdarzenia nie są niezależne. Rozważmy na przykład pojedynczy blok. Jeśli powiem ci, że (a) nie ma ósemki rozpoczynającej się od pozycji 6, (b) jest to bieg rozpoczynający się od pozycji 7 i (c) nie ma biegu rozpoczynającego się od pozycji 8, co to oznacza o prawdopodobieństwo uruchomienia rozpoczynającego się od pozycji, powiedzmy, od 9 do 15?

kardynał

14

Śledząc rzeczy, możesz uzyskać dokładną formułę .

Niech jest prawdopodobieństwo sukcesu i jest liczba sukcesów w rzędzie chcesz liczyć. Zostały one naprawione dla problemu. Zmienne wartości to , liczba prób pozostałych w bloku; oraz , liczba zaobserwowanych sukcesów. Niech szansa na osiągnięcie sukcesów z rzędu przed wyczerpaniem prób zostanie zapisana . Dążyć $p=1/3$ $k=8$ $m$ $j$ $k$ $m$ $f_{p,k}(j,m)$ . $f_{1/3,8}(0,25)$

Załóżmy, że właśnie widzieliśmy nasz sukces z rzędu z próbami do przejścia. Kolejna próba jest albo sukcesem, z prawdopodobieństwem którym przypadku wzrasta do -; albo jest to awaria, z prawdopodobieństwem --W tym przypadku jest resetowany do . W obu przypadkach zmniejsza się o . Skąd $j^\text{th}$ $m\gt0$ $p$ $j$ $j+1$ $1-p$ $j$ $0$ $m$ $1$

f_{p, k} (j, m) = p f_{p, k} (j + 1, m - 1) + (1 - p) f_{p, k} (0, m - 1) .

$f_{p,k}(j,m) = p f_{p,k}(j+1,m-1) + (1-p)f_{p,k}(0,m-1).$

Jako warunki początkowe mamy oczywiste wyniki dla ( tj. Widzieliśmy już w rzędzie) i dla ( tzn. nie ma wystarczającej liczby prób, aby uzyskać $f_{p,k}(k,m)=1$ $m \ge 0$ $k$ $f_{p,k}(j,m)=0$ $k-j \gt m$ $k$ z rzędu). Jest teraz szybki i prosty (przy użyciu programowania dynamicznego lub, ponieważ parametry tego problemu są tak małe, rekurencja) do obliczeń

f_{p, 8} (0, 25) = 18 p^{8} - 17 p^{9} - 45 p^{16} + 81 p^{17} - 36 p^{18} .

$f_{p,8}(0,25) = 18p^8 - 17p^9 - 45p^{16} + 81p^{17}-36p^{18}.$

Gdy Daje to . $p=1/3$ $80897 / 43046721 \approx 0.0018793$

Jest to stosunkowo szybki Rkod do symulacji

hits8 <- function() {
    x <- rbinom(26, 1, 1/3)                # 25 Binomial trials
    x[1] <- 0                              # ... and a 0 to get started with `diff`
    if(sum(x) >= 8) {                      # Are there at least 8 successes?
        max(diff(cumsum(x), lag=8)) >= 8   # Are there 8 successes in a row anywhere?
    } else {
        FALSE                              # Not enough successes for 8 in a row
    }
}
set.seed(17)
mean(replicate(10^5, hits8()))

Po 3 sekundach obliczeń wyjście wynosi . Chociaż wygląda to wysoko, to tylko 1,7 standardowych błędów jest wyłączonych. Przeprowadziłem kolejne iteracji, uzyskując : tylko błędy standardowe mniej niż oczekiwano. (Jako podwójną kontrolę, ponieważ wcześniejsza wersja tego kodu zawierała subtelny błąd, uruchomiłem również 400 000 iteracji w Mathematica, uzyskując szacunkową wartość .) $0.00213$ $10^6$ $0.001867$ $0.3$ $0.0018475$

Wynik ten jest mniejszy niż jedna dziesiąta Oszacowanie w pytaniu. Ale może jeszcze nie w pełni zrozumiałe go: kolejna interpretacja „masz 8 Wszystkie bloki ... aby uzyskać 8 prób skorygowania w rzędzie” jest to, że istota odpowiedź poszukiwane jest równa . $1-(1-(1/3)^8)^{136} \approx 0.0205$ $1 - (1 - f_{1/3,8}(0,25))^8) = 0.0149358...$

Whuber
źródło

13

Chociaż doskonałe rozwiązanie do programowania dynamicznego @ Whuber jest warte przeczytania, jego czas działania wynosi w odniesieniu do całkowitej liczby prób i pożądanej długości próby podczas gdy metoda potęgowania macierzy to . Jeśli jest znacznie większe niż , następująca metoda jest szybsza. $\mathcal O(k^2m)$ $m$ $k$ $\mathcal O(k^3\log(m))$ $m$ $k$

Oba rozwiązania traktują problem jako łańcuch Markowa ze stanami reprezentującymi do tej pory liczbę poprawnych prób na końcu łańcucha oraz stanem do osiągnięcia pożądanych poprawnych prób z rzędu. Macierz przejścia jest taka, że zobaczenie awarii z prawdopodobieństwem odsyła cię z powrotem do stanu 0, a w przeciwnym razie z prawdopodobieństwem przechodzi do następnego stanu (stanem ostatecznym jest stan pochłaniania). Podnosząc tę macierz do potęgi , wartość w pierwszym rzędzie i ostatniej kolumnie jest prawdopodobieństwem zobaczenia główek z rzędu. W Pythonie: $p$ $1-p$ $n$ $k=8$

import numpy as np

def heads_in_a_row(flips, p, want):
    a = np.zeros((want + 1, want + 1))
    for i in range(want):
        a[i, 0] = 1 - p
        a[i, i + 1] = p
    a[want, want] = 1.0
    return np.linalg.matrix_power(a, flips)[0, want]

print(heads_in_a_row(flips=25, p=1.0 / 3.0, want=8))

daje pożądane 0,00187928367413.

Neil G.
źródło

10

Zgodnie z tą odpowiedzią wyjaśnię podejście Markov-Chain autorstwa @Neil G i przedstawię ogólne rozwiązanie takich problemów w R. Oznaczmy pożądaną liczbę poprawnych prób z rzędu przez , liczbę prób jako oraz prawidłową próbę przez (wygrana) i niepoprawną próbę przez (niepowodzenie). W trakcie śledzenia prób chcesz wiedzieć, czy masz już serię 8 poprawnych prób i liczbę poprawnych prób na końcu bieżącej sekwencji. Istnieje 9 stanów ( ): $k$ $n$ $W$ $F$ $k+1$

: Nie mieliśmy poprawnych prób jeszcze w rzędzie, a ostatnia próba była . $A$ $8$ $F$

: Nie mieliśmy poprawnych prób jeszcze w rzędzie, a dwie ostatnie próby były . $B$ $8$ $FW$

: Nie miał właściwych prób w rzędzie jednak, a ostatnie trzy próby były . $C$ $8$ $FWW$

$\ldots$

: Nie miał właściwych prób w rzędzie a i osiem ostatnich badaniach były . $H$ $8$ $FWWWWWWW$

: Przeprowadziliśmy poprawnych prób z rzędu! $I$ $8$

Prawdopodobieństwo przejścia do stanu ze stanu jest i prawdopodobieństwo pobyt w stan . Od stanu , prawdopodobieństwo przejścia do stanu wynosi oraz z prawdopodobieństwem ruszamy z powrotem do . I tak dalej. Jeśli jesteśmy w stanie , zostajemy tam. $B$ $A$ $p=1/3$ $1-p=2/3$ $A$ $B$ $C$ $1/3$ $2/3$ $A$ $I$

Na tej podstawie możemy zbudować macierz przejściową (ponieważ każda kolumna sumuje się do a wszystkie wpisy są dodatnie, nazywa się lewą macierzą stochastyczną ): $9\times9$ $M$ $M$ $1$ $M$

M = (\begin{matrix} 2 / 3 & 2 / 3 & 2 / 3 & 2 / 3 & 2 / 3 & 2 / 3 & 2 / 3 & 2 / 3 & 0 \\ 1 / 3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 / 3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 / 3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 / 3 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 / 3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 / 3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 / 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 / 3 & 1 \end{matrix})

$M= \begin{pmatrix} 2/3 & 2/3 & 2/3 & 2/3 & 2/3 & 2/3 & 2/3 & 2/3 & 0 \\ 1/3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1/3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1/3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/3 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1/3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1/3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1/3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1/3 & 1 \end{pmatrix}$

$n$ $M^{n}$ $j$ $i$ $n$ $I$ $1$ $I$ $I$ $1$ $I$ $A$ $n=25$ $M^{25}$ $M^{25}_{91}$ $M^{25}$ Rexpm

library(expm)

k <- 8   # desired number of correct trials in a row
p <- 1/3 # probability of getting a correct trial
n <- 25  # Total number of trials 

# Set up the transition matrix M

M <- matrix(0, k+1, k+1)

M[ 1, 1:k ] <- (1-p)

M[ k+1, k+1 ] <- 1

for( i in 2:(k+1) ) {

  M[i, i-1] <- p

}

# Name the columns and rows according to the states (A-I)

colnames(M) <- rownames(M) <- LETTERS[ 1:(k+1) ]

round(M,2)

     A    B    C    D    E    F    G    H I
A 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0
B 0.33 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0
C 0.00 0.33 0.00 0.00 0.00 0.00 0.00 0.00 0
D 0.00 0.00 0.33 0.00 0.00 0.00 0.00 0.00 0
E 0.00 0.00 0.00 0.33 0.00 0.00 0.00 0.00 0
F 0.00 0.00 0.00 0.00 0.33 0.00 0.00 0.00 0
G 0.00 0.00 0.00 0.00 0.00 0.33 0.00 0.00 0
H 0.00 0.00 0.00 0.00 0.00 0.00 0.33 0.00 0
I 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.33 1

# Calculate M^25

Mn <- M%^%n
Mn[ (k+1), 1 ]
[1] 0.001879284

$A$ $I$ $0.001879284$

COOLSerdash
źródło

3

Oto kod R, który napisałem, aby to zasymulować:

tmpfun <- function() {
     x <- rbinom(25, 1, 1/3)  
     rx <- rle(x)
     any( rx$lengths[ rx$values==1 ] >= 8 )
}

tmpfun2 <- function() {
    any( replicate(8, tmpfun()) )
}

mean(replicate(100000, tmpfun2()))

Dostaję wartości nieco mniejsze niż twoja formuła, więc jedno z nas mogło gdzieś popełnić błąd.

Greg Snow
źródło

Czy twoja funkcja obejmuje próby, w których niemożliwe jest uzyskanie 8 z rzędu, np. Gdzie „bieg” rozpoczął się na próbie 20?

Michelle,

Najprawdopodobniej moja symulacja R daje mi również mniejsze wartości. Jestem ciekawy, czy istnieje rozwiązanie algebraiczne, które rozwiązałoby to jako prosty problem prawdopodobieństwa na wypadek, gdyby ktoś zakwestionował symulację.

AcidNynex

1

Myślę, że ta odpowiedź zostałaby poprawiona poprzez dostarczenie uzyskanych wyników, aby można je było porównać. Oczywiście dodanie czegoś takiego jak histogram byłoby jeszcze lepsze! Kod wygląda mi na pierwszy rzut oka. Twoje zdrowie. :)

kardynał

3

$1$ $0$

M = Table[e[i, j] /. {
    e[9, 1] :> 0,
    e[9, 9] :> 1,
    e[_, 1] :> (1 - p),
    e[_, _] /; j == i + 1 :> p,
    e[_, _] :> 0
  }, {i, 1, 9}, {j, 1, 9}];

x = MatrixPower[M, 25][[1, 9]] // Expand

18 p^{8} - 17 p^{9} - 45 p^{16} + 81 p^{17} - 36 p^{18}

$18 p^8 - 17 p^9 - 45 p^{16} + 81 p^{17} - 36 p^{18}$

$p=\frac{1.0}{3.0}$

x /. p -> 1/3 // N

$0.00187928$

Można to również ocenić bezpośrednio za pomocą funkcji wbudowanych Probabilityi DiscreteMarkovProcess Mathematica :

Probability[k[25] == 9, Distributed[k, DiscreteMarkovProcess[1, M /. p -> 1/3]]] // N

$0.00187928$

Hossam Karim
źródło

Prawdopodobieństwo serii k sukcesów w sekwencji n prób Bernoulliego

Odpowiedzi: