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 }
Generated on Tue Sep 7 2010 18:30:39 for RMOL by Doxygen 1.7.1