Problems modeling multivariate Gaussian distribution

# Problems modeling multivariate Gaussian distribution

• Friday, 27 April 2012 1:16 PM

Hello,

I'm trying to extend Bayesian PCA example to Bayesian mixture of FA, with full covariance matrix on the rows of matrix of factor loadings (vW in this case). But when I'm trying to define vW this way:

```vW[rC][rM] = Variable.Random(new VectorGaussian(Vector.Zero(vD),
PositiveDefiniteMatrix.IdentityScaledBy(vD, vAlpha[rC][rM])));```

compiler throws a bunch of errors like this: "cannot convert from 'MicrosoftResearch.Infer.Models.Variable<int>' to 'int' ". I understand what it means, but I can't figure out how to solve this. Could you please hint on how to do that?

Here is the full code:

```public class BayesianMFA
{
// Inference engine
public InferenceEngine engine = null;
// Model variables
public Variable<int> vN = null; // Number of data points
public Variable<int> vD = null; // Dimensionality of data
public Variable<int> vM = null; // Dimensionality of latent variables
public Variable<int> vC = null; // Number of components
public VariableArray2D<double> vData = null; // Data

public VariableArray<VariableArray<Vector>, Vector[][]> vW = null;
public VariableArray<VariableArray2D<double>, double[][,]> vZ = null;
public VariableArray<VariableArray2D<double>, double[][,]> vT = null;
public VariableArray<VariableArray2D<double>, double[][,]> vU = null;

public VariableArray<VariableArray<double>, double[][]> vMu = null;
public VariableArray<VariableArray<double>, double[][]> vPi = null;
public VariableArray<VariableArray<double>, double[][]> vAlpha = null;

// Priors - these are declared as distribution variables
// so that we can set them at run-time.
public VariableArray<Gamma> priorAlpha = null;
public VariableArray<Gaussian> priorMu = null;
public VariableArray<Gamma> priorPi = null;
// Index of mixture component for each data point
// and weights of the mixture
public VariableArray<int> index = null;
public Variable<Vector> weights = null;
// Ranges of vN, vD, vM, vC
public Range rN = null;
public Range rD = null;
public Range rM = null;
public Range rC = null;
/// <summary>
/// Model constructor
/// </summary>
public BayesianMFA()
{
// The various dimensions will be set externally...
vN = Variable.New<int>().Named("NumObs");
vD = Variable.New<int>().Named("NumFeats");
vM = Variable.New<int>().Named("MaxComponents");
vC = Variable.New<int>().Named("NumOfMixComponents");

rN = new Range(vN).Named("N");
rD = new Range(vD).Named("D");
rM = new Range(vM).Named("M");
rC = new Range(vC).Named("C");

// ... as will the data
vData = Variable.Array<double>(rN, rD).Named("data");

// ... and the priors
priorAlpha = Variable.Array<Gamma>(rC).Named("PriorAlpha");
priorMu = Variable.Array<Gaussian>(rC).Named("PriorMu");
priorPi = Variable.Array<Gamma>(rC).Named("PriorPi");

// Mixing matrix
vAlpha = Variable.Array(Variable.Array<double>(rM), rC).Named("Alpha");
vW = Variable.Array(Variable.Array<Vector>(rM), rC).Named("W");
vAlpha[rC][rM] = Variable.Random<double, Gamma>(priorAlpha[rC]).ForEach(rM);
vW[rC][rM] = Variable.Random(new VectorGaussian(Vector.Zero(vD),
PositiveDefiniteMatrix.IdentityScaledBy(vD, vAlpha[rC][rM])));

// Latent variables are drawn from a standard Gaussian
vZ = Variable.Array(Variable.Array<double>(rN,rM), rC).Named("Z");
vZ[rC][rN,rM] = Variable.GaussianFromMeanAndPrecision(0.0, 1.0).ForEach(rN, rM, rC);

// Multiply the latent variables with the mixing matrix...
vT = Variable.Array(Variable.Array<double>(rN, rD), rC).Named("T");
VariableArray<VariableArray<VariableArray<double>, double[][]>, double[][][]> vWMatrix;
vWMatrix = Variable.Array(Variable.Array(Variable.Array<double>(rD),rM),rC).Named("U");
vWMatrix[rC][rM] = Variable.ArrayFromVector(vW[rC][rM], rD);
// Explicit Matrix multiply
using (Variable.ForEach(rC))
{
using (Variable.ForEach(rN))
{
using (Variable.ForEach(rD))
{
var vProds = Variable.Array<double>(rM);
vProds[rM] = vZ[rC][rN,rM] * vWMatrix[rC][rM][rD];
vT[rC][rN,rD] = Variable.Sum(vProds);
}
}
}

// ... add in a bias ...
vMu = Variable.Array(Variable.Array<double>(rD), rC).Named("mu");
vMu[rC][rD] = Variable.Random<double, Gaussian>(priorMu[rC]).ForEach(rD);
vU = Variable.Array(Variable.Array<double>(rN, rD), rC).Named("U");
vU[rC][rN, rD] = vT[rC][rN, rD] + vMu[rC][rD];

// ... and add in some observation noise ...
vPi = Variable.Array(Variable.Array<double>(rD), rC).Named("pi");
vPi[rC][rD] = Variable.Random<double, Gamma>(priorPi[rC]).ForEach(rD);

// Working with the mixture
weights = Variable.New<Vector>().Named("weights");
weights = Variable.DirichletSymmetric(rC, 1);
index = Variable.Array<int>(rN).Named("index");
using (Variable.ForEach(rN))
{
index[rN] = Variable.Discrete(weights);
using (Variable.Switch(index[rN]))
{
vData[rN, rD] = Variable.GaussianFromMeanAndPrecision(vU[index[rN]][rN, rD], vPi[index[rN]][rD]);
}
}
// Inference engine
engine = new InferenceEngine(new VariationalMessagePassing());
return;
}
}```

Thanks a lot.

### All Replies

• Friday, 27 April 2012 6:50 PM

The issue is how to create a VectorGaussian variable with a diagonal precision matrix.  This is possible by first creating a uniform VectorGaussian and then applying constraints, as illustrated in the code below:
```public void VectorGaussianFromPrecisionDiagonalExample()
{
Range d = new Range(3);
VariableArray<double> prec = Variable.Array<double>(d);
prec[d] = Variable.GammaFromShapeAndRate(1,1).ForEach(d);
Variable<Vector> x = VectorGaussianFromPrecisionDiagonal(prec);
InferenceEngine engine = new InferenceEngine(new VariationalMessagePassing());
Console.WriteLine(engine.Infer(x));
}
public Variable<Vector> VectorGaussianFromPrecisionDiagonal(VariableArray<double> precisionDiagonal)
{
Range r = precisionDiagonal.Range;
int dimension = r.SizeAsInt;
Variable<Vector> x = Variable.VectorGaussianFromMeanAndPrecision(Vector.Zero(dimension), PositiveDefiniteMatrix.IdentityScaledBy(dimension, 1e-20));
using (var fb = Variable.ForEach(r)) {
Variable<double> xi = Variable.GetItem(x, fb.Index);
var y = Variable.GaussianFromMeanAndPrecision(xi, precisionDiagonal[r]);
y.ObservedValue = 0;
}
return x;
}
```
This technique allows essentially any structure to be imposed on the precision matrix.