Kami sedang menyelidiki uji statistik Bayesian, dan menemukan fenomena aneh (bagi saya setidaknya).
Pertimbangkan kasus berikut: kami tertarik untuk mengukur populasi mana, A atau B, yang memiliki tingkat konversi yang lebih tinggi. Untuk pemeriksaan kewarasan, kami menetapkan , yaitu probabilitas konversi sama di kedua grup. Kami menghasilkan data buatan menggunakan model binomial, misalnya
Kami kemudian mencoba memperkirakan menggunakan model beta-binomial Bayesian sehingga kami mendapatkan posisi untuk setiap tingkat konversi, mis.
Statistik pengujian kami dihitung dengan menghitung melalui monte carlo.
Yang mengejutkan saya adalah jika , maka S \ sim \ text {Uniform (0,1)} . Pikiranku adalah bahwa itu akan berpusat di sekitar 0,5, dan bahkan bertemu menjadi 0,5 ketika ukuran sampel, N , bertambah.
Pertanyaan saya adalah, mengapa saat ?
Berikut beberapa kode Python untuk diperagakan:
%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()
sumber
R
Jawaban:
TL; DR: Campuran dari distribusi normal mungkin terlihat seragam ketika ukuran bin besar.
Jawaban ini meminjam dari kode sampel @ whuber (yang saya pikir pertama adalah kesalahan, tetapi dalam retrospeksi mungkin merupakan petunjuk).
Proporsi yang mendasari dalam populasi adalah sama:
a = b = 0.5
.Masing-masing kelompok, A dan B memiliki 10.000 anggota:
N = 10000
.Kita akan melakukan 5000 ulangan dari simulasi:
for i in range(5000):
.Sebenarnya, apa yang kita lakukan adalah dari . Di masing-masing dari 5000 iterasi kita akan melakukan . s i msimulationprime s i m u l a t i o n p r i m e s i m u l a t i o n usimulationunderlying simulationprime simulationunderlying
Dalam setiap iterasi dari kita akan mensimulasikan nomor acak dari A dan B yang 'keberhasilan' (AKA dikonversi) diberikan proporsi yang mendasari yang sama ditetapkan sebelumnya: . Secara nominal ini akan menghasilkan A = 5000 dan B = 5000, tetapi A dan B bervariasi dari sim run ke sim run dan didistribusikan melintasi 5000 simulasi berjalan secara independen dan (kurang-lebih) secara normal (kami akan kembali ke sana).simulationprime
A = np.random.binomial(N, a); B = np.random.binomial(N, b)
Sekarang mari kita melangkahi untuk satu iterasi tunggal di mana A dan B telah mengambil jumlah keberhasilan yang sama (seperti rata-rata kasus). Dalam setiap iterasi dari kita akan, diberikan A dan B, membuat variasi acak dari distribusi beta untuk setiap grup. Kemudian kita akan membandingkannya dan mencari tahu apakah , menghasilkan TRUE atau FALSE (1 atau 0). Pada akhir jangka waktu , kami telah menyelesaikan 15000 iterasi dan memiliki 15000 nilai TRUE / FALSE. Rata-rata ini akan menghasilkan nilai tunggal dari distribusi sampling (sekitar normal) dari proporsi s i m u l a t i o n p r i m e s i m u l a t i o n u n d e r l y i n g B e t a A >simulationunderlying simulationprime simulationunderlying simulatio n u n d e r l y i n g B e t a A > B e t a BBetaA>BetaB simulationunderlying BetaA>BetaB .
Kecuali sekarang akan memilih nilai 5000 A dan B. A dan B jarang akan sama persis, tetapi perbedaan tipikal dalam jumlah keberhasilan A dan B dikerdilkan oleh ukuran sampel total A dan B. As dan B khas akan menghasilkan lebih banyak tarikan dari distribusi sampling mereka dengan proporsi , tetapi yang ada di tepi distribusi A / B juga akan ditarik.B e t a A > B e t a Bsimulationprime BetaA>BetaB
Jadi, yang pada intinya kita menarik banyak sim run adalah kombinasi distribusi sampling dari untuk kombinasi A dan B (dengan lebih banyak tarikan dari distribusi sampling yang dibuat dari nilai-nilai umum A dan B dari nilai-nilai yang tidak biasa dari A dan B). Ini menghasilkan campuran distribusi normal-ish. Ketika Anda menggabungkan mereka di atas ukuran nampan kecil (seperti default untuk fungsi histogram yang Anda gunakan dan ditentukan langsung dalam kode asli Anda), Anda berakhir dengan sesuatu yang tampak seperti distribusi yang seragam.BetaA>BetaB
Mempertimbangkan:
sumber
rbinom
, yang mengembalikan vektor. Panggilan berikutnya kerbeta
dalamreplicate
adalah vektor, sehingga loop (internal) internal menggunakan dan berbeda untuk masing-masing 15000 variabel acak yang dihasilkan (membungkus sekitar 5000 final sejak Anda ). Lihat lebih lanjut. Ini berbeda dari kode @ Cam dengan memiliki dan tetap tunggal yang digunakan di semua 15000 panggilan variate acak untuk masing-masing dari 5000 sampling ( ) loop. B A BNSIM = 10000
?rbeta
replicate
Untuk mendapatkan intuisi tentang apa yang sedang terjadi, mari kita merasa bebas untuk membuat sangat besar dan dengan demikian mengabaikan perilaku dan mengeksploitasi teorema asimtotik yang menyatakan bahwa distribusi Beta dan Binomial menjadi mendekati Normal. (Dengan beberapa masalah, semua ini bisa dibuat sulit.) Ketika kita melakukan ini, hasilnya muncul dari hubungan spesifik di antara berbagai parameter.O ( 1 / N )N O(1/N)
Karena kami berencana untuk menggunakan perkiraan Normal, kami akan memperhatikan harapan dan varian variabel:
Sebagai Binomial variates, dan memiliki harapan dan varians dari . Akibatnya dan memiliki harapan dari dan varians .n A n B p N p ( 1 - p ) N α = n A / N β = n B / N p p ( 1 - p ) / N(N,p) nA nB pN p(1−p)N α=nA/N β=nB/N p p(1−p)/N
Sebagai Beta , memiliki ekspektasi dan varian . Mendekati, kami menemukan bahwa memiliki harapanP A ( n A + 1 ) / ( N + 2 ) ( n A + 1 ) ( N + 1 - n A ) / [ ( N + 2 ) 2 ( N + 3 ) ](nA+1,N+1−nA) PA (nA+1)/(N+2) (nA+1)(N+1−nA)/[(N+2)2(N+3)] PA
dan varian dari
dengan hasil yang serupa untuk .PB
Oleh karena itu, mari kita distribusi dan dengan Normal dan Normal (di mana parameter kedua menentukan varians ). Distribusi akibatnya kira-kira Normal; yakni,PA PB (α,α(1−α)/N) (β,β(1−β)/N) PA−PB
Untuk sangat besar , ekspresi tidak akan berbeda dari kecuali dengan probabilitas yang sangat rendah ( istilah lainnya yang diabaikan ). Dengan demikian, membiarkan menjadi CDF normal standar,N α(1−α)+β(1−β) p(1−p)+p(1−p)=2p(1−p) O(1/N) Φ
Tetapi karena memiliki rata-rata nol dan varians adalah standar Normal variate (setidaknya sekitar). adalah probabilitas integral perubahannya ; adalah seragam .α−β 2p(1−p)/N, Z=α−β2p(1−p)/N√ Φ Φ(Z)
sumber