this code is to be optimized and made faster
This commit is contained in:
Parikshit Ram
@@ -0,0 +1,26 @@
|
||||
|
||||
librule(
|
||||
name = "mog_em", # this line can be safely omitted
|
||||
sources = ["mog.cc"], # files that must be compiled
|
||||
headers = ["mog.h","phi.h","math_functions.h"],# include files part of the 'lib'
|
||||
deplibs = ["fastlib:fastlib"], # depends on fastlib core
|
||||
#tests = ["mog_em_tests.cc"] # this file contains a main with test functions
|
||||
)
|
||||
|
||||
binrule(
|
||||
name = "mog_em_main", # the executable name
|
||||
sources = ["mog_em_main.cc"], # compile main.cc
|
||||
#headers = ["mog_em_main.h"], # no extra headers
|
||||
deplibs = [":mog_em", "fastlib:fastlib"] #
|
||||
)
|
||||
|
||||
# to build:
|
||||
# 1. make sure have environment variables set up:
|
||||
# $ source /full/path/to/fastlib/script/fl-env /full/path/to/fastlib
|
||||
# (you might want to put this in bashrc)
|
||||
# 2. fl-build main
|
||||
# - this automatically will assume --mode=check, the default
|
||||
# - type fl-build --help for help
|
||||
# 3. ./main
|
||||
# - to build same target again, type: make
|
||||
# - to force recompilation, type: make clean
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,138 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file math_functions.h
|
||||
*
|
||||
* This file has certain functions that find the
|
||||
* highest or lowest element in an array or
|
||||
* in a row of a matrix and returns them
|
||||
*
|
||||
*/
|
||||
|
||||
#include "fastlib/fastlib.h"
|
||||
#include "fastlib/fastlib_int.h"
|
||||
|
||||
|
||||
/**
|
||||
* Finds the index of the minimum element
|
||||
* in each row of a matrix
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* Matrix& mat;
|
||||
* index_t indices[mat.n_rows()];
|
||||
*
|
||||
* ...
|
||||
*
|
||||
* min_element(mat, indices);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
void min_element(Matrix& element, index_t *indices) {
|
||||
|
||||
index_t last = element.n_cols() - 1;
|
||||
index_t first, lowest;
|
||||
index_t i;
|
||||
|
||||
for (i = 0; i < element.n_rows(); i++) {
|
||||
|
||||
first = lowest = 0;
|
||||
if (first == last) {
|
||||
indices[i] = last;
|
||||
}
|
||||
while (++first <= last) {
|
||||
if (element.get(i, first) < element.get(i, lowest)) {
|
||||
lowest = first;
|
||||
}
|
||||
}
|
||||
indices[i] = lowest;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the index of the maximum element
|
||||
* in a float array 'array' of length 'length'.
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* index_t length, index;
|
||||
* float array[length];
|
||||
* ...
|
||||
* index = max_element_index(array, length);
|
||||
* @endcode
|
||||
*/
|
||||
int max_element_index(float *array, int length) {
|
||||
|
||||
int last = length - 1;
|
||||
int first = 0;
|
||||
int highest = 0;
|
||||
|
||||
if (first == last) {
|
||||
return last;
|
||||
}
|
||||
while (++first <= last) {
|
||||
if (array[first] > array[highest]) {
|
||||
highest = first;
|
||||
}
|
||||
}
|
||||
return highest;
|
||||
}
|
||||
|
||||
/**
|
||||
* Finds the index of the maximum element in an arraylist
|
||||
* of floats
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* index_t index;
|
||||
* ArrayList<float> array;
|
||||
* ...
|
||||
* index = max_element_index(array);
|
||||
* @endcode
|
||||
*/
|
||||
int max_element_index(ArrayList<float>& array){
|
||||
|
||||
int last = array.size() - 1;
|
||||
int first = 0;
|
||||
int highest = 0;
|
||||
|
||||
if (first == last) {
|
||||
return last;
|
||||
}
|
||||
while (++first <= last) {
|
||||
if (array[first] > array[highest]) {
|
||||
highest = first;
|
||||
}
|
||||
}
|
||||
return highest;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the index of the maximum element in
|
||||
* an arraylist of doubles
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* index_t index;
|
||||
* ArrayList<double> array;
|
||||
* ...
|
||||
* index = max_element_index(array);
|
||||
* @endcode
|
||||
*/
|
||||
int max_element_index(ArrayList<double>& array) {
|
||||
|
||||
int last = array.size() - 1;
|
||||
int first = 0;
|
||||
int highest = 0;
|
||||
|
||||
if (first == last) {
|
||||
return last;
|
||||
}
|
||||
while (++first <= last) {
|
||||
if(array[first] > array[highest]) {
|
||||
highest = first;
|
||||
}
|
||||
}
|
||||
return highest;
|
||||
}
|
||||
|
||||
@@ -0,0 +1,301 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog.cc
|
||||
*
|
||||
* Implementation for the loglikelihood function, the EM algorithm
|
||||
* and also computes the K-means for getting an initial point
|
||||
*
|
||||
*/
|
||||
|
||||
#include "mog.h"
|
||||
#include "phi.h"
|
||||
#include "math_functions.h"
|
||||
|
||||
void MoGEM::ExpectationMaximization(Matrix& data_points) {
|
||||
|
||||
// Declaration of the variables */
|
||||
index_t num_points;
|
||||
index_t dim, num_gauss;
|
||||
double sum, tmp;
|
||||
ArrayList<Vector> mu_temp, mu;
|
||||
ArrayList<Matrix> sigma_temp, sigma;
|
||||
Vector omega_temp, omega, x;
|
||||
Matrix cond_prob;
|
||||
long double l, l_old, best_l, INFTY = 99999, TINY = 1.0e-10;
|
||||
|
||||
// Initializing values
|
||||
dim = dimension();
|
||||
num_gauss = number_of_gaussians();
|
||||
num_points = data_points.n_cols();
|
||||
|
||||
// Initializing the number of the vectors and matrices
|
||||
// according to the parameters input
|
||||
mu_temp.Init(num_gauss);
|
||||
mu.Init(num_gauss);
|
||||
sigma_temp.Init(num_gauss);
|
||||
sigma.Init(num_gauss);
|
||||
omega_temp.Init(num_gauss);
|
||||
omega.Init(num_gauss);
|
||||
|
||||
// Allocating size to the vectors and matrices
|
||||
// according to the dimensionality of the data
|
||||
for(index_t i = 0; i < num_gauss; i++) {
|
||||
mu_temp[i].Init(dim);
|
||||
mu[i].Init(dim);
|
||||
sigma_temp[i].Init(dim, dim);
|
||||
sigma[i].Init(dim, dim);
|
||||
}
|
||||
x.Init(dim);
|
||||
cond_prob.Init(num_gauss, num_points);
|
||||
|
||||
best_l = -INFTY;
|
||||
index_t restarts = 0;
|
||||
// performing 5 restarts and choosing the best from them
|
||||
while (restarts < 1) {
|
||||
|
||||
// assign initial values to 'mu', 'sig' and 'omega' using k-means
|
||||
KMeans(data_points, &mu_temp, &sigma_temp, &omega_temp, num_gauss);
|
||||
|
||||
l_old = -INFTY;
|
||||
|
||||
// calculates the loglikelihood value
|
||||
l = Loglikelihood(data_points, mu_temp, sigma_temp, omega_temp);
|
||||
|
||||
// added a check here to see if any
|
||||
// significant change is being made
|
||||
// at every iteration
|
||||
while (l - l_old > TINY) {
|
||||
// calculating the conditional probabilities
|
||||
// of choosing a particular gaussian given
|
||||
// the data and the present theta value
|
||||
/**********
|
||||
for (index_t j = 0; j < num_points; j++) {
|
||||
x.CopyValues(data_points.GetColumnPtr(j));
|
||||
sum = 0;
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
tmp = phi(x, mu_temp[i], sigma_temp[i]) * omega_temp.get(i); // can be made faster by sending all the data at once instead of looping over
|
||||
cond_prob.set(i, j, tmp);
|
||||
sum += tmp;
|
||||
}
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
tmp = cond_prob.get(i, j);
|
||||
cond_prob.set(i, j, tmp / sum);
|
||||
}
|
||||
}
|
||||
******/
|
||||
|
||||
/*******FIX THIS SOON********/
|
||||
Vector phi_val;
|
||||
phi_val.Init(num_points);
|
||||
|
||||
for(index_t i = 0; i < num_gauss; i++) {
|
||||
phi_val.CopyValues(phi(data_points, mu_temp[i], sigma_temp[i]));
|
||||
la::Scale(omega_temp.get(i), &phi_val);
|
||||
*/
|
||||
|
||||
// calculating the new value of the mu
|
||||
// using the updated conditional probabilities
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
sum = 0;
|
||||
mu_temp[i].SetZero();
|
||||
for (index_t j = 0; j < num_points; j++) {
|
||||
x.CopyValues(data_points.GetColumnPtr(j));
|
||||
la::AddExpert(cond_prob.get(i, j), x, &mu_temp[i]);
|
||||
sum += cond_prob.get(i, j);
|
||||
}
|
||||
la::Scale((1.0 / sum), &mu_temp[i]);
|
||||
}
|
||||
|
||||
// calculating the new value of the sig
|
||||
// using the updated conditional probabilities
|
||||
// and the updated mu
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
sum = 0;
|
||||
sigma_temp[i].SetZero();
|
||||
for (index_t j = 0; j < num_points; j++) {
|
||||
Matrix co, ro, c;
|
||||
c.Init(dim, dim);
|
||||
x.CopyValues(data_points.GetColumnPtr(j));
|
||||
la::SubFrom(mu_temp[i] , &x);
|
||||
co.AliasColVector(x);
|
||||
ro.AliasRowVector(x);
|
||||
la::MulOverwrite(co, ro, &c);
|
||||
la::AddExpert(cond_prob.get(i, j), c, &sigma_temp[i]);
|
||||
sum += cond_prob.get(i, j);
|
||||
}
|
||||
la::Scale((1.0 / sum), &sigma_temp[i]);
|
||||
}
|
||||
|
||||
// calculating the new values for omega
|
||||
// using the updated conditional probabilities
|
||||
Vector identity_vector;
|
||||
identity_vector.Init(num_points);
|
||||
identity_vector.SetAll(1.0 / num_points);
|
||||
la::MulOverwrite(cond_prob, identity_vector, &omega_temp);
|
||||
|
||||
l_old = l;
|
||||
l = Loglikelihood(data_points, mu_temp, sigma_temp, omega_temp);
|
||||
}
|
||||
|
||||
// putting a check to see if the best one is chosen
|
||||
if(l > best_l){
|
||||
best_l = l;
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
mu[i].CopyValues(mu_temp[i]);
|
||||
sigma[i].CopyValues(sigma_temp[i]);
|
||||
}
|
||||
omega.CopyValues(omega_temp);
|
||||
}
|
||||
restarts++;
|
||||
}
|
||||
|
||||
for (index_t i = 0; i < num_gauss; i++) {
|
||||
set_mu(i, mu[i]);
|
||||
set_sigma(i, sigma[i]);
|
||||
}
|
||||
set_omega(omega);
|
||||
|
||||
NOTIFY("loglikelihood value of the estimated model: %Lf\n", best_l);
|
||||
return;
|
||||
}
|
||||
|
||||
long double MoGEM::Loglikelihood(Matrix& data_points, ArrayList<Vector>& means,
|
||||
ArrayList<Matrix>& covars, Vector& weights) {
|
||||
|
||||
index_t i, j;
|
||||
Vector x;
|
||||
long double likelihood, loglikelihood = 0;
|
||||
|
||||
x.Init(data_points.n_rows());
|
||||
|
||||
for (j = 0; j < data_points.n_cols(); j++) {
|
||||
x.CopyValues(data_points.GetColumnPtr(j));
|
||||
likelihood = 0;
|
||||
for(i = 0; i < number_of_gaussians() ; i++){
|
||||
likelihood += weights.get(i) * phi(x, means[i], covars[i]);
|
||||
}
|
||||
loglikelihood += log(likelihood);
|
||||
}
|
||||
return loglikelihood;
|
||||
}
|
||||
|
||||
void MoGEM::KMeans(Matrix& data, ArrayList<Vector> *means,
|
||||
ArrayList<Matrix> *covars, Vector *weights, index_t value_of_k){
|
||||
|
||||
|
||||
ArrayList<Vector> mu, mu_old;
|
||||
double* tmpssq;
|
||||
double* sig;
|
||||
double* sig_best;
|
||||
index_t *y;
|
||||
Vector x, diff;
|
||||
Matrix ssq;
|
||||
index_t i, j, k, n, t, dim;
|
||||
double score, score_old, sum;
|
||||
|
||||
n = data.n_cols();
|
||||
dim = data.n_rows();
|
||||
mu.Init(value_of_k);
|
||||
mu_old.Init(value_of_k);
|
||||
tmpssq = (double*)malloc(value_of_k * sizeof( double ));
|
||||
sig = (double*)malloc(value_of_k * sizeof( double ));
|
||||
sig_best = (double*)malloc(value_of_k * sizeof( double ));
|
||||
ssq.Init(n, value_of_k);
|
||||
|
||||
for( i = 0; i < value_of_k; i++){
|
||||
mu[i].Init(dim);
|
||||
mu_old[i].Init(dim);
|
||||
}
|
||||
x.Init(dim);
|
||||
y = (index_t*)malloc(n * sizeof(index_t));
|
||||
diff.Init(dim);
|
||||
|
||||
score_old = 999999;
|
||||
|
||||
// putting 5 random restarts to obtain the k-means
|
||||
for(i = 0; i < 5; i++){
|
||||
t = -1;
|
||||
for (k = 0; k < value_of_k; k++){
|
||||
t = (t + 1 + (rand()%((n - 1 - (value_of_k - k)) - (t + 1))));
|
||||
mu[k].CopyValues(data.GetColumnPtr(t));
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues( data.GetColumnPtr(j));
|
||||
la::SubOverwrite(mu[k], x, &diff);
|
||||
ssq.set( j, k, la::Dot(diff, diff));
|
||||
}
|
||||
}
|
||||
min_element(ssq, y);
|
||||
|
||||
do{
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
mu_old[k].CopyValues(mu[k]);
|
||||
}
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
index_t p = 0;
|
||||
mu[k].SetZero();
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues(data.GetColumnPtr(j));
|
||||
if(y[j] == k){
|
||||
la::AddTo(x, &mu[k]);
|
||||
p++;
|
||||
}
|
||||
}
|
||||
|
||||
if(p == 0){
|
||||
}
|
||||
else{
|
||||
double sc = 1 ;
|
||||
sc = sc / p;
|
||||
la::Scale(sc , &mu[k]);
|
||||
}
|
||||
for(j = 0; j < n; j++){
|
||||
x.CopyValues(data.GetColumnPtr(j));
|
||||
la::SubOverwrite(mu[k], x, &diff);
|
||||
ssq.set(j, k, la::Dot(diff, diff));
|
||||
}
|
||||
}
|
||||
min_element(ssq, y);
|
||||
|
||||
sum = 0;
|
||||
for(k = 0; k < value_of_k; k++) {
|
||||
la::SubOverwrite(mu[k], mu_old[k], &diff);
|
||||
sum += la::Dot(diff, diff);
|
||||
}
|
||||
}while(sum != 0);
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
index_t p = 0;
|
||||
tmpssq[k] = 0;
|
||||
for(j = 0; j < n; j++){
|
||||
if(y[j] == k){
|
||||
tmpssq[k] += ssq.get(j, k);
|
||||
p++;
|
||||
}
|
||||
}
|
||||
sig[k] = sqrt(tmpssq[k] / p);
|
||||
}
|
||||
|
||||
score = 0;
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
score += tmpssq[k];
|
||||
}
|
||||
score = score / n;
|
||||
|
||||
if (score < score_old) {
|
||||
score_old = score;
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
(*means)[k].CopyValues(mu[k]);
|
||||
sig_best[k] = sig[k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(k = 0; k < value_of_k; k++){
|
||||
x.SetAll(sig_best[k]);
|
||||
(*covars)[k].SetDiagonal(x);
|
||||
}
|
||||
double tmp = 1;
|
||||
(*weights).SetAll(tmp / value_of_k);
|
||||
return;
|
||||
}
|
||||
@@ -0,0 +1,244 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog.h
|
||||
*
|
||||
* Defines a Gaussian Mixture model and
|
||||
* estimates the parameters of the model
|
||||
*/
|
||||
|
||||
#ifndef MOGEM_H
|
||||
#define MOGEM_H
|
||||
|
||||
#include <fastlib/fastlib.h>
|
||||
|
||||
/**
|
||||
* A Gaussian mixture model class.
|
||||
*
|
||||
* This class uses maximum likelihood loss functions to
|
||||
* estimate the parameters of a gaussian mixture
|
||||
* model on a given data via the EM algorithm.
|
||||
*
|
||||
*
|
||||
* Example use:
|
||||
*
|
||||
* @code
|
||||
* MoGEM mog;
|
||||
* ArrayList<double> results;
|
||||
*
|
||||
* mog.Init(number_of_gaussians, dimension);
|
||||
* mog.ExpectationMaximization(data, &results, optim_flag);
|
||||
* @endcode
|
||||
*/
|
||||
class MoGEM {
|
||||
|
||||
private:
|
||||
|
||||
// The parameters of the mixture model
|
||||
ArrayList<Vector> mu_;
|
||||
ArrayList<Matrix> sigma_;
|
||||
Vector omega_;
|
||||
index_t number_of_gaussians_;
|
||||
index_t dimension_;
|
||||
|
||||
public:
|
||||
|
||||
MoGEM() {
|
||||
mu_.Init(0);
|
||||
sigma_.Init(0);
|
||||
}
|
||||
|
||||
~MoGEM() {
|
||||
}
|
||||
|
||||
void Init(index_t num_gauss, index_t dimension) {
|
||||
|
||||
// Initialize the private variables
|
||||
number_of_gaussians_ = num_gauss;
|
||||
dimension_ = dimension;
|
||||
|
||||
// Resize the ArrayList of Vectors and Matrices
|
||||
mu_.Resize(number_of_gaussians_);
|
||||
sigma_.Resize(number_of_gaussians_);
|
||||
}
|
||||
|
||||
void Init(datanode *mog_em_module) {
|
||||
|
||||
index_t num_gauss = fx_param_int_req(mog_em_module, "K");
|
||||
index_t dim = fx_param_int_req(mog_em_module, "D");
|
||||
Init(num_gauss, dim);
|
||||
}
|
||||
// The get functions
|
||||
|
||||
ArrayList<Vector>& mu() {
|
||||
return mu_;
|
||||
}
|
||||
|
||||
ArrayList<Matrix>& sigma() {
|
||||
return sigma_;
|
||||
}
|
||||
|
||||
Vector& omega() {
|
||||
return omega_;
|
||||
}
|
||||
|
||||
index_t number_of_gaussians() {
|
||||
return number_of_gaussians_;
|
||||
}
|
||||
|
||||
index_t dimension() {
|
||||
return dimension_;
|
||||
}
|
||||
|
||||
Vector& mu(index_t i) {
|
||||
return mu_[i] ;
|
||||
}
|
||||
|
||||
Matrix& sigma(index_t i) {
|
||||
return sigma_[i];
|
||||
}
|
||||
|
||||
double omega(index_t i) {
|
||||
return omega_.get(i);
|
||||
}
|
||||
|
||||
// The set functions
|
||||
|
||||
void set_mu(index_t i, Vector& mu) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(mu.length() == dimension());
|
||||
mu_[i].Copy(mu);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_mu(index_t i, index_t length, const double *mu) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(length == dimension());
|
||||
mu_[i].Copy(mu, length);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_sigma(index_t i, Matrix& sigma) {
|
||||
DEBUG_ASSERT(i < number_of_gaussians());
|
||||
DEBUG_ASSERT(sigma.n_rows() == dimension());
|
||||
DEBUG_ASSERT(sigma.n_cols() == dimension());
|
||||
sigma_[i].Copy(sigma);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_omega(Vector& omega) {
|
||||
DEBUG_ASSERT(omega.length() == number_of_gaussians());
|
||||
omega_.Copy(omega);
|
||||
return;
|
||||
}
|
||||
|
||||
void set_omega(index_t length, const double *omega) {
|
||||
DEBUG_ASSERT(length == number_of_gaussians());
|
||||
omega_.Copy(omega, length);
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* This function outputs the parameters of the model
|
||||
* to an arraylist of doubles
|
||||
*
|
||||
* @code
|
||||
* ArrayList<double> results;
|
||||
* mog.OutputResults(&results);
|
||||
* @endcode
|
||||
*/
|
||||
void OutputResults(ArrayList<double> *results) {
|
||||
|
||||
// Initialize the size of the output array
|
||||
(*results).Init(number_of_gaussians_ * (1 + dimension_*(1 + dimension_)));
|
||||
|
||||
// Copy values to the array from the private variables of the class
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
(*results)[i] = omega(i);
|
||||
for (index_t j = 0; j < dimension_; j++) {
|
||||
(*results)[number_of_gaussians_ + i*dimension_ + j] = mu(i).get(j);
|
||||
for (index_t k = 0; k < dimension_; k++) {
|
||||
(*results)[number_of_gaussians_*(1 + dimension_)
|
||||
+ i*dimension_*dimension_ + j*dimension_
|
||||
+ k] = sigma(i).get(j, k);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* This function prints the parameters of the model
|
||||
*
|
||||
* @code
|
||||
* mog.Display();
|
||||
* @endcode
|
||||
*/
|
||||
void Display(){
|
||||
|
||||
// Output the model parameters as the omega, mu and sigma
|
||||
printf(" Omega : [ ");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
printf("%lf ", omega(i));
|
||||
}
|
||||
printf("]\n");
|
||||
printf(" Mu : \n[");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
for (index_t j = 0; j < dimension_ ; j++) {
|
||||
printf("%lf ", mu(i).get(j));
|
||||
}
|
||||
printf(";");
|
||||
if (i == (number_of_gaussians_ - 1)) {
|
||||
printf("\b]\n");
|
||||
}
|
||||
}
|
||||
printf("Sigma : ");
|
||||
for (index_t i = 0; i < number_of_gaussians_; i++) {
|
||||
printf("\n[");
|
||||
for (index_t j = 0; j < dimension_ ; j++) {
|
||||
for(index_t k = 0; k < dimension_ ; k++) {
|
||||
printf("%lf ",sigma(i).get(j, k));
|
||||
}
|
||||
printf(";");
|
||||
}
|
||||
printf("\b]");
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
/**
|
||||
* This function calculates the parameters of the model
|
||||
* using the Maximum Likelihood function via the
|
||||
* Expectation Maximization (EM) Algorithm.
|
||||
*
|
||||
* @code
|
||||
* MoG mog;
|
||||
* Matrix data = "the data on which you want to fit the model";
|
||||
* ArrayList<double> results;
|
||||
* mog.ExpectationMaximization(data, &results);
|
||||
* @endcode
|
||||
*/
|
||||
void ExpectationMaximization(Matrix& data_points);
|
||||
|
||||
/**
|
||||
* This function computes the loglikelihood of model.
|
||||
* This function is used by the 'ExpectationMaximization'
|
||||
* function.
|
||||
*
|
||||
*/
|
||||
long double Loglikelihood(Matrix& data_points, ArrayList<Vector>& means,
|
||||
ArrayList<Matrix>& covars, Vector& weights);
|
||||
|
||||
/**
|
||||
* This function computes the k-means of the data and stores
|
||||
* the calculated means and covariances in the ArrayList
|
||||
* of Vectors and Matrices passed to it. It sets the weights
|
||||
* uniformly.
|
||||
*
|
||||
* This function is used to obtain a starting point for
|
||||
* the optimization
|
||||
*/
|
||||
void KMeans(Matrix& data, ArrayList<Vector> *means,
|
||||
ArrayList<Matrix> *covars, Vector *weights, index_t value_of_k);
|
||||
};
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,72 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file mog_l2e_main.cc
|
||||
*
|
||||
* This program test drives the L2 estimation
|
||||
* of a Gaussian Mixture model.
|
||||
*
|
||||
* PARAMETERS TO BE INPUT:
|
||||
*
|
||||
* --data
|
||||
* This is the file that contains the data on which
|
||||
* the model is to be fit
|
||||
*
|
||||
* --mog_em/K
|
||||
* This is the number of gaussians we want to fit
|
||||
* on the data, defaults to '1'
|
||||
*
|
||||
* --output
|
||||
* This file will contain the parameters estimated,
|
||||
* defaults to 'ouotput.csv'
|
||||
*
|
||||
*/
|
||||
|
||||
#include "mog.h"
|
||||
|
||||
|
||||
int main(int argc, char* argv[]) {
|
||||
|
||||
fx_init(argc, argv);
|
||||
|
||||
////// READING PARAMETERS AND LOADING DATA //////
|
||||
|
||||
const char *data_filename = fx_param_str_req(NULL, "data");
|
||||
|
||||
Matrix data_points;
|
||||
data::Load(data_filename, &data_points);
|
||||
|
||||
////// MIXTURE OF GAUSSIANS USING EM //////
|
||||
|
||||
MoGEM mog;
|
||||
|
||||
struct datanode* mog_em_module = fx_submodule(NULL, "mog_em", "mog_em");
|
||||
fx_param_int(mog_em_module, "K", 1);
|
||||
fx_format_param(mog_em_module, "D", "%d", data_points.n_rows());
|
||||
|
||||
////// Timing the initialization of the mixture model //////
|
||||
fx_timer_start(mog_em_module, "model_init");
|
||||
mog.Init(mog_em_module);
|
||||
fx_timer_stop(mog_em_module, "model_init");
|
||||
|
||||
////// Computing the parameters of the model using the EM algorithm //////
|
||||
ArrayList<double> results;
|
||||
|
||||
fx_timer_start(mog_em_module, "EM");
|
||||
mog.ExpectationMaximization(data_points);
|
||||
fx_timer_stop(mog_em_module, "EM");
|
||||
|
||||
mog.Display();
|
||||
mog.OutputResults(&results);
|
||||
|
||||
////// OUTPUT RESULTS //////
|
||||
|
||||
const char *output_filename = fx_param_str(NULL, "output", "output.csv");
|
||||
|
||||
FILE *output_file = fopen(output_filename, "w");
|
||||
|
||||
ot::Print(results, output_file);
|
||||
fclose(output_file);
|
||||
fx_done();
|
||||
|
||||
return 1;
|
||||
}
|
||||
@@ -0,0 +1,142 @@
|
||||
/**
|
||||
* @author Parikshit Ram (pram@cc.gatech.edu)
|
||||
* @file phi.h
|
||||
*
|
||||
* This file computes the Gaussian probability
|
||||
* density function
|
||||
*/
|
||||
#include "fastlib/fastlib.h"
|
||||
#include "fastlib/fastlib_int.h"
|
||||
#include <cmath>
|
||||
|
||||
/**
|
||||
* Calculates the multivariate Gaussian probability density function
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* Vector x, mean;
|
||||
* Matrix cov;
|
||||
* ....
|
||||
* long double f = phi(x, mean, cov);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(Vector& x , Vector& mean , Matrix& cov) {
|
||||
|
||||
long double det, f;
|
||||
double exponent;
|
||||
index_t dim;
|
||||
Matrix inv;
|
||||
Vector diff, tmp;
|
||||
|
||||
dim = x.length();
|
||||
la::InverseInit(cov, &inv);
|
||||
det = la::Determinant(cov);
|
||||
|
||||
if( det < 0){
|
||||
det = -det;
|
||||
}
|
||||
la::SubInit(mean,x,&diff);
|
||||
la::MulInit(inv, diff, &tmp);
|
||||
exponent = la::Dot(diff, tmp);
|
||||
long double tmp1, tmp2, tmp3;
|
||||
tmp1 = 1;
|
||||
tmp2 = dim;
|
||||
tmp2 = tmp2/2;
|
||||
tmp2 = pow((2*(math::PI)),tmp2);
|
||||
tmp1 = tmp1/tmp2;
|
||||
tmp3 = 1;
|
||||
tmp2 = sqrt(det);
|
||||
tmp3 = tmp3/tmp2;
|
||||
tmp2 = -exponent;
|
||||
tmp2 = tmp2 / 2;
|
||||
f = (tmp1*tmp3*exp(tmp2));
|
||||
|
||||
return f;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates the univariate Gaussian probability density function
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* double x, mean, var;
|
||||
* ....
|
||||
* long double f = phi(x, mean, var);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(double x, double mean, double var) {
|
||||
|
||||
long double f;
|
||||
|
||||
f = exp(-1.0*((x-mean)*(x-mean)/(2*var)))/sqrt(2*math::PI*var);
|
||||
return f;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates the multivariate Gaussian probability density function
|
||||
* and also the gradients with respect to the mean and the variance
|
||||
*
|
||||
* Example use:
|
||||
* @code
|
||||
* Vector x, mean, g_mean, g_cov;
|
||||
* ArrayList<Matrix> d_cov; // the dSigma
|
||||
* ....
|
||||
* long double f = phi(x, mean, cov, d_cov, &g_mean, &g_cov);
|
||||
* @endcode
|
||||
*/
|
||||
|
||||
long double phi(Vector& x, Vector& mean, Matrix& cov, ArrayList<Matrix>& d_cov, Vector *g_mean, Vector *g_cov){
|
||||
|
||||
long double det, f;
|
||||
double exponent;
|
||||
index_t dim;
|
||||
Matrix inv;
|
||||
Vector diff, tmp;
|
||||
|
||||
dim = x.length();
|
||||
la::InverseInit(cov, &inv);
|
||||
det = la::Determinant(cov);
|
||||
|
||||
if( det < 0){
|
||||
det = -det;
|
||||
}
|
||||
la::SubInit(mean,x,&diff);
|
||||
la::MulInit(inv, diff, &tmp);
|
||||
exponent = la::Dot(diff, tmp);
|
||||
long double tmp1, tmp2, tmp3;
|
||||
tmp1 = 1;
|
||||
tmp2 = dim;
|
||||
tmp2 = tmp2/2;
|
||||
tmp2 = pow((2*(math::PI)),tmp2);
|
||||
tmp1 = tmp1/tmp2;
|
||||
tmp3 = 1;
|
||||
tmp2 = sqrt(det);
|
||||
tmp3 = tmp3/tmp2;
|
||||
tmp2 = -exponent;
|
||||
tmp2 = tmp2 / 2;
|
||||
f = (tmp1*tmp3*exp(tmp2));
|
||||
|
||||
// Calculating the g_mean values which would be a (1 X dim) vector
|
||||
la::ScaleInit(f,tmp,g_mean);
|
||||
|
||||
// Calculating the g_cov values which would be a (1 X (dim*(dim+1)/2)) vector
|
||||
double *g_cov_tmp;
|
||||
g_cov_tmp = (double*)malloc(d_cov.size()*sizeof(double));
|
||||
for(index_t i = 0; i < d_cov.size(); i++){
|
||||
Vector tmp_d;
|
||||
Matrix inv_d;
|
||||
long double tmp_d_cov_d_r;
|
||||
|
||||
la::MulInit(d_cov[i],tmp,&tmp_d);
|
||||
tmp_d_cov_d_r = la::Dot(tmp_d,tmp);
|
||||
la::MulInit(inv,d_cov[i],&inv_d);
|
||||
for(index_t j = 0; j < dim; j++)
|
||||
tmp_d_cov_d_r += inv_d.get(j,j);
|
||||
g_cov_tmp[i] = f*tmp_d_cov_d_r/2;
|
||||
}
|
||||
g_cov->Copy(g_cov_tmp,d_cov.size());
|
||||
|
||||
return f;
|
||||
}
|
||||
@@ -51,7 +51,7 @@ void MoGEM::ExpectationMaximization(Matrix& data_points) {
|
||||
best_l = -INFTY;
|
||||
index_t restarts = 0;
|
||||
// performing 5 restarts and choosing the best from them
|
||||
while (restarts < 5) {
|
||||
while (restarts < 1) {
|
||||
|
||||
// assign initial values to 'mu', 'sig' and 'omega' using k-means
|
||||
KMeans(data_points, &mu_temp, &sigma_temp, &omega_temp, num_gauss);
|
||||
|
||||
Reference in New Issue
Block a user