如何代碼中使用PYMC

Question

我成功實施使用PyMC 3個法線混合多元法線的分層混合模型（在https://drive.google.com/file/d/0Bwnmbh6ueWhqSkUtV1JFZDJwLWc所示，並且類似於在How to model a mixture of 3 Normals in PyMC?問的問題）

我的下一步是嘗試代碼混合多元法線。

然而，數據有一個額外的複雜性 - 一個層次結構，觀察組屬於父親觀察。聚類是在父親觀察中完成的，而不是在個人觀察本身上進行的。該第一步驟生成的代碼（60個父母，每每個父50個觀察值），並且工作正常。

import numpy as np 
import pymc as mc 
n = 3 #mixtures 
B = 5 #Bias between those at different mixtures 
tau = 3 #Variances 
nprov = 60 #number of parent observations 
mu = [[0,0],[0,B],[-B,0]] 
true_cov0 = np.array([[1.,0.],[0.,1.]]) 
true_cov1 = np.array([[1.,0.],[0.,tau**(2)]]) 
true_cov2 = np.array([[tau**(-2),0],[0.,1.]]) 
trueprobs = [.4, .3, .3] #probability of being in each of the three mixtures 

prov = np.random.multinomial(1, trueprobs, size=nprov) 
v = prov[:,1] + (prov[:,2])*2 
numtoeach = 50 
n_obs = nprov*numtoeach 
vAll = np.tile(v,numtoeach) 
ndata = numtoeach*nprov 
p1 = range(nprov) 
prov1 = np.tile(p1,numtoeach) 

data = (vAll==0)*(np.random.multivariate_normal(mu[0],true_cov0,ndata)).T \ 
    + (vAll==1)*(np.random.multivariate_normal(mu[1],true_cov1,ndata)).T \ 
    + (vAll==2)*(np.random.multivariate_normal(mu[2],true_cov2,ndata)).T 
data=data.T 

然而，當我嘗試使用PyMC做取樣，我運行前奏麻煩（「錯誤：失敗第三參數`tau蛋白轉化」 flib.prec_mvnorm到C/Fortran數組'）

p = 2 #covariates 
prior_mu1=np.ones(p) 
prior_mu2=np.ones(p) 
prior_mu3=np.ones(p) 
post_mu1 = mc.Normal("returns1",prior_mu1,1,size=p) 
post_mu2 = mc.Normal("returns2",prior_mu2,1,size=p) 
post_mu3 = mc.Normal("returns3",prior_mu3,1,size=p) 
post_cov_matrix_inv1 = mc.Wishart("cov_matrix_inv1",n_obs,np.eye(p)) 
post_cov_matrix_inv2 = mc.Wishart("cov_matrix_inv2",n_obs,np.eye(p)) 
post_cov_matrix_inv3 = mc.Wishart("cov_matrix_inv3",n_obs,np.eye(p)) 

#Combine prior means and variance matrices 
meansAll= np.array([post_mu1,post_mu2,post_mu3]) 
precsAll= np.array([post_cov_matrix_inv1,post_cov_matrix_inv2,post_cov_matrix_inv3]) 

dd = mc.Dirichlet('dd', theta=(1,)*n) 
category = mc.Categorical('category', p=dd, size=nprov) 

#This step accounts for the hierarchy: observations' means are equal to their parents mean 
#Parent is labeled prov1 

@mc.deterministic 
def mean(category=category, meansAll=meansAll): 
lat = category[prov1] 
new = meansAll[lat] 
return new 

@mc.deterministic 
def prec(category=category, precsAll=precsAll): 
lat = category[prov1] 
return precsAll[lat] 

obs = mc.MvNormal("observed returns", mean, prec, observed = True, value = data) 

我知道這個問題是不是與模擬觀測數據的格式，因爲這一步就做工精細，代替上述的：

obs = mc.MvNormal("observed returns", post_mu3, post_cov_matrix_inv3, observed = True, value = data) 

因此，我認爲這個問題是怎麼意思是向量（'mean'）和協方差矩陣（'prec'）被輸入，我只是不知道如何。就像我說的，這與正態分佈的混合正常工作，但多元正態分佈的混合正在增加我無法弄清楚的複雜性。

Answer 1

這是PyMC與多元變量向量難度很好的一個例子。並不是說它很難 - 只是沒有它應該的那麼優雅。您應該創建MVN節點列表理解和包裝，作爲一個觀察隨機的。

@mc.observed 
def obs(value=data, mean=mean, prec=prec): 
    return sum(mc.mv_normal_like(v, m, T) for v,m,T in zip(data, mean, prec)) 

Here is the IPython notebook