Dlaczego ta dystrybucja jest jednolita?

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. $p_A = p_B$

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

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

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

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

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

Nasza statystyka testowa jest obliczana przez obliczenie za pomocą Monte Carlo. $S = P(P_A > P_B\; |\; N, n_A, n_B)$

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, ,. $p_A = p_B$ $S \sim \text{Uniform(0,1)}$ $N$

Moje pytanie brzmi: dlaczego kiedy ? $S \sim \text{Uniform(0,1)}$ $p_A = p_B$

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.

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 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

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, 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)$ R

N

$N$

— whuber

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ę.

S

$S$

0

$0$

1

$1$

N^{2}

$N^2$

S

$S$

— whuber

@ whuber tak, masz rację, doskonała obserwacja. Nadal utknąłem na tym, dlaczego to wygląda jednolicie.

— Cam.Davidson.Pilon

Odpowiedzi:

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 . $\rm simulation_\rm{prime}$ $\rm simulation_\rm{underlying}$ $\rm simulation_\rm{prime}$ $\rm simulation_\rm{underlying}$

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). $\rm simulation_\rm{prime}$ A = 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 $\rm simulation_\rm {underlying}$ $\rm simulation_\rm{prime}$ $\rm simulation_\rm{underlying}$ ${\rm Beta}_A > {\rm Beta}_B$ $\rm simulation_\rm {underlying}$ ${\rm Beta}_A > {\rm Beta}_B$ .

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. $\rm simulation_\rm{prime}$ ${\rm Beta}_A > {\rm Beta}_B$

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. ${\rm Beta}_A > {\rm Beta}_B$

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

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

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

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

używany we wszystkich 15000 wywołań losowych dla każdej z 5000 pętli próbkowania ( ).

A

$A$

B

$B$ NSIM = 10000?rbeta

A

$A$

B

$B$ replicate

— kardynał

oto wyniki dla tych ciekawych: imgur.com/ryvWbJO

— Cam.Davidson.Pilon

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

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

Aby uzyskać intuicję w tym, co się dzieje, nie krępujmy się, aby bardzo duży, a tym samym ignorując zachowanie 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. $N$ $O(1/N)$

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

Jako dwumianowego zmiennych towarzyszących, i mają oczekiwania a wariancje . W związku z tym i mieć oczekiwania i wariancji . $(N, p)$ $n_A$ $n_B$ $pN$ $p(1-p)N$ $\alpha=n_A/N$ $\beta=n_B/N$ $p$ $p(1-p)/N$
Ponieważ zmienia się Beta , ma oczekiwane i wariancję $(n_A+1, N+1-n_A)$ $P_A$ $(n_A+1)/(N+2)$ . Przybliżając, stwierdzamy, że ma oczekiwania $(n_A+1)(N+1-n_A) / [(N+2)^2(N+3)]$ $P_A$

$E (P_{A}) = α + O (1 / N)$ $\mathbb{E}(P_A) = \alpha+O(1/N)$
i wariant

$Var (P_{A}) = α (1 - α) / N + O (1 / N^{2}),$ $\text{Var}(P_A) = \alpha(1-\alpha)/N + O(1/N^2),$
z podobnymi wynikami dla . $P_B$

Niech się więc w przybliżeniu rozkładów i z normalnymi i Normal dystrybucji (gdzie drugi parametr wyznacza odchylenie ) . W związku z tym rozkład jest w przybliżeniu normalny; dowcip $P_A$ $P_B$ $(\alpha, \alpha(1-\alpha)/N)$ $(\beta,\beta(1-\beta)/N)$ $P_A-P_B$

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).$

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

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).$

Ale ponieważ ma zerową średnią i wariancję $\alpha-\beta$ $2p(1-p)/N,$ $Z=\frac{\alpha-\beta}{\sqrt{2p(1-p)/N}}$ $\Phi$ $\Phi(Z)$

— Whuber
źródło

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

$P_A - P_B \approx Normal$

Φ

$\Phi$

P_{A} - P_{B}

$P_A-P_B$

P_{A}

$P_A$

P_{B}

$P_B$

X

$X$

F

$F$

F (X)

$F(X)$

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

$Pr(P_A>P_B)$

X = P_{A} - P_{B}

$X = P_A-P_B$

μ = α - β

$\mu=\alpha-\beta$

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

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

X

$X$

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 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