Warum ist diese Verteilung einheitlich?

Wir untersuchen statistische Tests nach Bayes und stoßen auf ein merkwürdiges (zumindest für mich) Phänomen.

Betrachten Sie den folgenden Fall: Wir sind daran interessiert zu messen, welche Population A oder B eine höhere Conversion-Rate aufweist. Für eine Plausibilitätsprüfung setzen wir $p_A = p_B$ , dh die Konversionswahrscheinlichkeit ist in beiden Gruppen gleich. Wir erzeugen künstliche Daten mit einem Binomialmodell, zB

n_{A} \sim Binomial (N, p_{A})

$n_A \sim \text{Binomial}(N, p_A)$

Wir versuchen dann, das $p_A, p_B$ Verwendung eines Bayes'schen Beta-Binomial-Modells zu schätzen sodass wir Posterioren für jede Conversion-Rate erhalten, z. B.

P_{A} \sim Beta (1 + n_{A}, N - n_{A} + 1)

$P_A \sim \text{Beta}(1 + n_A, N - n_A +1 )$

Unsere Teststatistik wird durch Berechnung von berechnet $S = P(P_A > P_B\; |\; N, n_A, n_B)$ über Monte Carlo.

Was mich überraschte, war, dass wenn $p_A = p_B$ , dann $S \sim \text{Uniform(0,1)}$ . Meine Gedanken waren, dass es um 0,5 zentriert sein und sogar auf 0,5 konvergieren würde, wenn die Stichprobengröße $N$ wächst.

Meine Frage ist, warum ist $S \sim \text{Uniform(0,1)}$ wenn $p_A = p_B$ ?

Hier ist ein Python-Code zur Veranschaulichung:

%pylab
from scipy.stats import beta
import numpy as np
import pylab as P

a = b = 0.5
N = 10000
samples = [] #collects the values of S
for i in range(5000):
    assert a==b
    A = np.random.binomial(N, a); B = np.random.binomial(N, b)
    S = (beta.rvs(A+1, N-A+1, size=15000) > beta.rvs(B+1, N-B+1, size=15000)).mean() 
    samples.append(S)

P.hist(samples)
P.show()

hypothesis-testing bayesian binomial python beta-distribution Cam.Davidson.Pilon
quelle

Beachten Sie, dass

nicht exakt gleichförmig sein kann, da es sich um eine diskrete Variable handelt. Sie fragen daher nach asymptotischem Verhalten. Darüber hinaus ist die Verteilung für kleines

(weniger als

, ungefähr mit

) nicht einmal annähernd gleichförmig.

S

$S$

N

$N$

100 / min (p, 1 - p)

$100/\min(p,1-p)$

p = p_{A} = p_{B}

$p=p_A=p_B$

whuber

@whuber S ist nicht diskret, es ist eine Wahrscheinlichkeit, die zwischen 0 und 1 liegen kann. Auch bei niedrigem N beobachte ich ein einheitliches Verhalten.

Cam.Davidson.Pilon

Ich muss also dein Setup missverstehen. Soweit ich das beurteilen kann

der Wert von

für alle gegebenen Werte von

eine Zahl. Daher Annahme , dass

und

sind für den Moment festgelegt (wie sie im Code sind),

ist eine Funktion der

. Letzteres kann jedoch als Realisierung von zwei Binomialverteilungen nur eine diskrete Menge von Werten erreichen. Wenn ich Ihren Code wiedergebe

N, n_{A}, n_{B},

$N,n_A,n_B,$

S

$S$

N, p_{A},

$N, p_A,$

p_{B}

$p_B$

S

$S$

(n_{A}, n_{B})

$(n_A,n_B)$ RIch entschieden nicht-einheitliche Histogramme für kleine bekommen

N

$N$

whuber

Obwohl Ihr

tatsächlich Werte zwischen

und

, verwechseln Sie das nicht mit nicht-diskret: Es kann höchstens

