test_gaussian.cpp

一个很全的matlab工具包

// $Id: test_gaussian.cpp,v 2.22 2005/01/28 16:05:36 kgadeyne Exp $
// Copyright (C) 2002 Klaas Gadeyne <klaas dot gadeyne at mech dot kuleuven dot ac dot be>
//  
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//  
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//  
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//  

// TEST program to check the Gaussian class
#include <pdf/gaussian.h>
#include <sample/sample.h>
#include <matrix_wrapper.h>
#include <vector_wrapper.h>

#include <iostream>
#include <fstream>

using namespace std;
using namespace BFL;
using namespace MatrixWrapper;


#define NUM_SAMPLES 500
#define DIMENSION 2

#define PLOT_LIMIT 2
#define PLOT_INTERVAL 0.1

#define MU_0 2.8
#define MU_1 3.4
#define SIGMA 0.01

ofstream functionfile("data_visualisation/function.dat", ios::out);
ofstream samplesfile("data_visualisation/samples.dat", ios::out);

int main()
{
  cerr << "==================================================\n"
       << "Test of the Gaussian Class\n"
       << "=================================================="
       << endl;
  int rows;
  
  // Creating a Gaussian with mu and sigma
  ColumnVector mu(DIMENSION);  SymmetricMatrix sigma(DIMENSION);
  mu[0] = MU_0; mu[1] = MU_1; sigma = 0.0;
  for (rows=0; rows < DIMENSION; rows++){ sigma[rows][rows]=SIGMA; }
  cerr << "Creating the Gaussian with mu = \n" << 
    mu << "and sigma = \n" << sigma << endl;
  Gaussian My_first_gaussian(mu,sigma);

  // Creating a second Gaussian with mu2 and sigma2
#define MU2 3
#define SIGMA2 0.06
  ColumnVector mu2(DIMENSION);  SymmetricMatrix sigma2(DIMENSION);
  mu2 = MU2; sigma2 = 0.0;
  for (rows=0; rows < DIMENSION; rows++){ sigma2[rows][rows]=SIGMA2; }
  sigma2[1][0] = 0.05;
  cerr << "Creating a second Gaussian with mu = \n" << 
    mu2 << "and sigma = \n" << sigma2 << endl;
  Gaussian My_second_gaussian(mu2,sigma2);

  // Plotting the function
//   ColumnVector state(DIMENSION); double pdf_value; 
//   rows = 0;
//   if (functionfile)
//     {
//       for (double state_x = (0-PLOT_LIMIT) ; 
// 	   state_x <= PLOT_LIMIT; 
// 	   state_x += PLOT_INTERVAL)
// 	{
// 	  for (double state_y = (0-PLOT_LIMIT) ; 
// 	       state_y <= PLOT_LIMIT; 
// 	       state_y += PLOT_INTERVAL)
// 	    {
// 	      for (double state_z = (0-PLOT_LIMIT) ; 
// 		   state_z <= PLOT_LIMIT; 
// 		   state_z += PLOT_INTERVAL)
// 		{
// 		  state[0] = state_x; state[1] = state_y; state[2] = state_z; 
// 		  pdf_value = My_first_gaussian.ProbabilityGet(state);
// 		  functionfile << state[0] << " " << state[1] << " " << state[2] << " " << pdf_value << "\n";
// 		}
// 	    }
// 	}
//       functionfile.close();
//     }
//   else
//     { 
//       cerr << "Error writing to file\n";
//     }

  // Generating Samples with cholesky sampling
  cerr << "Generating " << NUM_SAMPLES << " Samples with cholesky sampling" << endl;
  list<Sample<ColumnVector> > los(NUM_SAMPLES);
  list<Sample<ColumnVector> >::iterator it;
  los = My_first_gaussian.SampleFrom(NUM_SAMPLES,CHOLESKY,NULL);

  cerr << "Generating " << NUM_SAMPLES << " Samples with cholesky sampling from second Gaussian" << endl;
  list<Sample<ColumnVector> > los2(NUM_SAMPLES/2);
  list<Sample<ColumnVector> >::iterator it2;
  los2 = My_second_gaussian.SampleFrom(NUM_SAMPLES,CHOLESKY,NULL);

  // Plot the list of samples
  if (samplesfile)
    {
      for (it = los.begin(); it != los.end(); it++)
	{
	  samplesfile << (it->ValueGet()).transpose();
	}

      for (it2 = los2.begin(); it2 != los2.end(); it2++)
	{
	  samplesfile << (it2->ValueGet()).transpose();
	}

      samplesfile.close();
    }
  
  // Calculating the mean
  cerr << "Calculating Sample Mean and Covariance" << endl;
  
  ColumnVector mean(DIMENSION); ColumnVector sum(DIMENSION); 
  for (int k = 0; k < DIMENSION; k++){ mean[k] = sum[k] = 0; }
  
  for (it = los.begin(); it != los.end(); it++)
    {
      sum += it->ValueGet();
    }
  mean = sum * (1.0 / NUM_SAMPLES);
  
  // Calculating the covariance matrix
  ColumnVector diff(DIMENSION); // Temporary storage
  Matrix diffsum(DIMENSION,DIMENSION);
  diffsum = 0.0;
  
  for (it = los.begin(); it != los.end(); it++)
    {
      diff = it->ValueGet() - mean;
      diffsum += diff * diff.transpose();
    }

  Matrix Covariance(DIMENSION,DIMENSION);
  Covariance = (Matrix) (diffsum / (NUM_SAMPLES-1));

  cout << "mean = \n" << mean << endl;
  cout << "Covariance Matrix = \n " << Covariance  << endl;

  // Trying to change matrix and covariance
  cerr << "Trying to change mean and covariance of gaussian" << endl;
  
  mu = 3; sigma = 1.0;
  for (rows=0; rows < DIMENSION; rows++){ sigma[rows][rows]=2.0; }

  My_first_gaussian.CovarianceSet(sigma);
  My_first_gaussian.ExpectedValueSet(mu);
  cout << "Covariance:\n" << My_first_gaussian.CovarianceGet() << endl;
  cout << "Mu:\n" << My_first_gaussian.ExpectedValueGet() << endl;
  cout << "Done" << endl;

  // Test the copy constructor and assignment operator
  cout << "Test the copy constructor and assignment operator" << endl;
  
  Gaussian My_third_gaussian(My_first_gaussian);
  My_third_gaussian = My_first_gaussian;
  
  return 1;
}