Model pas untuk dua distribusi normal di PyMC

10

Karena saya seorang insinyur perangkat lunak yang mencoba mempelajari lebih banyak statistik, Anda harus memaafkan saya bahkan sebelum saya mulai, ini adalah wilayah newb yang serius ...

Saya telah belajar PyMC dan bekerja melalui beberapa contoh sederhana. Satu masalah yang saya tidak dapat mulai bekerja (dan tidak dapat menemukan contoh terkait) adalah memasang model ke data yang dihasilkan dari dua distribusi normal.

Katakanlah saya memiliki 1000 nilai; 500 dihasilkan dari a Normal(mean=100, stddev=20)dan 500 lainnya dihasilkan dari a Normal(mean=200, stddev=20).

Jika saya ingin mencocokkan model dengan mereka, yaitu menentukan dua cara dan standar deviasi tunggal, menggunakan PyMC. Aku tahu itu sesuatu di sepanjang garis ...

mean1 = Uniform('mean1', lower=0.0, upper=200.0)
mean2 = Uniform('mean2', lower=0.0, upper=200.0)
precision = Gamma('precision', alpha=0.1, beta=0.1)

data = read_data_from_file_or_whatever()

@deterministic(plot=False)
def mean(m1=mean1, m2=mean2):
    # but what goes here?

process = Normal('process', mu=mean, tau=precision, value=data, observed=True)

yaitu, proses menghasilkan adalah Normal, tetapi mu adalah salah satu dari dua nilai. Saya hanya tidak tahu bagaimana mewakili "keputusan" antara apakah suatu nilai berasal m1atau tidak m2.

Mungkin saya benar-benar mengambil pendekatan yang salah untuk pemodelan ini? Adakah yang bisa menunjukkan saya pada contoh? Saya dapat membaca BUGS dan JAGS sehingga semuanya benar-benar baik-baik saja.

mat kelcey
sumber

Jawaban:

11

Apakah Anda benar-benar yakin bahwa separuh berasal dari satu distribusi dan setengah lainnya dari yang lain? Jika tidak, kita dapat memodelkan proporsi sebagai variabel acak (yang merupakan hal yang sangat bayesian untuk dilakukan).

Berikut ini yang akan saya lakukan, beberapa tips tertanam.

from pymc import *

size = 10
p = Uniform( "p", 0 , 1) #this is the fraction that come from mean1 vs mean2

ber = Bernoulli( "ber", p = p, size = size) # produces 1 with proportion p.

precision = Gamma('precision', alpha=0.1, beta=0.1)

mean1 = Normal( "mean1", 0, 0.001 ) #better to use normals versus Uniforms (unless you are certain the value is  truncated at 0 and 200 
mean2 = Normal( "mean2", 0, 0.001 )

@deterministic
def mean( ber = ber, mean1 = mean1, mean2 = mean2):
    return ber*mean1 + (1-ber)*mean2


#generate some artificial data   
v = np.random.randint( 0, 2, size)
data = v*(10+ np.random.randn(size) ) + (1-v)*(-10 + np.random.randn(size ) )


obs = Normal( "obs", mean, precision, value = data, observed = True)

model = Model( {"p":p, "precision": precision, "mean1": mean1, "mean2":mean2, "obs":obs} )
Cam.Davidson.Pilon
sumber
2
Promosi memalukan: Saya baru saja menulis artikel blog tentang Bayes dan pyMC secara harfiah 1 menit sebelum Anda memposting ini, jadi saya mengundang Anda untuk memeriksanya. The Awesome Power of Bayes - Bagian 1
Cam.Davidson.Pilon
luar biasa! ini pendekatan untuk pencampuran dari dua cara persis apa yang saya coba untuk mendapatkan kepalaku.
mat kelcey
Tidak yakin saya sepenuhnya memahami manfaat pemodelan sebenarnya dari mengatakan mean1 & mean2 didistribusikan secara normal, bukan Uniform (Sama berlaku untuk ketepatan jujur, saya telah menggunakan Gamma sejak "orang lain melakukannya"). Saya harus banyak belajar :)
mat kelcey
Menggunakan Seragam, seperti dalam contoh asli Anda, menyiratkan bahwa Anda tahu dengan kepastian absolut bahwa rata-rata tidak melebihi nilai tertentu. Ini agak patologis. Lebih baik menggunakan normal, karena memungkinkan semua bilangan real untuk dipertimbangkan.
Cam.Davidson.Pilon
1
Pilihan gamma memiliki alasan matematis. Gamma adalah konjugat sebelum ketepatan, lihat tabel di sini
Cam.Davidson.Pilon
6

Beberapa poin, terkait dengan diskusi di atas:

  1. Pilihan difus normal vs seragam cukup akademis kecuali (a) Anda khawatir tentang konjugasi, dalam hal ini Anda akan menggunakan normal atau (b) ada kemungkinan yang masuk akal bahwa nilai sebenarnya bisa berada di luar titik akhir seragam. . Dengan PyMC, tidak ada alasan untuk khawatir tentang konjugasi, kecuali jika Anda secara khusus ingin menggunakan sampler Gibbs.

  2. Gamma sebenarnya bukan pilihan tepat untuk informasi sebelum parameter varians / presisi. Itu bisa menjadi lebih informatif menurut Anda. Pilihan yang lebih baik adalah meletakkan seragam sebelum deviasi standar, kemudian mengubahnya dengan kuadrat terbalik. Lihat Gelman 2006 untuk detailnya.

fonnesbeck
sumber
1
ah fonnesbeck adalah salah satu pengembang inti pymc! Bisakah Anda menunjukkan kepada kami contoh cara membuat kode titik 2?
Cam.Davidson.Pilon
terima kasih fonnesbeck dan, ya tolong! ke eg cepat dari poin 2 :)
mat kelcey
1
sebenarnya saya kira Anda berarti sesuatu di sepanjang baris ... gist.github.com/4404631 ?
mat kelcey
Ya persis. Anda dapat melakukan transformasi sedikit lebih ringkas:tau = std_dev**-2
fonnesbeck
apa tempat yang tepat untuk membaca tentang dari mana hubungan antara presisi dan std_dev ini berasal?
user979