verschiedene Werte haben (und tatsächlich weniger als diese). Dies ist für Sie möglicherweise nicht ganz klar, da Ihre Simulation Schätzungen von

anstelle der korrekten Werte generiert und die Schätzungen im Wesentlichen eine kontinuierliche Verteilung aufweisen.

S

$S$

0

$0$

1

$1$

N^{2}

$N^2$

S

$S$

whuber

@whuber ja, du hast recht, ausgezeichnete Beobachtung. Ich bin immer noch gespannt, warum es dann einheitlich aussieht .

Cam.Davidson.Pilon

Antworten:

TL; DR: Gemische mit Normalverteilungen können bei großen Behältern einheitlich aussehen.

Diese Antwort stammt aus @ whubers Beispielcode (was ich zuerst für einen Fehler hielt, aber im Nachhinein wahrscheinlich für einen Hinweis).

Die zugrunde liegenden Verhältnisse in der Bevölkerung sind gleich: a = b = 0.5.
Jede Gruppe, A und B beträgt 10000 Mitglieder: N = 10000.
Wir werden 5000 Replikaten einer Simulation durchzuführen: for i in range(5000):.

Tatsächlich das , was wir tun , ist a von a . In jedem der 5000 Iterationen wir tun $\rm simulation_\rm{prime}$ $\rm simulation_\rm{underlying}$ $\rm simulation_\rm{prime}$ . $\rm simulation_\rm{underlying}$

In jeder Iteration von werden wir eine Zufallszahl von A und B simulieren , die 'Erfolge' (AKA umgewandelt) gegeben sind die gleichen zugrundeliegenden Anteile zuvor definiert: . Normalerweise ergibt dies A = 5000 und B = 5000, aber A und B variieren von Simulationslauf zu Simulationslauf und werden unabhängig und (ungefähr) normal auf die 5000 Simulationsläufe verteilt (wir werden darauf zurückkommen). $\rm simulation_\rm{prime}$ A = np.random.binomial(N, a); B = np.random.binomial(N, b)

Lassen Sie sich nun Schritt für Schritt durch für eine einzige Iteration , in der A und B gleich viele Erfolge erzielt (wie im Durchschnitt der Fall sein wird). In jeder Iteration von $\rm simulation_\rm {underlying}$ $\rm simulation_\rm{prime}$ werden wir, gegeben A und B erstellen Zufallszahl der BetaVerteilung für jede Gruppe. Dann werden wir sie vergleichen und herausfinden, ob , was WAHR oder FALSCH (1 oder 0) ergibt. Am Ende eines Laufs von $\rm simulation_\rm{underlying}$ ${\rm Beta}_A > {\rm Beta}_B$ $\rm simulation_\rm {underlying}$ haben wir 15000 Iterationen abgeschlossen und 15000 TRUE / FALSE-Werte. Der Durchschnitt dieser einen einzigen Wert aus dem (ungefähr normal) Stichprobenverteilung des Anteils der Fließ . ${\rm Beta}_A > {\rm Beta}_B$

Aber jetzt wird ausgewählt 5000 A und B - Werte. A und B werden selten exakt gleich sein, aber die typischen Unterschiede in der Anzahl der A- und B-Erfolge werden durch die Gesamtstichprobengröße von A und B in den Schatten gestellt. Typische As und Bs ergeben aus ihrer Stichprobenverteilung der Anteile von mehr Zugkraft , sondern auch erhalten diejenigen an den Rändern der A / B - Verteilung gezogen. $\rm simulation_\rm{prime}$ ${\rm Beta}_A > {\rm Beta}_B$

