Sprawdzanie, czy moneta jest uczciwa

Przyjaciel zadał mi następujące pytanie. Nie mogłem jej pomóc, ale mam nadzieję, że ktoś mi to wyjaśni. Nie mogłem znaleźć podobnego przykładu. Dzięki za wszelką pomoc i wyjaśnienia.

P: Wyniki 100 eksperymentów rzucania monetą są rejestrowane jako 0 = „Ogon” i 1 = „Głowa”. Wyjściowy x jest łańcuchem zer i 1 o długości 100. A liczba razy, gdy otrzymujemy 1-0-0 w x jest obliczana i wynosi 20 (np .: if x = (001001110100), 1-0-0 występuje 2 razy). Czy uważasz, że to uczciwa moneta?

Rozwiązanie problemu przez symulację

Moją pierwszą próbą byłoby symulowanie tego na komputerze, który potrafi bardzo szybko rzucić wiele uczciwych monet. Poniżej znajduje się przykład z milionem prób. Zdarzenie „, w którym liczba wzorca„ 1-0-0 ”występuje w rzutach monetą, wynosi 20 lub więcej” występuje mniej więcej raz na trzy tysiące prób, więc to, co zaobserwowałeś, nie jest bardzo prawdopodobne (uczciwie moneta).$X$$n=100$

Zauważ, że histrogram służy do symulacji, a linia jest dokładnym obliczeniem wyjaśnionym poniżej.

set.seed(1)

# number of trials
n <- 10^6

# flip coins
q <- matrix(rbinom(100*n, 1, 0.5),n)

# function to compute number of 100 patterns
npattern <- function(x) {
  sum((1-x[-c(99,100)])*(1-x[-c(1,100)])*x[-c(1,2)])
}

# apply function on data 
counts <- sapply(1:n, function(x) npattern(q[x,]))
hist(counts, freq = 0) 

# estimated probability
sum(counts>=20)/10^6
10^6/sum(counts>=20)

---

Rozwiązanie problemu za pomocą dokładnych obliczeń

W przypadku podejścia analitycznego można wykorzystać fakt, że „prawdopodobieństwo zaobserwowania 20 lub więcej sekwencji„ 1-0-0 ”na 100 rzutów monetą jest równe 1 minus prawdopodobieństwo, że wykonanie więcej niż 100 rzutów zajmie więcej niż 100 rzutów” . Rozwiązano to w następujących krokach:

Czas oczekiwania na prawdopodobieństwo przewrócenia „1-0-0”

Dystrybucja, $f_{N,x=1}(n)$, z ilu razy musisz przerzucić, dopóki nie otrzymasz dokładnie jednej sekwencji „1-0-0”, można obliczyć w następujący sposób:

Przeanalizujmy sposoby dotarcia do „1-0-0” jako łańcucha Markowa. Śledzimy stany opisane sufiksem ciągu flips: „1”, „1-0” lub „1-0-0”. Na przykład, jeśli masz następujące osiem rzutów 10101100, zdałeś kolejno następujące osiem stanów: „1”, „1-0”, „1”, „1-0”, „1”, „1”, „1-0”, „1-0-0”, a osiągnięcie „1-0-0” wymagało ośmiu rzutów. Pamiętaj, że nie ma jednakowego prawdopodobieństwa osiągnięcia stanu „1-0-0” przy każdym rzucie. Dlatego nie można modelować tego jako rozkładu dwumianowego . Zamiast tego powinieneś podążać za drzewem prawdopodobieństw. Stan „1” może przejść do „1” i „1-0”, stan „1-0” może przejść do „1” i „1-0-0”, a stan „1-0-0” jest stanem absorbującym. Możesz zapisać to jako:

           number of flips
           1   2   3   4   5   6   7   8   9   ....   n

'1'        1   1   2   3   5   8  13  21  34   ....   F_n
'1-0'      0   1   1   2   3   5   8  13  21          F_{n-1}
'1-0-0'    0   0   1   2   4   7   12 20  33          sum_{x=1}^{n-2} F_{x}

i prawdopodobieństwo osiągnięcia wzoru „1-0-0”, po wyrzuceniu pierwszego „1” (zaczynasz od stanu „0”, nie odwracając jeszcze głowy), w ciągu $n$ przewrócenie jest o połowę mniejsze niż prawdopodobieństwo, że znajdzie się w stanie „1-0” $n-1$ klapki:

$$f_{N_c,x=1}(n) = \frac{F_{n-2}}{2^{n-1}}$$

