RMOL Logo  0.25.2
C++ library of Revenue Management and Optimisation classes and functions
mog.cpp
Go to the documentation of this file.
00001 #include <itpp/itstat.h>
00002 
00003 #include <fstream>
00004 #include <iostream>
00005 #include <iomanip>
00006 #include <ios>
00007 
00008 using std::cout;
00009 using std::endl;
00010 using std::fixed;
00011 using std::setprecision;
00012 
00013 using namespace itpp;
00014 
00015 int main() {
00016 
00017   bool print_progress = false;
00018 
00019   //
00020   // first, let's generate some synthetic data
00021 
00022   int N = 100000;  // number of vectors
00023   int D = 3;       // number of dimensions
00024   int K = 5;       // number of Gaussians
00025 
00026   Array<vec> X(N); for(int n=0;n<N;n++) { X(n).set_size(D); X(n) = 0.0; }
00027 
00028   // the means
00029 
00030   Array<vec> mu(K);
00031   mu(0) = "-6, -4, -2";
00032   mu(1) = "-4, -2,  0";
00033   mu(2) = "-2,  0,  2";
00034   mu(3) = " 0, +2, +4";
00035   mu(4) = "+2, +4, +6";
00036 
00037 
00038   // the diagonal variances
00039 
00040   Array<vec> var(K);
00041   var(0) = "0.1, 0.2, 0.3";
00042   var(1) = "0.2, 0.3, 0.1";
00043   var(2) = "0.3, 0.1, 0.2";
00044   var(3) = "0.1, 0.2, 0.3";
00045   var(4) = "0.2, 0.3, 0.1";
00046   
00047   cout << fixed << setprecision(3);
00048   cout << "user configured means and variances:" << endl;
00049   cout << "mu = " << mu << endl;
00050   cout << "var = " << var << endl;
00051   
00052   // randomise the order of Gaussians "generating" the vectors
00053   I_Uniform_RNG rnd_uniform(0, K-1);
00054   ivec gaus_id = rnd_uniform(N); 
00055 
00056   ivec gaus_count(K); gaus_count = 0;
00057   Array<vec> mu_test(K);  for(int k=0;k<K;k++) { mu_test(k).set_size(D); mu_test(k) = 0.0; }
00058   Array<vec> var_test(K);  for(int k=0;k<K;k++) { var_test(k).set_size(D); var_test(k) = 0.0; }
00059 
00060   Normal_RNG rnd_normal;
00061   for(int n=0;n<N;n++) {
00062     
00063     int k = gaus_id(n);
00064     gaus_count(k)++;
00065    
00066     for(int d=0;d<D;d++) {
00067       rnd_normal.setup( mu(k)(d), var(k)(d) );
00068       double tmp = rnd_normal();
00069       X(n)(d) = tmp; 
00070       mu_test(k)(d) += tmp;
00071     }
00072   }
00073 
00074   //
00075   // find the stats for the generated data
00076 
00077   for(int k=0;k<K;k++) mu_test(k) /= gaus_count(k);
00078 
00079   for(int n=0;n<N;n++) {
00080     int k = gaus_id(n);
00081     
00082     for(int d=0;d<D;d++) {
00083       double tmp = X(n)(d) - mu_test(k)(d);
00084       var_test(k)(d) += tmp*tmp;  
00085     }
00086   }
00087 
00088   for(int k=0;k<K;k++)  var_test(k) /= (gaus_count(k)-1.0); 
00089 
00090   cout << endl << endl;
00091   cout << fixed << setprecision(3);
00092   cout << "stats for X:" << endl;
00093   
00094   for(int k=0;k<K;k++) {
00095     cout << "k = " << k << "  count = " << gaus_count(k) << "  weight = " << gaus_count(k)/double(N) << endl; 
00096     for(int d=0;d<D;d++) cout << "  d = " << d << "  mu_test = " << mu_test(k)(d) << "  var_test = " << var_test(k)(d) << endl; 
00097     cout << endl;
00098   }
00099  
00100   
00101   // make a model with initial values (zero mean and unit variance)
00102   // the number of gaussians and dimensions of the model is specified here
00103 
00104   MOG_diag mog(K,D);
00105 
00106   cout << endl;
00107   cout << fixed << setprecision(3);
00108   cout << "mog.avg_log_lhood(X) = " << mog.avg_log_lhood(X) << endl;
00109 
00110   //
00111   // find initial parameters via k-means (which are then used as seeds for EM based optimisation)
00112 
00113   cout << endl << endl;
00114   cout << "running kmeans optimiser" << endl << endl;
00115   
00116   MOG_diag_kmeans(mog, X, 10, 0.5, true, print_progress);
00117   
00118   cout << fixed << setprecision(3);
00119   cout << "mog.get_means() = " << endl << mog.get_means() << endl; 
00120   cout << "mog.get_diag_covs() = " << endl << mog.get_diag_covs() << endl;
00121   cout << "mog.get_weights() = " << endl << mog.get_weights() << endl;
00122 
00123   cout << endl;
00124   cout << "mog.avg_log_lhood(X) = " << mog.avg_log_lhood(X) << endl;
00125 
00126   
00127   //
00128   // EM ML based optimisation
00129 
00130   cout << endl << endl;
00131   cout << "running ML optimiser" << endl << endl;
00132   
00133   MOG_diag_ML(mog, X, 10, 0.0, 0.0, print_progress);
00134 
00135   cout << fixed << setprecision(3);
00136   cout << "mog.get_means() = " << endl << mog.get_means() << endl;
00137   cout << "mog.get_diag_covs() = " << endl << mog.get_diag_covs() << endl;
00138   cout << "mog.get_weights() = " << endl << mog.get_weights() << endl;
00139 
00140   cout << endl;
00141   cout << "mog.avg_log_lhood(X) = " << mog.avg_log_lhood(X) << endl;
00142 
00143   return 0;
00144 }