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

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?

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

Niech jest prawdopodobieństwo sukcesu i k = 8 jest liczba sukcesów w rzędzie chcesz liczyć. Zostały one naprawione dla problemu. Zmienne wartości to m , liczba prób pozostałych w bloku; oraz j , liczba zaobserwowanych sukcesów. Niech szansa na osiągnięcie k sukcesów z rzędu przed wyczerpaniem m prób zostanie zapisana f p , k ( j , m ) . Dążyć F 1 / 3 , 8$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 m > 0 próbami do przejścia. Kolejna próba jest albo sukcesem, z prawdopodobieństwem p - w którym przypadku j wzrasta do j + 1 -; albo jest to awaria, z prawdopodobieństwem 1 - p --W tym przypadku j jest resetowany do 0 . W obu przypadkach m zmniejsza się o 1 . 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).$$

Jako warunki początkowe mamy oczywiste wyniki dla m ≥ 0 ( tj. Widzieliśmy już k w rzędzie) i f p , k ( j , m ) = 0 dla k - j > m ( 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) = 18p^8 - 17p^9 - 45p^{16} + 81p^{17}-36p^{18}.$$

Gdy Daje to 80.897 / 43.046.721 ≈ 0,0018793 .$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 10 6 iteracji, uzyskując 0,001867 : tylko 0,3 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,0018475 .)$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 - f 1 / 3 , 8 ( 0 , 25 ) ) 8 ) = 0,0149358 ... .$1-(1-(1/3)^8)^{136} \approx 0.0205$$1 - (1 - f_{1/3,8}(0,25))^8) = 0.0149358...$

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 m i pożądanej długości próby k, podczas gdy metoda potęgowania macierzy to O ( k 3 log ( m ) ) . Jeśli m jest znacznie większe niż k , 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 1 - p przechodzi do następnego stanu (stanem ostatecznym jest stan pochłaniania). Podnosząc tę ​​macierz do potęgi n , wartość w pierwszym rzędzie i ostatniej kolumnie jest prawdopodobieństwem zobaczenia k = 8 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.

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 n oraz prawidłową próbę przez W (wygrana) i niepoprawną próbę przez F (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 + 1 ):$k$$n$$W$$F$$k+1$

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

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

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

$\ldots$

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

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

Prawdopodobieństwo przejścia do stanu ze stanu A jest P = 1 / 3 i prawdopodobieństwo 1 - P = 2 / 3 pobyt w stan A . Od stanu B , prawdopodobieństwo przejścia do stanu C wynosi 1 / 3 oraz z prawdopodobieństwem 2 / 3 ruszamy z powrotem do A . I tak dalej. Jeśli jesteśmy w stanie I , 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ą M (ponieważ każda kolumna M sumuje się do 1, a wszystkie wpisy są dodatnie, M nazywa się lewą macierzą stochastyczną ):$9\times9$ $M$$M$$1$$M$

$$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$

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.

$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}$$

$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$