How to code a hierarchical mixture model of multivariate normals using PYMC

Question

I successfully implemented a mixture of 3 normals using PyMC (shown at https://drive.google.com/file/d/0Bwnmbh6ueWhqSkUtV1JFZDJwLWc, and similar to the question asked at How to model a mixture of 3 Normals in PyMC?)

My next step is to try and code mixtures of multivariate normals.

There is, however, an additional complexity to the data - a hierarchy, with sets of observations belonging to a parent observation. The clustering is done on the parent observations, and not on the individual observations themselves. This first step generates the code (60 parents, with 50 observations per each parent), and works fine.

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

However, when I try and use PyMC to do the sampling, I run intro trouble ('error: failed in converting 3rd argument `tau' of flib.prec_mvnorm to C/Fortran array')

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)

I know the problem is not with the format of the simulated observed data, because this step would work fine, in place of the above:

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

As a result, I think the issue is how the mean vector ('mean') and the covariance matrix ('prec') are entered, I just don't know how. Like I said, this worked fine with mixtures of normal distributions, but mixtures of multivariate normals is adding a complexity I can't figure out.

Answer 1

This is a good example of the difficulty PyMC has with vectors of multivariate variables. Not that its difficult--just not as elegant as it should be. You should create a list comprehension of the MVN nodes and wrap that as an observed stochastic.

@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