Warum ist diese Verteilung einheitlich?

12

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 pA=pB , dh die Konversionswahrscheinlichkeit ist in beiden Gruppen gleich. Wir erzeugen künstliche Daten mit einem Binomialmodell, zB

nABinomial(N,pA)

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

PABeta(1+nA,NnA+1)

Unsere Teststatistik wird durch Berechnung von S = P ( P A > P B berechnetS=P(PA>PB|N,nA,nB) über Monte Carlo.

Was mich überraschte, war, dass wenn pA=pB , dann SUniform(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 SUniform(0,1) wenn pA=pB ?


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()
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 N (weniger als 100 / min ( p , 1 - p ) , ungefähr mit p = p A = p B ) nicht einmal annähernd gleichförmig. SN100/min(p,1p)p=pA=pB
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
2
Ich muss also dein Setup missverstehen. Soweit ich das beurteilen kann , ist der Wert von S für alle gegebenen Werte von eine Zahl. Daher Annahme , dass N , P A , und p B sind für den Moment festgelegt (wie sie im Code sind), S ist eine Funktion der ( n A , n B ) . Letzteres kann jedoch als Realisierung von zwei Binomialverteilungen nur eine diskrete Menge von Werten erreichen. Wenn ich Ihren Code wiedergebeN,nA,nB,SN,pA,pBS(nA,nB)RIch entschieden nicht-einheitliche Histogramme für kleine bekommen . N
whuber
1
Obwohl Ihr tatsächlich Werte zwischen 0 und 1 hat , verwechseln Sie das nicht mit nicht-diskret: Es kann höchstens N 2 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 S anstelle der korrekten Werte generiert und die Schätzungen im Wesentlichen eine kontinuierliche Verteilung aufweisen. S01N2S
whuber
1
@whuber ja, du hast recht, ausgezeichnete Beobachtung. Ich bin immer noch gespannt, warum es dann einheitlich aussieht .
Cam.Davidson.Pilon

Antworten:

11

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 s i m u l a t i o n u n d e r l y i n g . In jedem der 5000 Iterationen s ich bin u l a t i o n p r i m e wir tun isimulationprimesimulationunderlyingsimulationprime .simulationunderlying

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).simulationprimeA = 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 s i m u l a t i o n p r i m e , in der A und B gleich viele Erfolge erzielt (wie im Durchschnitt der Fall sein wird). In jeder Iteration von s i m u l a t i o n usimulationunderlyingsimulationprime werden wir, gegeben A und B erstellen Zufallszahl der BetaVerteilung für jede Gruppe. Dann werden wir sie vergleichen und herausfinden, ob B e t a A > B e t a B ist , was WAHR oder FALSCH (1 oder 0) ergibt. Am Ende eines Laufs vonsimuleinemtio n u n d e r l y i n gsimulationunderlyingBetaA>BetaBsimulationunderlyinghaben 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ß .BetaA>BetaB

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 B e mehr Zugkraft t a A > B e t ein B , sondern auch erhalten diejenigen an den Rändern der A / B - Verteilung gezogen.simulationprimeBetaA>BetaB

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.BetaA>BetaB

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
1
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
1
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 and 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 and B used in all 15000 random-variate calls for each of the 5000 sampling (replicate) loops.
cardinal
1
here's the output for those curious: imgur.com/ryvWbJO
Cam.Davidson.Pilon
1
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
1
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
16

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, nA and nB have expectations of pN and variances of p(1p)N. Consequently α=nA/N and β=nB/N have expectations of p and variance p(1p)/N.

  • As a Beta(nA+1,N+1nA) variate, PA has an expectation of (nA+1)/(N+2) and a variance of (nA+1)(N+1nA)/[(N+2)2(N+3)]. Approximating, we find that PA has an expectation of

    E(PA)=α+O(1/N)

    and a variance of

    Var(PA)=α(1α)/N+O(1/N2),

    with similar results for PB.

Let us therefore approximate the distributions of PA and PB with Normal(α,α(1α)/N) and Normal(β,β(1β)/N) distributions (where the second parameter designates the variance). The distribution of PAPB consequently is approximately Normal; to wit,

PAPBNormal(αβ,α(1α)+β(1β)N).

For very large N, the expression α(1α)+β(1β) will not vary appreciably from p(1p)+p(1p)=2p(1p) except with very low probability (another neglected O(1/N) term). Accordingly, letting Φ be the standard normal CDF,

Pr(PA>PB)=Pr(PAPB>0)Φ(αβ2p(1p)/N).

But since αβ has zero mean and variance 2p(1p)/N, Z=αβ2p(1p)/N is a standard Normal variate (at least approximately). Φ is its probability integral transform; Φ(Z) is uniform.

whuber
quelle
1
I'mn with you up until PAPBNormal... then you go off another direction that I didn't quite follow. Is Φ 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
1
@rpierce (1) The difference PAPB is approximately normal because PA and PB 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 with continuous distribution F, that F(X) is uniform.
whuber
1
Oh I got 1, it was the stuff after it where I got lost. This will be mindbogglingly dumb, but why is Pr(PA>PB) the same as the CDF?
russellpierce
1
@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=PAPB assumed to have an expectation of μ=αβ and variance σ2=2p(1p)/N. Standardizing X, it is natural to rewrite the probability as
Pr(X>0)=Pr((Xμ)/σ>(0μ)/σ)=1Φ(μ/σ)=Φ(μ/σ).
whuber
3
@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