test_mcpdf.cpp

一个很全的matlab工具包

// $Id: test_mcpdf.cpp,v 2.21 2004/11/10 12:21:15 kgadeyne Exp $

#include <pdf/mcpdf.h>
#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_ORIG 20000  // Exact samples drawn from analytic gaussian
#define NUM_SAMPLES_MC 10000 // Samples drawn from Monte Carlo approximation
#define DIMENSION 2

int row=0; // row index for for loops
int column=0; // column index for for loops
int sample=0; // index for sample number 

int main()
{
  cerr << "==================================================\n"
       << "Test of the MCPdf Class\n"
       << "=================================================="
       << endl;

  // Creating a (continuous) Gaussian with mu and sigma
  ColumnVector mu(DIMENSION);  SymmetricMatrix sigma(DIMENSION);
  mu = 1.5; sigma = 0.0;

  for (row=0; row < DIMENSION; row++)
    {
      for (column=0; column < DIMENSION; column++)
	  {
	    if (column == row) sigma[row][column]=1.0;
	  }
    }
  Gaussian My_first_gaussian(mu,sigma);
  cout << "Sampling experiment from a continuous Gaussian with MU =\n" << 
    mu << "and SIGMA =\n" << sigma << endl;

  // Generating (exact) Samples from this continuous density with 
  // cholesky sampling
  cerr << "Generating Samples from this gaussian with cholesky sampling" << endl;
  
  list<Sample<ColumnVector> > exact_samples(NUM_SAMPLES_ORIG);
  list<Sample<ColumnVector> >::iterator it;
  exact_samples = My_first_gaussian.SampleFrom(NUM_SAMPLES_ORIG,CHOLESKY,NULL);

  // Calculating the mean of these samples
  ColumnVector orig_mean(DIMENSION); ColumnVector orig_sum(DIMENSION); 
  orig_mean = 0.0; orig_sum = 0.0;

  for (it = exact_samples.begin(); it != exact_samples.end(); it++)
    {
      orig_sum += it->ValueGet();
    }
  orig_mean = orig_sum * (1.0 / NUM_SAMPLES_ORIG);
  
  // Calculating the covariance matrix
  ColumnVector orig_diff(DIMENSION); // Temporary storage of
				     // difference sample-mean
  Matrix difforig_sum(DIMENSION,DIMENSION); // Sum of the previous ones
  difforig_sum = 0.0;

  for (it = exact_samples.begin(); it != exact_samples.end(); it++)
    {
      orig_diff = it->ValueGet() - orig_mean;
      difforig_sum += orig_diff * orig_diff.transpose();
    }

  Matrix orig_Covariance(DIMENSION,DIMENSION);
  orig_Covariance = (Matrix) (difforig_sum / (NUM_SAMPLES_ORIG-1));

  // Part 2
  // Creating a MONTE CARLO density with these samples
  cerr << "Creating a MONTE CARLO density with these samples" << endl;
  
  MCPdf<ColumnVector> My_MonteCarlo_gaussian(NUM_SAMPLES_ORIG, DIMENSION);
  My_MonteCarlo_gaussian.ListOfSamplesSet(exact_samples);

  // Get Expected Value from this pdf
  ColumnVector EV(2);
  EV = My_MonteCarlo_gaussian.ExpectedValueGet();
  cerr << "EV of My_MonteCarlo_gaussian = \n" << EV << endl;
  
  // Now sample from this Monte Carlo distribution and make new stats
  cerr << "Sampling " << NUM_SAMPLES_MC << " samples from MONTE CARLO density with default method (O (N log N) )" << endl;
  list<Sample<ColumnVector> > mc_samples(NUM_SAMPLES_MC);
  mc_samples = My_MonteCarlo_gaussian.SampleFrom(NUM_SAMPLES_MC,DEFAULT,NULL);
  
  cerr << "And now faster (ordered uniform samples : O(N) )" << endl;
  mc_samples = My_MonteCarlo_gaussian.SampleFrom(NUM_SAMPLES_MC,RIPLEY,NULL);

  // Calculating the mean of the MC samples from the MC pdf
  cerr << "Calculating Mean and Covariance of the samples from the MC density" << endl;
  ColumnVector MC_mean(DIMENSION); ColumnVector MC_sum(DIMENSION); 
  MC_mean = 0.0 ; MC_sum = 0.0;
  
  for (it = mc_samples.begin(); it != mc_samples.end(); it++)
    { 
      MC_sum += it->ValueGet();
    }
  MC_mean = MC_sum * (1.0 / NUM_SAMPLES_MC);
  
  // Calculating the covariance matrix
  ColumnVector MC_diff(DIMENSION); // Temporary storage
  Matrix MC_diffsum(DIMENSION,DIMENSION);
  MC_diffsum = 0.0;

  for (it = mc_samples.begin(); it != mc_samples.end(); it++)
    {
      MC_diff = it->ValueGet() - MC_mean;
      MC_diffsum += MC_diff * MC_diff.transpose();
    }

  Matrix MC_Covariance(DIMENSION,DIMENSION);
  MC_Covariance = (Matrix) (MC_diffsum / (NUM_SAMPLES_MC-1));

  cout << "Number of original samples from continuous gaussian = " << NUM_SAMPLES_ORIG << endl;
  cout << "Number of MC Samples drawn from MC gaussian = " << NUM_SAMPLES_MC << "\n" << endl;

  cout << "Mean of original samples =\n" << orig_mean;
  cout << "Covariance of original samples =\n" <<  orig_Covariance << endl;

  cout << "mean of the MC experiment = \n" << MC_mean;
  cout << "Covariance Matrix of MC experiment =\n" << MC_Covariance  << endl;

  // Testing the ripley sampling method now!
  mc_samples = My_MonteCarlo_gaussian.SampleFrom(NUM_SAMPLES_MC,RIPLEY,NULL);
  // Calculate mean
  MC_mean = 0.0 ; MC_sum = 0.0;
  for (it = mc_samples.begin(); it != mc_samples.end(); it++)
    { 
      MC_sum += it->ValueGet();
    }
  MC_mean = MC_sum / NUM_SAMPLES_MC;

  // Calculating the covariance matrix
  MC_diffsum = 0.0;
  for (it = mc_samples.begin(); it != mc_samples.end(); it++)
    {
      MC_diff = it->ValueGet() - MC_mean;
      MC_diffsum += MC_diff * MC_diff.transpose();
    }
  MC_Covariance = (Matrix) (MC_diffsum / (NUM_SAMPLES_MC-1));

  cout << "RIPLEY Sampling experiment from a continuous Gaussian with MU =\n" << 
    mu << "and SIGMA =\n" << sigma << endl;

  cout << "Number of original samples = " << NUM_SAMPLES_ORIG << endl;
  cout << "Number of RIPLEY MC Samples = " << NUM_SAMPLES_MC << "\n" << endl;

  cout << "Mean of original samples =\n" << orig_mean;
  cout << "Covariance of original samples =\n" <<  orig_Covariance << endl;

  cout << "mean of the MC experiment = \n" << MC_mean;
  cout << "Covariance Matrix of MC experiment =\n" << MC_Covariance  << endl;
  return 0;
}