Dlaczego ta dystrybucja jest jednolita?

12

Badamy bayesowskie testy statystyczne i natrafiliśmy na dziwne (przynajmniej dla mnie) zjawisko.

Rozważ następujący przypadek: interesuje nas pomiar, która populacja, A lub B, ma wyższy współczynnik konwersji. Dla kontroli poczytalności ustawiamy , to znaczy prawdopodobieństwo konwersji jest równe w obu grupach. Generujemy sztuczne dane przy użyciu modelu dwumianowego, np.pA=pB

nABinomial(N,pA)

Następnie próbujemy oszacować za pomocą bayesowskiego modelu dwumianowego beta, aby uzyskać tylne współczynniki dla każdego współczynnika konwersji, np.pA,pB

PABeta(1+nA,NnA+1)

Nasza statystyka testowa jest obliczana przez obliczenie za pomocą Monte Carlo.S=P(PA>PB|N,nA,nB)

Zaskoczyło mnie to, że jeśli , to . Myślałem, że będzie on wyśrodkowany wokół 0,5, a nawet zbiegnie się do 0,5 w miarę wzrostu wielkości próbki, ,. pA=pBSUniform(0,1)N

Moje pytanie brzmi: dlaczego kiedy ?SUniform(0,1)pA=pB


Oto kod Python do zademonstrowania:

%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
źródło
Zauważ, że nie może być dokładnie jednorodne, ponieważ jest zmienną dyskretną. Pytasz zatem o zachowanie asymptotyczne. Co więcej, dla małych (mniej niż , w przybliżeniu, przy ) rozkład nie jest nawet w przybliżeniu zbliżony do jednorodności. SN100/min(p,1p)p=pA=pB
whuber
@ Whuber S nie jest dyskretny, jest prawdopodobne, że może on spaść między 0 a 1. Również, nawet dla niskiej N, obserwuję jednolite zachowanie.
Cam.Davidson.Pilon
2
Więc chyba nie rozumiem twojej konfiguracji. O ile mogę stwierdzić, dla każdej podanej wartości wartość jest liczbą. Dlatego, akceptując, że i są ustalone (tak jak są w kodzie), jest funkcją . Ale ta ostatnia, będąca realizacją dwóch rozkładów dwumianowych, może osiągnąć jedynie dyskretny zestaw wartości. Kiedy reprodukować swój kod , mam zdecydowanie niejednolitych histogramy dla małych . N,nA,nB,SN,pA,pBS(nA,nB)RN
whuber
1
Chociaż faktycznie twoje ma wartości od do , nie myl tego z niedyskretnymi: może mieć co najwyżej odrębne wartości (i faktycznie ma mniej niż to). To może nie być całkowicie jasne dla ciebie, ponieważ symulacja generuje szacunki dotyczące , a nie jego poprawne wartości i szacunki zasadniczo mają ciągłą dystrybucję. S01N2S
whuber
1
@ whuber tak, masz rację, doskonała obserwacja. Nadal utknąłem na tym, dlaczego to wygląda jednolicie.
Cam.Davidson.Pilon

Odpowiedzi:

11

TL; DR: Mieszaniny normalnych rozkładów mogą wyglądać jednolicie, gdy rozmiary pojemników są duże.

Ta odpowiedź zapożycza z przykładowego kodu @ Whuber (który, jak sądzę, był początkowo błędem, ale z perspektywy czasu prawdopodobnie był wskazówką).

Bazowe proporcje w populacji są równe: a = b = 0.5.
Każda z grup A i B ma 10000 członków: N = 10000.
Mamy zamiar przeprowadzić 5000 powtórzeń symulacji: for i in range(5000):.

Właściwie to, co robimy, to z . W każdej z 5000 iteracji zrobimy . s i m u l a t i o n u n d e r l y i n g s i m u l a t i o n p r i m e s I m u ł a t i o n usimulationprimesimulationunderlyingsimulationprimesimulationunderlying

