title: How to code a hierarchical mixture model of multivariate normals using PYMC
body: <p>I successfully implemented a mixture of 3 normals using PyMC (shown at <a href="https://drive.google.com/file/d/0Bwnmbh6ueWhqSkUtV1JFZDJwLWc" rel="nofollow noreferrer">https://drive.google.com/file/d/0Bwnmbh6ueWhqSkUtV1JFZDJwLWc</a>, and similar to the question asked at <a href="https://stackoverflow.com/questions/18987697/how-to-model-a-mixture-of-3-normals-in-pymc">How to model a mixture of 3 Normals in PyMC?</a>)</p>
  
  <p>My next step is to try and code mixtures of multivariate normals. </p>
  
  <p>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.</p>
  
  <pre><code>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
  </code></pre>
  
  <p>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>
  
  <pre><code>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)
  </code></pre>
  
  <p>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:</p>
  
  <pre><code>obs = mc.MvNormal( "observed returns", post_mu3, post_cov_matrix_inv3, observed = True, value = data )
  </code></pre>
  
  <p>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. </p>