PyMcの最初の使用は失敗します

Question

私はPyMcを初めて使用しています。なぜこのコードが機能しないのかを知りたいと思います。私はすでにこれに数時間を費やしましたが、何かが恋しいです。誰か助けてくれますか？私が対処したい

質問：

• 私は3個のバンプを示しNPTS対策のセットを持っているので、私は3つのガウス分布の和としてこれをモデル化したい（対策を想定してうるさいとガウスです約8つのパラメータ：バンプの相対的な重み（すなわち2つのパラメータ）、3つの平均と3つの分散を見積もります。

• 私は、このアプローチが同じバンプを持たない可能性のある他のセットにも適用できるようにしたいので、フラットなプライヤーを取ります。

問題： 以下のコードは、私には間違った推測を与えます。どうしましたか ？ thx

""" 
hypothesis: multimodal distrib sum of 3 gaussian distributions 

model description: 
* p1, p2, p3 are the probabilities for a point to belong to gaussian 1, 2 or 3 
==> p1, p2, p3 are the relative weights of the 3 gaussians 

* once a point is associated with a gaussian, 
it is distributed normally according to the parameters mu_i, sigma_i of the gaussian 
but instead of considering sigma, pymc prefers considering tau=1/sigma**2 

* thus, PyMc must guess 8 parameters: p1, p2, mu1, mu2, mu3, tau1, tau2, tau3 

* priors on p1, p2 are flat between 0.1 and 0.9 ==> 'pm.Uniform' variables 
with the constraint p2<=1-p1. p3 is deterministic ==1-p1-p2 

* the 'assignment' variable assigns each point to a gaussian, according to probabilities p1, p2, p3 

* priors on mu1, mu2, mu3 are flat between 40 and 120 ==> 'pm.Uniform' variables 

* priors on sigma1, sigma2, sigma3 are flat between 4 and 12 ==> 'pm.Uniform' variables 
""" 

    import numpy as np 
    import pymc as pm 

    data = np.loadtxt('distrib.txt') 
    Npts = len(data) 

    mumin = 40 
    mumax = 120 
    sigmamin=4 
    sigmamax=12 

    p1 = pm.Uniform("p1",0.1,0.9) 
    p2 = pm.Uniform("p2",0.1,1-p1) 
    p3 = 1-p1-p2 
    assignment = pm.Categorical('assignment',[p1,p2,p3],size=Npts) 
    mu = pm.Uniform('mu',[mumin,mumin,mumin],[mumax,mumax,mumax]) 
    sigma = pm.Uniform('sigma',[sigmamin,sigmamin,sigmamin], 
         [sigmamax,sigmamax,sigmamax]) 
    tau = 1/sigma**2 

    @pm.deterministic 
    def assign_mu(assi=assignment,mu=mu): 
     return mu[assi] 

    @pm.deterministic 
    def assign_tau(assi=assignment,sig=tau): 
     return sig[assi] 

    hypothesis = pm.Normal("obs", assign_mu, assign_tau, value=data, observed=True) 

    model = pm.Model([hypothesis, p1, p2, tau, mu]) 
    test = pm.MCMC(model) 
    test.sample(50000,burn=20000) # conservative values, let's take a coffee... 

    print('\nguess\n* p1, p2 = ', 
      np.mean(test.trace('p1')[:]),' ; ', 
      np.mean(test.trace('p2')[:]),' ==> p3 = ', 
      1-np.mean(test.trace('p1')[:])-np.mean(test.trace('p2')[:]), 
      '\n* mu = ', 
      np.mean(test.trace('mu')[:,0]),' ; ', 
      np.mean(test.trace('mu')[:,1]),' ; ', 
      np.mean(test.trace('mu')[:,2])) 

    print('why does this guess suck ???!!!')  

私はデータファイル 'distrib.txt'を送信できます。それは約500kbであり、データは以下にプロットされている。

p1, p2 = 0.366913192214 ; 0.583816452532 ==> p3 = 0.04927035525400003 
mu = 77.541619286 ; 75.3371615466 ; 77.2427165073 

が明らかにされている周りのバンプながら〜の周りの確率で55〜75及び〜90、〜0.2、〜0.5と〜0.3

Answer 1

ます：たとえば最後の実行は、私を与えましたNegative Binomial Mixture in PyMC

問題は、3つのコンポーネントの分布が近すぎるためにカテゴリ変数が収束するには遅すぎるという問題があります。

まず、我々はあなたのテストデータを生成します。

data1 = np.random.normal(55,5,2000) 
data2 = np.random.normal(75,5,5000) 
data3 = np.random.normal(90,5,3000) 
data=np.concatenate([data1, data2, data3]) 
np.savetxt("distrib.txt", data) 

その後、我々は事後グループの割り当てによって着色、ヒストグラムをプロットします

tablebyassignment = [data[np.nonzero(np.round(test.trace("assignment")[:].mean(axis=0)) == i)] for i in range(0,3) ] 
plt.hist(tablebyassingment, bins=30, stacked = True) 

これは、最終的には十分にすぐに収束しませんが、あなたに役立つこと。

あなたはMCMCを開始する前に、割り当ての値を推測して、この問題を解決することができます：あなたに与え

from sklearn.cluster import KMeans 
kme = KMeans(3) 
kme.fit(np.atleast_2d(data).T) 
assignment = pm.Categorical('assignment',[p1,p2,p3],size=Npts, value=kme.labels_) 

： 時間のすべてをすることができないカテゴリを初期化するために、K-手段を用いるが、収束しないよりも優れています。