W każdej iteracji my symulować liczbę losową A i B, które są „osiągnięcia” (czyli przekształcone) podane równe proporcje bazowe zdefiniowane wcześniej: . Nominalnie da to A = 5000 i B = 5000, ale A i B różnią się w zależności od przebiegu karty SIM i są rozdzielone na 5000 przebiegów symulacji niezależnie i (w przybliżeniu) normalnie (do tego wrócimy).simulationprimeA = np.random.binomial(N, a); B = np.random.binomial(N, b)

Przejdźmy teraz przez dla pojedynczej iteracji w której A i B odniosły taką samą liczbę sukcesów (co będzie średnią w przypadku). W każdej iteracji będziemy, biorąc pod uwagę A i B, tworzyć losowe warianty rozkładu beta dla każdej grupy. Następnie porównamy je i dowiemy się, czy , dając wartość PRAWDA lub FAŁSZ (1 lub 0). Pod koniec serii zakończyliśmy 15000 iteracji i otrzymaliśmy 15000 wartości PRAWDA / FAŁSZ. Średnia z nich da pojedynczą wartość z (w przybliżeniu normalnego) rozkładu próbkowania proporcji s i m U L t i o n s r i m e s I m u ł a t i o n u n d e r l y i n g B E t a A >simulationunderlyingsimulationprimesjamulzatjaonunremirlyjansol simulatio n u n d e r l y i n g B E t a A > B e t a BBetaA>BetaBsimulationunderlyingBetaA>BetaB .

Z wyjątkiem teraz wybierze 5000 wartości A i B. A i B rzadko są dokładnie równe, ale typowe różnice w liczbie sukcesów A i B są zmniejszone przez całkowitą wielkość próby A i B. Typowe As i Bs przyniosą więcej ściągnięć z ich rozkładu próbkowania proporcji , ale te na krawędziach dystrybucji A / B również zostaną pociągnięte.B e t a A > B e t a BsimulationprimeBetaA>BetaB

Tak więc, w zasadzie wyciągamy wiele przebiegów SIM, jest to kombinacja rozkładów próbkowania dla kombinacji A i B (z większą ilością ściągnięć z rozkładów próbkowania wykonanych ze wspólnych wartości A i B niż nietypowe wartości A i B). Daje to w wyniku mieszaniny rozkładów normalnych. Kiedy połączysz je w małym pojemniku (co jest domyślną funkcją histogramu, której użyłeś i która została określona bezpośrednio w oryginalnym kodzie), otrzymujesz coś, co wygląda jak jednolity rozkład.BetaA>BetaB

Rozważać:

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
źródło
1
Jest więc różnica między moim a twoim kodem. Próbuję A i B w każdej pętli, próbkujesz je raz i obliczasz S 5000 razy.
Cam.Davidson.Pilon
1
Rozbieżność polega na twoich połączeniach z rbinom, które zwracają wektor. Późniejsze wywołanie do rbetawewnątrz replicatejest wektoryzowane, więc wewnętrzna (wewnętrzna) pętla używa różnych i B dla każdej z 15000 wygenerowanych zmiennych losowych (owija się wokół końcowych 5000 od twojego ). Zobacz więcej. Różni się od kodu @ Cam tym, że ma jeden stały A i B używany we wszystkich 15000 wywołań losowych dla każdej z 5000 pętli próbkowania ( ). ABNSIM = 10000?rbetaABreplicate
kardynał
1
oto wyniki dla tych ciekawych: imgur.com/ryvWbJO
Cam.Davidson.Pilon
1
Jedyne, o czym wiem, że są potencjalnie istotne na poziomie pojęciowym, to: a) oczekiwany rozkład wyników jest symetryczny, b) rozmiar pojemnika 1 jest zawsze równomierny, c) rozmiar pojemnika 2 dla rozkładu symetrycznego będzie również zawsze wyglądał jednolicie, d) liczba możliwych rozkładów próbkowania, które można wyciągnąć ze wzrostu przy N, e) wartości S nie mogą się kumulować na samym 0 lub 1, ponieważ beta jest niezdefiniowana, gdy w każdej grupie jest 0 sukcesów oraz f) próbki są ograniczone od 0 do 1.
russellpierce
1
Na podstawie samej obserwacji możemy zauważyć, że odległości między centroidami rozkładów próbkowania zmniejszają się, gdy centroidy rozkładów próbkowania oddalają się od 0,5 (prawdopodobnie związane z punktem f powyżej). Efekt ten ma tendencję do przeciwdziałania tendencji do wysokich częstotliwości obserwacji częstszych, niemal równych sukcesów w przypadku grupy A i grupy B. Jednak podanie matematycznego rozwiązania, dlaczego tak się dzieje lub dlaczego powinno dawać normalne rozkłady dla niektórych rozmiarów pojemników, nie jest w pobliżu mojego terytorium.
russellpierce
16

