Using Mixture of Gaussians (MOG) module to model data

This example demonstrates how to find the parameters of a MOG model via using the kmeans and EM based optimisers. Synthetic data is utilised.

#include <itpp/itstat.h>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <ios>
using std::cout;
using std::endl;
using std::fixed;
using std::setprecision;
using namespace itpp;
int main()
{
bool print_progress = false;
//
// first, let's generate some synthetic data
int N = 100000; // number of vectors
int D = 3; // number of dimensions
int K = 5; // number of Gaussians
Array<vec> X(N);
for (int n = 0;n < N;n++) { X(n).set_size(D); X(n) = 0.0; }
// the means
Array<vec> mu(K);
mu(0) = "-6, -4, -2";
mu(1) = "-4, -2, 0";
mu(2) = "-2, 0, 2";
mu(3) = " 0, +2, +4";
mu(4) = "+2, +4, +6";
// the diagonal variances
Array<vec> var(K);
var(0) = "0.1, 0.2, 0.3";
var(1) = "0.2, 0.3, 0.1";
var(2) = "0.3, 0.1, 0.2";
var(3) = "0.1, 0.2, 0.3";
var(4) = "0.2, 0.3, 0.1";
cout << fixed << setprecision(3);
cout << "user configured means and variances:" << endl;
cout << "mu = " << mu << endl;
cout << "var = " << var << endl;
// randomise the order of Gaussians "generating" the vectors
I_Uniform_RNG rnd_uniform(0, K - 1);
ivec gaus_id = rnd_uniform(N);
ivec gaus_count(K);
gaus_count = 0;
Array<vec> mu_test(K);
for (int k = 0;k < K;k++) { mu_test(k).set_size(D); mu_test(k) = 0.0; }
Array<vec> var_test(K);
for (int k = 0;k < K;k++) { var_test(k).set_size(D); var_test(k) = 0.0; }
Normal_RNG rnd_normal;
for (int n = 0;n < N;n++) {
int k = gaus_id(n);
gaus_count(k)++;
for (int d = 0;d < D;d++) {
rnd_normal.setup(mu(k)(d), var(k)(d));
double tmp = rnd_normal();
X(n)(d) = tmp;
mu_test(k)(d) += tmp;
}
}
//
// find the stats for the generated data
for (int k = 0;k < K;k++) mu_test(k) /= gaus_count(k);
for (int n = 0;n < N;n++) {
int k = gaus_id(n);
for (int d = 0;d < D;d++) {
double tmp = X(n)(d) - mu_test(k)(d);
var_test(k)(d) += tmp * tmp;
}
}
for (int k = 0;k < K;k++) var_test(k) /= (gaus_count(k) - 1.0);
cout << endl << endl;
cout << fixed << setprecision(3);
cout << "stats for X:" << endl;
for (int k = 0;k < K;k++) {
cout << "k = " << k << " count = " << gaus_count(k) << " weight = " << gaus_count(k) / double(N) << endl;
for (int d = 0;d < D;d++) cout << " d = " << d << " mu_test = " << mu_test(k)(d) << " var_test = " << var_test(k)(d) << endl;
cout << endl;
}
// make a model with initial values (zero mean and unit variance)
// the number of gaussians and dimensions of the model is specified here
MOG_diag mog(K, D);
cout << endl;
cout << fixed << setprecision(3);
cout << "mog.avg_log_lhood(X) = " << mog.avg_log_lhood(X) << endl;
//
// find initial parameters via k-means (which are then used as seeds for EM based optimisation)
cout << endl << endl;
cout << "running kmeans optimiser" << endl << endl;
MOG_diag_kmeans(mog, X, 10, 0.5, true, print_progress);
cout << fixed << setprecision(3);
cout << "mog.get_means() = " << endl << mog.get_means() << endl;
cout << "mog.get_diag_covs() = " << endl << mog.get_diag_covs() << endl;
cout << "mog.get_weights() = " << endl << mog.get_weights() << endl;
cout << endl;
cout << "mog.avg_log_lhood(X) = " << mog.avg_log_lhood(X) << endl;
//
// EM ML based optimisation
cout << endl << endl;
cout << "running ML optimiser" << endl << endl;
MOG_diag_ML(mog, X, 10, 0.0, 0.0, print_progress);
cout << fixed << setprecision(3);
cout << "mog.get_means() = " << endl << mog.get_means() << endl;
cout << "mog.get_diag_covs() = " << endl << mog.get_diag_covs() << endl;
cout << "mog.get_weights() = " << endl << mog.get_weights() << endl;
cout << endl;
cout << "mog.avg_log_lhood(X) = " << mog.avg_log_lhood(X) << endl;
return 0;
}

Generated on Thu Mar 20 2014 12:26:44 for IT++ by Doxygen 1.8.6