Also, was im Grunde ziehen wir über viele SIM - Läufe ist eine Kombination aus Stichprobenverteilungen von für Kombinationen von A und B (mit mehr Zügen aus den Stichprobenverteilungen aus den gemeinsamen Werten von A und B als die ungewöhnlichen Werte von A und B). Dies führt zu Gemischen normaler Verteilungen. Wenn Sie sie über eine kleine Bin-Größe kombinieren (wie es die Standardeinstellung für die verwendete Histogrammfunktion ist, die direkt in Ihrem ursprünglichen Code angegeben wurde), erhalten Sie etwas, das wie eine gleichmäßige Verteilung aussieht. ${\rm Beta}_A > {\rm Beta}_B$

Erwägen:

a = b = 0.5
N = 10
samples = [] #collects the values of S
for i in range(5000):
    assert a==b
    A = np.random.binomial(N, a); B = np.random.binomial(N, b)
    S = (beta.rvs(A+1, N-A+1, size=15000) > beta.rvs(B+1, N-B+1, size=15000)).mean() 
    samples.append(S)

P.hist(samples,1000)
P.show()

russellpierce
quelle

So there is a difference between mine and your code. I sample A and B in each loop, you sample it once and calculate S 5000 times.

Cam.Davidson.Pilon

The discrepancy lies in your calls to rbinom, which returns a vector. The subsequent call to rbeta inside replicate is vectorized, so the inner (internal) loop is using a different

A

$A$ and

B

$B$ for each of the 15000 random variables generated (wrapping around for the final 5000 since your NSIM = 10000). See ?rbeta for more. This differs from @Cam's code with has a single fixed

A

$A$ and

B

$B$ used in all 15000 random-variate calls for each of the 5000 sampling (replicate) loops.

cardinal

here's the output for those curious: imgur.com/ryvWbJO

Cam.Davidson.Pilon

The only things I'm aware of that are potentially pertinent at a conceptual level are that a) the expected distribution of results are symmetric, b) a bin size of 1 is always uniform, c) a bin size of 2 for a symmetric distribution will also always appear uniform, d) the number of possible sampling distributions that can be drawn from increases with N, e) values of S can't stack up on 0 or 1 alone because beta is undefined when there are 0 successes in either group, and f) the samples are restricted between 0 and 1.

russellpierce

Betrachtet man nur, so kann man feststellen, dass sich die Abstände zwischen den Schwerpunkten der Stichprobenverteilungen verringern, wenn sich die Schwerpunkte der Stichprobenverteilungen von 0,5 entfernen (wahrscheinlich im Zusammenhang mit dem obigen Punkt f). Dieser Effekt wirkt tendenziell der Tendenz zu hohen Beobachtungsfrequenzen für die häufiger auftretenden nahezu gleichen Erfolge in Gruppe A- und Gruppe B-Fällen entgegen. Eine mathematische Lösung, warum dies der Fall ist oder warum es zu Normalverteilungen für bestimmte Behältergrößen führen sollte, ist jedoch nicht in der Nähe meines Hoheitsgebiets.

Russellpierce

To get some intuition for what is going on, let us feel free to make $N$ very large and in so doing ignore $O(1/N)$ behavior and exploit asymptotic theorems that state both Beta and Binomial distributions become approximately Normal. (With some trouble, all this can be made rigorous.) When we do this, the result emerges from a specific relationship among the various parameters.

Because we plan to use Normal approximations we will pay attention to the expectations and variances of the variables:

As Binomial $(N, p)$ variates, $n_A$ and $n_B$ have expectations of $pN$ and variances of $p(1-p)N$ . Consequently $\alpha=n_A/N$ and $\beta=n_B/N$ have expectations of $p$ and variance $p(1-p)/N$ .
As a Beta $(n_A+1, N+1-n_A)$ variate, $P_A$ has an expectation of $(n_A+1)/(N+2)$ and a variance of $(n_A+1)(N+1-n_A) / [(N+2)^2(N+3)]$ . Approximating, we find that $P_A$ has an expectation of

$E (P_{A}) = α + O (1 / N)$ $\mathbb{E}(P_A) = \alpha+O(1/N)$
and a variance of

