Mit welcher Wahrscheinlichkeit ist eine Münze besser als die andere?

Nehmen wir an, wir haben zwei voreingenommene Münzen C1und C2beide haben eine unterschiedliche Wahrscheinlichkeit, den Kopf zu drehen.

Wir werfen C1 n1Zeiten und bekommen H1Köpfe, C2 n2Zeiten und bekommen H2Köpfe. Und wir stellen fest, dass das Verhältnis der Köpfe für eine Münze höher ist als für die andere.

Mit welcher Wahrscheinlichkeit können wir sagen, dass eine Münze besser ist als die andere? (Besser bedeutet hier eine höhere tatsächliche Wahrscheinlichkeit, den Kopf zu drehen).

Es ist einfach, die Wahrscheinlichkeit für diese Beobachtung zu berechnen, da die beiden Münzen gleich sind. Dies kann durch einen genauen Fishers-Test erfolgen . Angesichts dieser Beobachtungen

$$\begin{array} {r|c|c} &\text{coin }1 &\text{coin }2 \\ \hline \text{heads} &H_1 &H_2\\ \hline \text{tails} &n_1-H_1 &n_2-H_2\\\end{array}$$

die Wahrscheinlichkeit , diese Zahlen zu beobachten , während die Münzen gleich die Anzahl von Versuchen gegeben sind $n_1$ , $n_2$ und die Gesamtmenge der Köpfe $H_1+H_2$ ist ,

$$p(H_1, H_2|n_1, n_2, H_1+H_2) = \frac{(H_1+H_2)!(n_1+n_2-H_1-H_2)!n_1!n_2!}{H_1!H_2!(n_1-H_1)!(n_2-H_2)!(n_1+n_2)!}.$$

