// Author : Saket Soni, MTech (CS) - 2004, Sem I, Roll - mtc0420
// Topic : Probability Assignment 2
// Description : To emperically test general multivariate normal distribution.
///////////////
// Header files
///////////////
#include
#include
#include
#include
#include
#include
#include
///////////////
// Type declarations
///////////////
typedef long double ProbType;
typedef long Counter;
///////////////
// Constant declarations
///////////////
const ProbType PI = 3.141592653589793238462643382795;
const long MAXARR = 100;
const int FORMATWIDTH = 19;
const ProbType ALLOWEDERROR = 0.0000001;
const int FLOATINGPRECISION = 8;
///////////////
// Global variable declarations
///////////////
Counter NumOfIterations; // for number of iterations
ProbType CovMat[MAXARR][MAXARR]; // for true covariance matrix
ProbType LMat[MAXARR][MAXARR]; // for cholesky factor of covar matrix
ProbType Mean[MAXARR]; // for true mean vector
ProbType X[MAXARR]; // for Multivariate samples
ProbType N[MAXARR]; // for samples from N(0,1)
ProbType ECovMat[MAXARR][MAXARR]; // for emperical covariance matrix
ProbType EMean[MAXARR]; // for emperical mean vector
ProbType ErrorCovMat[MAXARR][MAXARR]; // for error in covariance matrix
ProbType ErrorMean[MAXARR]; // for error in mean vector
Counter MatSize; // for sample size, no of random vars
///////////////
// Function declarations
///////////////
inline ProbType Unifrnd(); // for generating uniform random nos. from (0,1]
void choldc(ProbType a[MAXARR][MAXARR], Counter n); // To decompose a given
// symmetric positive definite matrix into its factors
// using Cholesky Decompostion
void genSample(); // to generate a new multivariate sample in the vector X
void printMatrix(ProbType p[MAXARR][MAXARR], Counter n);
// to print a square matrix of size n
void printVector(ProbType p[MAXARR], Counter n);
// to print a vector of size n
///////////////
// Application Entry Point
///////////////
int main()
{
Counter i,j,k;
// initialize the random number generator
srand(time(0));
// set format specifications
cout.precision(FLOATINGPRECISION);
//
// reading values from the user ...
//
cout << "Enter the num of iterations required : "; // read how many iterations
cin >> NumOfIterations;
cout << "Enter the no of random variables required (> MatSize; // read sample size
cout << "Enter the co-variance matrix (what ever is ..." << endl
<< "inputed in lower triangular matrix is skipped) : " << endl;
for ( i = 0 ; i < MatSize ; i++ ) // read the covar matrix
{
for ( j = 0 ; j < MatSize ; j++ )
{
cin >> CovMat[i][j];
if ( i >= j )
LMat[i][j] = CovMat[i][j];
}
}
cout << "Enter the mean vector : " << endl; // read the mean vector
for ( i = 0 ; i < MatSize ; i++ )
cin >> Mean[i];
for ( i = 0 ; i < MatSize ; i++ ) // to copy upper triangular part of the
for ( j = MatSize-1 ; j > i ; j-- ) // covariance matrix onto the lower
CovMat[j][i] = CovMat[i][j]; // triangular part
//
// printing actual starting values ...
//
cout << "Covariance Matrix : " << endl; // print covar matrix
printMatrix(CovMat,MatSize);
cout << "Mean Vector : " << endl; // print mean vector
printVector(Mean,MatSize);
choldc(LMat,MatSize); // find choleskey decomposition of covar
// matrix and store it in matrix LMat
//
// Experiment starts from here ...
//
for ( k = 0 ; k < NumOfIterations ; k++ ) // Iterate NumOfIteratins times.
{
genSample(); // To generate new multivariate sample in X;
for ( i = 0 ; i < MatSize ; i++ )
{
EMean[i] = EMean[i] + X[i]; // keep running total for finding E(X)
for ( j = i ; j < MatSize ; j++ ) // keep running total for findind E(X*Y)
ECovMat[i][j] = ECovMat[i][j] + X[i]*X[j];
}
}
//
// Experiment finishes here ...
//
//
// Computing results ...
//
for ( i = 0 ; i < MatSize ; i++ )
{
EMean[i] = EMean[i]/ProbType(NumOfIterations); // find E(X) = sum(X)/total
for ( j = i ; j < MatSize ; j++ ) // find E(X^2) = sum(X*Y)/total
ECovMat[i][j] = (ECovMat[i][j])/ProbType(NumOfIterations);
}
for ( i = 0 ; i < MatSize ; i++ )
for ( j = i ; j < MatSize ; j++ ) // find Covar(XY) = E(X*Y)-E(X)*E(Y)
ECovMat[i][j] = (ECovMat[i][j]) - ((EMean[i])*(EMean[j]));
for ( i = 0 ; i < MatSize ; i++ ) // to copy upper triangular of emperical
for ( j = MatSize-1 ; j > i ; j-- ) // covar matrix onto the lower triangular part
ECovMat[j][i] = ECovMat[i][j];
//
// Printing the emperical results ...
//
cout << "Emperical Covariance Matrix : " << endl;
printMatrix(ECovMat,MatSize); // print the emperical covar matrix
cout << "Emperical Mean Vector : " << endl;
printVector(EMean,MatSize); // print the emperical mean vector
//
// Finding error in the emperical test...
//
for ( i = 0 ; i < MatSize ; i++ ) {
ErrorMean[i] = Mean[i] - EMean[i];
ErrorMean[i] = ErrorMean[i]/Mean[i]; // relative error in mean
if ( fabs(ErrorMean[i]) <= ALLOWEDERROR ) // if error is less than allowed error
ErrorMean[i] = 0; // then set error to zero
for ( j = 0 ; j < MatSize ; j++ ) {
ErrorCovMat[i][j] = CovMat[i][j] - ECovMat[i][j];
ErrorCovMat[i][j] = ErrorCovMat[i][j]/CovMat[i][j]; // relative error in covar mat
if ( fabs(ErrorCovMat[i][j]) <= ALLOWEDERROR ) // if error is less than allowed
ErrorCovMat[i][j] = 0; // error then set error to zero
}
}
//
// Printing the error found ...
///
cout << "Error in Covariance Matrix : " << endl;
printMatrix(ErrorCovMat,MatSize);
cout << "Error in Emperical Mean Vector : " << endl;
printVector(ErrorMean,MatSize);
return 0; // return successfully.
}
///////////////
// Fnction Definitions
///////////////
inline ProbType unifrnd()
{
return (rand() / ProbType(RAND_MAX)); // generate random nos from uniform (0,1]
}
void choldc(ProbType a[MAXARR][MAXARR], Counter n)
{
ProbType p[MAXARR], sum;
Counter i,j,k;
for ( i = 0 ; i < n ; i++ )
{
for ( j = i ; j < n ; j++ )
{
sum = CovMat[i][j]; // Initialize L(j,i) with Covar(i,j)
for ( k = i-1 ; k >= 0 ; k-- )
sum = sum - a[i][k] * a[j][k];
if (i == j) // if i == j then diagonal element
{
if (sum <= 0.0) // if sum <=0 then not +ve definite
{
cerr << endl << endl
<< "Error : Covariance Matrix is not positive definite !"
<< endl << endl;
exit(1);
}
p[i]=sqrt(sum);
}
else
a[j][i] = sum / p[i]; // for non diagonal elements
}
}
for (i = 0 ; i < n ; i++ ) // erase upper diagonal part
for ( j = n-1 ; j > i ; j-- )
a[i][j] = 0;
for ( i = 0 ; i < n ; i++ ) // copy diagonal elements
a[i][i] = p[i];
}
void genSample()
{
static ProbType u1, u2, l1, l2;
static Counter i,j;
for ( i = 0 ; i < MatSize ; )
{
u1 = unifrnd(); // generate uniform sample u1, u2
u2 = unifrnd();
if ( (u1 > 1e-32 || u1 < -1e-32) )
{
l1 = sqrt(-2.0*log(u1));
l2 = 2.0*PI*u2;
N[i++] = l1*cos(l2); // generate normal sample N[i] & N[i+1] from N(0,1)
N[i++] = l1*sin(l2);
}
}
for ( i = 0 ; i < MatSize ; i++ ) { // X = L * N + Mean
X[i] = 0.0;
for ( j = 0 ; j <= i ; j++ )
X[i] = X[i] + LMat[i][j]*N[j];
X[i] = X[i] + Mean[i]; // multivariate sample is now ready in X
}
}
void printMatrix(ProbType p[MAXARR][MAXARR], Counter n)
{
Counter i, j;
for ( i = 0, j ; i < MatSize ; i++ ) { // print the matrix row by row
for ( j = 0 ; j < MatSize ; j++ )
cout << setw(FORMATWIDTH) << p[i][j];
cout << endl;
}
}
void printVector(ProbType p[MAXARR], Counter n)
{
for ( Counter i = 0 ; i < MatSize ; i++ ) // prints the transpose of the col vector
cout << setw(FORMATWIDTH) << p[i];
cout << endl;
}