Aby uzyskać intuicję w tym, co się dzieje, nie krępujmy się, aby bardzo duży, a tym samym ignorując zachowanie O ( 1 / N ) i wykorzystując asymptotyczne twierdzenia, które stwierdzają, że zarówno rozkład beta, jak i dwumianowy stają się w przybliżeniu normalne. (Przy pewnym problemie wszystko to można uczynić rygorystycznym.) Kiedy to robimy, wynik wyłania się z określonej zależności między różnymi parametrami.NO(1/N)


Ponieważ planujemy stosować przybliżenia normalne, zwrócimy uwagę na oczekiwania i wariancje zmiennych:

  • Jako dwumianowego zmiennych towarzyszących, n i n B mają oczekiwania p N a wariancje p ( 1 - p ) N . W związku z tym α = n / N i β = N B / N mieć oczekiwania p i wariancji p ( 1 - p ) / N .(N,p)nAnBpNp(1p)Nα=nA/Nβ=nB/Npp(1p)/N

  • Ponieważ zmienia się Beta , P A ma oczekiwane ( n A + 1 ) / ( N + 2 ) i wariancję ( n A + 1 ) ( N + 1 - n A ) / [ ( N + 2 ) 2 ( N + 3)(nA+1,N+1nA)PA(nA+1)/(N+2) . Przybliżając, stwierdzamy, że P A ma oczekiwania(nA+1)(N+1nA)/[(N+2)2(N+3)]PA

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

    i wariant

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

    z podobnymi wynikami dla .PB

Niech się więc w przybliżeniu rozkładów i P B z normalnymi ( alfa , alfa ( 1 - alfa ) / N ) i Normal ( β , β ( 1 - P ) / N ) dystrybucji (gdzie drugi parametr wyznacza odchylenie ) . W związku z tym rozkład P A - P B jest w przybliżeniu normalny; dowcipPAPB(α,α(1α)/N)(β,β(1β)/N)PAPB

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

W przypadku bardzo dużego wyrażenie α ( 1 - α ) + β ( 1 - β ) nie różni się znacznie od p ( 1 - p ) + p ( 1 - p ) = 2 p ( 1 - p ), z wyjątkiem bardzo niskiej prawdopodobieństwo (inny niedoceniany termin O ( 1 / N ) ). Odpowiednio, pozwalając Φ być standardowym normalnym CDF,Nα(1α)+β(1β)p(1p)+p(1p)=2p(1p)O(1/N)Φ

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

Ale ponieważ ma zerową średnią i wariancję 2 p ( 1 - p ) / N , Z = α - βαβ2p(1p)/N, Z=αβ2p(1p)/NΦΦ(Z)

Whuber
źródło
1
PAPBNormalΦ
1
PAPBPAPBXFF(X)
1
Pr(PA>PB)
1
X=PAPBμ=αβσ2=2p(1p)/NX
Pr(X>0)=Pr((Xμ)/σ>(0μ)/σ)=1Φ(μ/σ)=Φ(μ/σ).
3
@ whuber to jest całkiem niesamowite. Jesteś wspaniałym nauczycielem. Doceniam odpowiedź zarówno twoją, jak i rpierce'a. Nadal mu to doceniam, ponieważ rozwiązało to nasz problem, a ty pokazałeś, dlaczego tak się dzieje. Ty!
Cam.Davidson.Pilon