Aber was Sie fragen, ist die Wahrscheinlichkeit, dass eine Münze besser ist. Da wir uns darüber streiten, wie voreingenommen die Münzen sind, müssen wir einen Bayes'schen Ansatz verwenden , um das Ergebnis zu berechnen. Bitte beachten Sie, dass in der Bayes'schen Inferenz der Begriff Glaube als Wahrscheinlichkeit modelliert wird und die beiden Begriffe austauschbar verwendet werden (s. Bayes'sche Wahrscheinlichkeit ). Wir nennen die Wahrscheinlichkeit, dass die Münze $i$ Köpfe wirft, $p_i$ . Die posteriore Verteilung nach Beobachtung für dieses $p_i$ ist durch den Satz von Bayes gegeben :

$$f(p_i|H_i,n_i)= \frac{f(H_i|p_i,n_i)f(p_i)}{f(n_i,H_i)}$$ Die Wahrscheinlichkeitsdichtefunktion (pdf)$f(H_i|p_i,n_i)$ist durch die Binomialwahrscheinlichkeit gegeben, da die einzelnen Versuche Bernoulli-Experimente sind: Ichdavon aus, dassdas Vorwissen überbesteht, dassmit gleicher Wahrscheinlichkeitirgendwo zwischenundliegen könnte, daher. Der Nominator ist also $$f(H_i|p_i,n_i) = \binom{n_i}{H_i}p_i^{H_i}(1-p_i)^{n_i-H_i}$$ f(pi)pi01f(pi)=1f(Hi|pi,ni)f(pi)=f(Hi|pi,ni$f(p_i)$$p_i$$0$$1$$f(p_i) = 1$$f(H_i|p_i,n_i)f(p_i)= f(H_i|p_i,n_i)$ .

Um zu berechnen wir die Tatsache, dass das Integral über einem PDF eins sein muss . Der Nenner wird also ein konstanter Faktor sein, um genau das zu erreichen. Es ist ein PDF bekannt, das sich vom Nominator nur durch einen konstanten Faktor unterscheidet, nämlich die Beta-Verteilung . Daher ist $f(n_i,H_i)$$\int_0^1f(p|H_i,n_i)\mathrm dp = 1$

$$f(p_i|H_i,n_i) = \frac{1}{B(H_i+1, n_i-H_i+1)}p_i^{H_i}(1-p_i)^{n_i-H_i}.$$

Das pdf für das Wahrscheinlichkeitspaar unabhängiger Münzen ist

$$f(p_1,p_2|H_1,n_1,H_2,n_2) = f(p_1|H_1,n_1)f(p_2|H_2,n_2).$$

Jetzt müssen wir dies über die Fälle integrieren, in denen , um herauszufinden, wie wahrscheinlich Münze besser ist als Münze : $p_1>p_2$$1$$2$

$$\begin{align} \mathbb P(p_1>p_2) &= \int_0^1 \int_0^{p‘_1} f(p‘_1,p‘_2|H_1,n_1,H_2,n_2)\mathrm dp‘_2 \mathrm dp‘_1\\ &=\int_0^1 \frac{B(p‘_1;H_2+1,n_2-H_2+1)}{B(H_2+1,n_2-H_2+1)} f(p‘_1|H_1,n_1)\mathrm dp‘_1 \end{align}$$

Ich kann dieses letzte Integral nicht analytisch lösen, aber man kann es mit einem Computer numerisch lösen, nachdem man die Zahlen eingegeben hat. ist die Beta-Funktion und ist die unvollständige Beta-Funktion. Beachten Sie, dass da eine fortlaufende Variable ist und niemals genau mit identisch ist .$B(\cdot,\cdot)$$B(\cdot;\cdot,\cdot)$$\mathbb P(p_1=p_2) = 0$$p_1$$p_2$

---

In Bezug auf die vorherige Annahme zu und Anmerkungen dazu: Eine gute Alternative zum Modell, von der viele glauben, ist die Verwendung einer Beta-Verteilung . Dies würde zu einer endgültigen Wahrscheinlichkeit führen: Auf diese Weise könnte man eine starke Tendenz zu regulären Münzen im Großen und Ganzen modellieren, aber gleich , . Es wäre gleichbedeutend damit, die Münze zusätzliche Male zu werfen und Köpfe zu empfangen, was gleichbedeutend wäre, nur mehr Daten zu haben. ist die Anzahl der Würfe, die wir nicht machen müssten$f(p_i)$$Beta(a_i+1,b_i+1)$

$$\mathbb P(p_1>p_2) =\int_0^1 \frac{B(p‘_1;H_2+1+a_2,n_2-H_2+1+b_2)}{B(H_2+1+a_2,n_2-H_2+1+b_2)} f(p‘_1|H_1+a_1,n_1+a_1+b_1)\mathrm dp‘_1.$$ $a_i$$b_i$$a_i+b_i$$a_i$$a_i + b_i$ wenn wir dies vorher einschließen.

Das OP gab an, dass beide Münzen in unbekanntem Maße voreingenommen sind. Ich habe also verstanden, dass alles Wissen aus den Beobachtungen abgeleitet werden muss. Aus diesem Grund habe ich mich für einen nicht informativen Vorgänger entschieden , der das Ergebnis nicht in Richtung normaler Münzen verzerrt.

Alle Informationen können in Form von pro Münze übermittelt werden . Das Fehlen eines informativen Prior bedeutet nur, dass mehr Beobachtungen erforderlich sind, um zu entscheiden, welche Münze mit hoher Wahrscheinlichkeit besser ist.$(H_i, n_i)$

---

Hier ist der Code in R, der eine Funktion Verwendung des einheitlichen Prior liefert : P(n1, H1, n2, H2) $=\mathbb P(p_1>p_2)$$f(p_i)=1$

mp <- function(p1, n1, H1, n2, H2) {
    f1 <- pbeta(p1, H2 + 1, n2 - H2 + 1)
    f2 <- dbeta(p1, H1 + 1, n1 - H1 + 1)
    return(f1 * f2)
}

P <- function(n1, H1, n2, H2) {
    return(integrate(mp, 0, 1, n1, H1, n2, H2))
}

Sie können für verschiedene experimentelle Ergebnisse zeichnen und , z. B. mit diesem Code festlegen:$P(p_1>p_2)$$n_1$$n_2$$n_1=n_2=4$

library(lattice)
n1 <- 4
n2 <- 4
heads <- expand.grid(H1 = 0:n1, H2 = 0:n2)
heads$P <- apply(heads, 1, function(H) P(n1, H[1], n2, H[2])$value)
levelplot(P ~ H1 + H2, heads, main = "P(p1 > p2)")

Möglicherweise müssen Sie install.packages("lattice")zuerst.

Man kann sehen, dass selbst mit dem einheitlichen Prior und einer kleinen Stichprobengröße die Wahrscheinlichkeit oder der Glaube, dass eine Münze besser ist, ziemlich solide werden kann, wenn sich und genug unterscheiden. Ein noch kleinerer relativer Unterschied ist erforderlich, wenn und noch größer sind. Hier ist ein Diagramm für und :$H_1$$H_2$$n_1$$n_2$$n_1=100$$n_2=200$

---

Martijn Weterings schlug vor , die posteriore Wahrscheinlichkeitsverteilung für die Differenz zwischen und zu berechnen . Dies kann durch Integrieren des PDFs des Paares über die Menge : $p_1$$p_2$$S(d)=\{(p_1,p_2)\in[0,1]^2|d=|p_1-p_2|\}$

$$\begin{align} f(d|H_1,n_1,H_2,n_2) &= \int_{S(d)}f(p_1,p_2|H_1,n_1,H_2,n_2) \mathrm d\gamma\\ &= \int_0^{1-d} f(p,p+d|H_1,n_1,H_2,n_2) \mathrm dp + \int_d^1 f(p,p-d|H_1,n_1,H_2,n_2) \mathrm dp\\ \end{align}$$

Auch hier kann ich kein Integral analytisch lösen, aber der R-Code wäre:

d1 <- function(p, d, n1, H1, n2, H2) {
    f1 <- dbeta(p, H1 + 1, n1 - H1 + 1)
    f2 <- dbeta(p + d, H2 + 1, n2 - H2 + 1)
    return(f1 * f2)
}

d2 <- function(p, d, n1, H1, n2, H2) {
    f1 <- dbeta(p, H1 + 1, n1 - H1 + 1)
    f2 <- dbeta(p - d, H2 + 1, n2 - H2 + 1)
    return(f1 * f2)
}

fd <- function(d, n1, H1, n2, H2) {
    if (d==1) return(0)
    s1 <- integrate(d1, 0, 1-d, d, n1, H1, n2, H2)
    s2 <- integrate(d2, d, 1, d, n1, H1, n2, H2)
    return(s1$value + s2$value)
}

Ich habe für , , und alle Werte von aufgetragen :$f(d|n_1,H_1,n_2,H_2)$$n_1=4$$H_1=3$$n_2=4$$H_2$

n1 <- 4
n2 <- 4
H1 <- 3
d <- seq(0, 1, length = 500)

get_f <- function(H2) sapply(d, fd, n1, H1, n2, H2)
dat <- sapply(0:n2, get_f)

matplot(d, dat, type = "l", ylab = "Density",
        main = "f(d | 4, 3, 4, H2)")
legend("topright", legend = paste("H2 =", 0:n2),
       col = 1:(n2 + 1), pch = "-")

Sie können die Wahrscheinlichkeit von berechnen über einem Wert von liegen . Beachten Sie, dass die doppelte Anwendung des numerischen Integrals mit einem numerischen Fehler verbunden ist. ZB sollte immer gleich da immer einen Wert zwischen und annimmt . Das Ergebnis weicht jedoch häufig geringfügig ab.$|p_1-p_2|$$d$1 d 0 1integrate(fd, d, 1, n1, H1, n2, H2)integrate(fd, 0, 1, n1, H1, n2, H2)$1$$d$$0$$1$

Ich habe eine numerische Simulation mit durchgeführt R, wahrscheinlich suchen Sie nach einer analytischen Antwort, aber ich dachte, dies könnte interessant zu teilen sein.

set.seed(123)
# coin 1
N1 = 20
theta1 = 0.7

toss1 <- rbinom(n = N1, size = 1, prob = theta1)

# coin 2
N2 = 25
theta2 = 0.5

toss2 <- rbinom(n = N2, size = 1, prob = theta2)

# frequency
sum(toss1)/N1 # [1] 0.65
sum(toss2)/N2 # [1] 0.52

In diesem ersten Code simuliere ich einfach zwei Münzwürfe. Hier kann man natürlich sehen, dass es theta1 > theta2dann natürlich die Frequenz von H1höher sein wird als H2. Beachten Sie die verschiedenen N1, N2Größen.

Mal sehen, was wir mit verschiedenen machen können thetas. Beachten Sie, dass der Code nicht optimal ist. Überhaupt.

simulation <- function(N1, N2, theta1, theta2, nsim = 100) {
  count1 <- count2 <- 0

  for (i in 1:nsim) {
    toss1 <- rbinom(n = N1, size = 1, prob = theta1)
    toss2 <- rbinom(n = N2, size = 1, prob = theta2)

    if (sum(toss1)/N1 > sum(toss2)/N2) {count1 = count1 + 1} 
    #if (sum(toss1)/N1 < sum(toss2)/N2) {count2 = count2 + 1} 
  }

  count1/nsim

}
set.seed(123)
simulation(20, 25, 0.7, 0.5, 100)
#[1] 0.93

0,93 ist also die Häufigkeit (von 100), bei der die erste Münze mehr Köpfe hatte. Dies scheint in Ordnung zu sein, betrachtet theta1und theta2verwendet.

Mal sehen mit zwei Vektoren von thetas.

theta1_v <- seq(from = 0.1, to = 0.9, by = 0.1)
theta2_v <- seq(from = 0.9, to = 0.1, by = -0.1)

res_v <- c()
for (i in 1:length(theta1_v)) {

  res <- simulation(1000, 1500, theta1_v[i], theta2_v[i], 100)
  res_v[i] <- res

}

plot(theta1_v, res_v, type = "l")

Denken Sie daran, dass dies res_vdie Frequenzen H1 > H2von 100 Simulationen sind.

Wenn also theta1zunimmt, H1steigt natürlich die Wahrscheinlichkeit , höher zu sein.

Ich habe einige andere Simulationen gemacht und es scheint , dass die Größen N1, N2sind weniger wichtig.

Wenn Sie mit dem RCode vertraut sind , können Sie ihn mit diesem Code beleuchten. Mir ist bewusst, dass dies keine vollständige Analyse ist und verbessert werden kann.