gdzie $F_i$ jest $i$-ty numer Fibonnaci. Prawdopodobieństwo bezwarunkowe jest sumą

$$f_{N,x=1}(n) = \sum_{k=1}^{n-2} 0.5^{k} f_{N_c,x=1}(1+(n-k)) = 0.5^{n} \sum_{k=1}^{n-2} F_{k}$$

Czas oczekiwania na prawdopodobieństwo przewrócenia $k$ razy „1-0-0”

Można to obliczyć za pomocą splotu.

$$f_{N,x=k}(n) = \sum_{l=1}^{n} f_{N,x=1}(l)f_{N,x=1}(n-l)$$

z prawdopodobieństwem zaobserwujesz 20 lub więcej wzorów „1-0-0” (w oparciu o hipotezę, że moneta jest uczciwa)

> # exact computation
> 1-Fx[20]
[1] 0.0003247105
> # estimated from simulation
> sum(counts>=20)/10^6
[1] 0.000337

Oto kod R do jego obliczenia:

# fibonacci numbers
fn <- c(1,1)
for (i in 3:99) {
  fn <- c(fn,fn[i-1]+fn[i-2])
}

# matrix to contain the probabilities
ps <- matrix(rep(0,101*33),33)

# waiting time probabilities to flip one pattern
ps[1,] <- c(0,0,cumsum(fn))/2^(c(1:101))

#convoluting to get the others
for (i in 2:33) {
  for (n in 3:101) {
     for (l in c(1:(n-2))) {
       ps[i,n] = ps[i,n] + ps[1,l]*ps[i-1,n-l]
     }  
  }
}

# cumulative probabilities to get x patterns in n flips
Fx <- 1-rowSums(ps[,1:100])

# probabilities to get x patterns in n flips
fx <- Fx[-1]-Fx[-33]

#plot in the previous histogram
lines(c(1:32)-0.5,fx)

---

Obliczanie nieuczciwych monet

Możemy uogólnić powyższe obliczenie prawdopodobieństwa zaobserwowania $x$ wzory w $n$ odwraca się, gdy prawdopodobieństwo „1 = głowa” wynosi $p$ a klapki są niezależne.

Teraz używamy uogólnienia liczb Fibonacciego:

$$F_{n}(x) = \begin{cases} 1 & \quad \text{if $n=1$} \\ x & \quad \text{if $n=2$} \\ x(F_{n-1}+ F_{n-2}) & \quad \text{if $n>2$} \end{cases}$$

prawdopodobieństwo jest teraz następujące:

$$f_{N_c,x=1,p}(n) = (1-p)^{n-1} F_{n-2}((1-p)^{-1}-1)$$

i

$$f_{N,x=1,p}(n) = \sum_{k=1}^{n-2} p(1-p)^{k-1} f_{N_c,x=1,p}(1+n-k) = p(1-p)^{n-1}\sum_{k=1}^{n-2} F_{k}((1-p)^{-1}-1)$$

Kiedy knujemy to, otrzymujesz:

Chociaż wartość p jest niewielka dla uczciwej monety 0,0003247, musimy zauważyć, że nie jest ona znacznie lepsza (tylko pojedyncze zamówienie) dla różnych nieuczciwych monet. Współczynnik prawdopodobieństwa lub współczynnik Bayesa wynosi około 11, gdy hipoteza zerowa ($p=0.5$) porównuje się z hipotezą alternatywną $p=0.33$. Oznacza to, że tylny iloraz szans jest tylko dziesięć razy wyższy niż wcześniejszy iloraz szans.

Jeśli więc przed eksperymentem myślałeś, że moneta jest mało prawdopodobna, to teraz powinieneś nadal myśleć, że moneta jest nieprawdziwa.

---

Moneta z $p_{heads} = p_{tails}$ ale niesprawiedliwość w odniesieniu do zdarzeń „1-0-0”

O wiele łatwiej byłoby przetestować prawdopodobieństwo uczciwej monety, licząc liczbę głów i ogonów, i zastosować rozkład dwumianowy do modelowania tych obserwacji i sprawdzić, czy obserwacja jest szczególna, czy nie.

Może się jednak zdarzyć, że moneta przerzuca średnio taką samą liczbę głów i ogonów, ale nie jest sprawiedliwa w odniesieniu do niektórych wzorów. Na przykład moneta może mieć pewną korelację dla kolejnych rzutów monetą (wyobrażam sobie jakiś mechanizm z zagłębieniami wewnątrz metalu monety wypełnionymi piaskiem, który popłynie jak klepsydra w kierunku przeciwnego końca poprzedniego rzutu monetą, który ładuje monetę bardziej prawdopodobne, że upadną po tej samej stronie, co poprzednia strona).

