PyMC 中兩個正態分佈的擬合模型

由於我是一名試圖了解更多統計數據的軟件工程師，所以在我開始之前你必須原諒我，這是嚴重的新手領域……

我一直在學習PyMC並研究了一些非常（非常）簡單的例子。我無法開始工作（也找不到任何相關示例）的一個問題是將模型擬合到從兩個正態分佈生成的數據。

假設我有 1000 個值；500 從 a 生成Normal(mean=100, stddev=20)，另外 500 從 a 生成Normal(mean=200, stddev=20)。

如果我想為它們擬合模型，即使用 PyMC 確定兩個均值和單個標準差。我知道這有點像……

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)

即，生成過程是正常的，但 mu 是兩個值之一。我只是不知道如何表示一個值是否來自m1or之間的“決定” m2。

也許我只是完全採取了錯誤的方法來建模這個？誰能給我舉個例子？我可以閱讀 BUGS 和 JAGS，所以一切都很好。

您絕對確定一半來自一個發行版，另一半來自另一個發行版嗎？如果不是，我們可以將比例建模為隨機變量（這是一個非常貝葉斯的事情）。

以下是我會做的，嵌入了一些提示。

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

引用自：https://stats.stackexchange.com/questions/46626