$Var (P_{A}) = α (1 - α) / N + O (1 / N^{2}),$ $\text{Var}(P_A) = \alpha(1-\alpha)/N + O(1/N^2),$
with similar results for $P_B$ .

Let us therefore approximate the distributions of $P_A$ and $P_B$ with Normal $(\alpha, \alpha(1-\alpha)/N)$ and Normal $(\beta,\beta(1-\beta)/N)$ distributions (where the second parameter designates the variance). The distribution of $P_A-P_B$ consequently is approximately Normal; to wit,

P_{A} - P_{B} \approx Normal (α - β, \frac{α (1 - α) + β (1 - β)}{N}) .

$P_A-P_B \approx \text{Normal}\left(\alpha-\beta, \frac{\alpha(1-\alpha) + \beta(1-\beta)}{N}\right).$

For very large $N$ , the expression $\alpha(1-\alpha) + \beta(1-\beta)$ will not vary appreciably from $p(1-p)+p(1-p)=2p(1-p)$ except with very low probability (another neglected $O(1/N)$ term). Accordingly, letting $\Phi$ be the standard normal CDF,

Pr (P_{A} > P_{B}) = Pr (P_{A} - P_{B} > 0) \approx Φ (\frac{α - β}{\sqrt{2 p (1 - p) / N}}) .

$\Pr(P_A\gt P_B) =\Pr(P_A-P_B\gt 0) \approx \Phi\left(\frac{\alpha-\beta}{\sqrt{2p(1-p)/N}}\right).$

But since $\alpha-\beta$ has zero mean and variance $2p(1-p)/N,$ $Z=\frac{\alpha-\beta}{\sqrt{2p(1-p)/N}}$ is a standard Normal variate (at least approximately). $\Phi$ is its probability integral transform; $\Phi(Z)$ is uniform.

whuber
quelle

I'mn with you up until

P_{A} - P_{B} \approx N o r m a l

$P_A - P_B \approx Normal$ ... then you go off another direction that I didn't quite follow. Is

Φ

$\Phi$ defined twice, once as the standard normal CDF and then as the probability integral transform? I'm hoping you can expand your description around these steps and relate them to the initial code/problem. Maybe loop back around and restate which specific parameters produce the uniform result.

russellpierce

@rpierce (1) The difference

P_{A} - P_{B}

$P_A-P_B$ is approximately normal because

P_{A}

$P_A$ and

P_{B}

$P_B$ are independent and each is approximately normal. The mean is the difference of the means and the variance is the sum of the variances. (2) The probability integral transform is the CDF: it is the case for any random variable

X

$X$ with continuous distribution

F

$F$ , that

F (X)

$F(X)$ is uniform.

whuber

Oh I got 1, it was the stuff after it where I got lost. This will be mindbogglingly dumb, but why is

P r (P_{A} > P_{B})

$Pr(P_A>P_B)$ the same as the CDF?

russellpierce

@rpierce That follows rather directly from the definition, but there's a slight twist in which the symmetry of the Normal distribution is invoked. We're dealing with a Normal variate

X = P_{A} - P_{B}

$X = P_A-P_B$ assumed to have an expectation of

μ = α - β

$\mu=\alpha-\beta$ and variance

σ^{2} = 2 p (1 - p) / N

$\sigma^2 = 2p(1-p)/N$ . Standardizing

X

$X$ , it is natural to rewrite the probability as

Pr (X > 0) = Pr ((X - μ) / σ > (0 - μ) / σ) = 1 - Φ (- μ / σ) = Φ (μ / σ) .

$\Pr(X\gt 0) = \Pr((X-\mu)/\sigma \gt (0-\mu)/\sigma) = 1-\Phi(-\mu/\sigma) = \Phi(\mu/\sigma).$

whuber

@whuber this is pretty amazing. You are a wonderful teacher. I appreciate both yours and rpierce's answer, I'll still give him credit as it did solve our problem, and you have shown why the behaviour occurs. Ty!

Cam.Davidson.Pilon