Niech pierwsza moneta będzie równa prawdopodobieństwu głów i reszek, a kolejne rzuty są prawdopodobne $p$ta sama strona, co wcześniej flip. Następnie podobna symulacja jak na początku tego postu da następujące prawdopodobieństwa, ile razy wzorzec „1-0-0” przekroczy 20:

Widać, że można nieco zmniejszyć prawdopodobieństwo zaobserwowania wzoru „1-0-0” (gdzieś w pobliżu $p=0.45$moneta, która ma pewną ujemną korelację), ale bardziej dramatyczne jest to, że można znacznie zmniejszyć prawdopodobieństwo przestrzegania wzoru „1-0-0”. Za niskie$p$dostajesz wielokrotnie ogony za głowami, pierwszą część „1-0” wzoru „1-0-0”, ale nie dostajesz tak często dwóch ogonów z rzędu część „0-0” wzór. Odwrotnie jest w przypadku wysokich$p$ wartości.

# number of trials
set.seed(1)
n <- 10^6

p <- seq(0.3,0.6,0.02)
np <- length(p)
mcounts <- matrix(rep(0,33*np),33)

pb <- txtProgressBar(title = "progress bar", min = 0,
                     max = np, style=3)
for (i in 1:np) {
  # flip first coins
  qfirst <- matrix(rbinom(n, 1, 0.5),n)*2-1
  # flip the changes of the sign of the coin
  qrest <- matrix(rbinom(99*n, 1, p[i]),n)*2-1
  # determining the sign of the coins
  qprod <- t(sapply(1:n, function(x) qfirst[x]*cumprod(qrest[x,])))
  # representing in terms of 1s and 0s
  qcoins <- cbind(qfirst,qprod)*0.5+0.5
  counts <- sapply(1:n, function(x) npattern(qcoins[x,]))

  mcounts[,i] <- sapply(1:33, function(x) sum(counts==x))
  setTxtProgressBar(pb, i)
}
close(pb)

plot(p,colSums(mcounts[c(20:33),]),
     type="l", xlab="p same flip", ylab="counts/million trials", 
     main="observation of 20 or more times '1-0-0' pattern \n for coin with correlated flips")
points(p,colSums(mcounts[c(20:33),]))

---

Wykorzystanie matematyki w statystyce

Wszystko powyższe jest w porządku, ale nie jest to bezpośrednia odpowiedź na pytanie

„Czy uważasz, że to uczciwa moneta?”

Aby odpowiedzieć na to pytanie, można użyć powyższej matematyki, ale naprawdę należy najpierw bardzo dobrze opisać sytuację, cele, definicję uczciwości itp. Bez znajomości tła i okoliczności jakiekolwiek obliczenia będą jedynie ćwiczeniem matematycznym, a nie odpowiedzią na wyraźne pytanie.

Jedno otwarte pytanie brzmi: dlaczego i jak szukamy wzoru „1-0-0”.

• Na przykład może ten wzór nie był celem, o którym zdecydowano przed przeprowadzeniem dochodzenia. Być może było to po prostu coś, co „wyróżniało się” w danych i przyciągnęło uwagę po eksperymencie. W takim przypadku należy wziąć pod uwagę, że skutecznie dokonuje się wielu porównań .
• Inną kwestią jest to, że obliczone powyżej prawdopodobieństwo jest wartością p. Znaczenie wartości p należy dokładnie rozważyć. To nie prawdopodobieństwo, że moneta jest fair. Zamiast tego istnieje prawdopodobieństwo zaobserwowania określonego wyniku, jeśli moneta jest uczciwa. Jeśli ktoś ma otoczenie, w którym zna się pewną dystrybucję uczciwości monet, lub można przyjąć rozsądne założenie, wówczas można wziąć to pod uwagę i zastosować wyrażenie bayesowskie .
• Co jest sprawiedliwe, a co niesprawiedliwe. W końcu, biorąc pod uwagę wystarczającą liczbę prób, można znaleźć odrobinę niesprawiedliwości. Ale czy jest to istotne i czy takie wyszukiwanie nie jest stronnicze? Trzymając się podejścia częstokroć, należy opisać coś w rodzaju granicy, powyżej której uważamy monetę za uczciwą (pewien odpowiedni rozmiar efektu). Następnie można zastosować coś podobnego do dwustronnego testu t , aby zdecydować, czy moneta jest uczciwa, czy nie (w odniesieniu do wzoru „1-0-0”).