#define NSAMPLE 1e6

void RandomVars()
{
  gSystem->Load("libMathCore");
  gSystem->Load("libMathMore");
  gSystem->Load("libHist");
  gSystem->Load("libGraf");
  gSystem->Load("libGpad");
  gSystem->Load("libMatrix");

  TRandom3 *ranx = new TRandom3(1234567);
  TRandom3 *rany = new TRandom3(7654321);
  TPrincipal *pr = new TPrincipal(2, "ND");

  TMatrixD MeanX(2,1);
  TMatrixD CovX(2,2), TX(2,1), TSX(1,2);
  TMatrixD **X = new TMatrixD[NSAMPLE];
  Option_t *opt = "F";

  MeanX.Zero();
  CovX.Zero();
  TSX.Zero();

  Long_t it=-1;
  Double_t x[2];
  while(it++<NSAMPLE){
    X[it] = new TMatrixD(2,1);
    x[0] = ranx->Gaus(0,5);
    x[1] = rany->Gaus(0,1);
    X[it]->SetMatrixArray(x,opt);
    MeanX = MeanX + *X[it];
    pr->AddRow(x);
  }
  pr->MakePrincipals();
  cout << "FROM TPRINCIPAL ::::::";
  pr->GetCovarianceMatrix()->Print();

  MeanX *= 1.0/NSAMPLE;
  for(it=0;it<NSAMPLE;it++) {
    TX = *X[it] - MeanX;
    TSX.Transpose(*X[it] - MeanX);
    TMatrixD samp = TX*TSX;
    CovX = CovX + samp;
  }
  CovX *= 1.0/NSAMPLE;
  CovX *= 1.0/CovX.Determinant();
  cout << "COVARIANCE FROM DATA :::::::" ;
  CovX.Print();
  
  delete pr;
  delete ranx;
  delete rany;
  for(it=0;it<NSAMPLE;it++)
    delete X[it];
  //delete [] X; 
}