PyMC を使用して 3 つの法線の混合を正常に実装しました ( https://drive.google.com/file/d/0Bwnmbh6ueWhqSkUtV1JFZDJwLWcに表示され、PyMC で 3 つの法線の混合物をモデル化する方法で尋ねられた質問と同様です) 。
私の次のステップは、多変量法線の混合を試してコーディングすることです。
ただし、データには追加の複雑さがあります。つまり、一連の観測が親観測に属する階層です。クラスタリングは、個々の観測自体ではなく、親の観測に対して行われます。この最初のステップでコード (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 を使用してサンプリングを実行しようとすると、イントロ トラブルが発生します (「エラー: flib.prec_mvnorm の 3 番目の引数 `tau' を 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') の入力方法だと思いますが、方法がわかりません。私が言ったように、これは正規分布の混合ではうまくいきましたが、多変量正規分布の混合は私が理解できない複雑さを追加しています.