Statistics/src/RandomMultiGauss.cc

0001 //#include "CommonDet/DetUtilities/interface/DetExceptions.h"
0002 #include "CommonTools/Statistics/interface/RandomMultiGauss.h"
0003 #include "CLHEP/Random/RandGauss.h"
0004
0005 #include <cfloat>
0006 //
0007 // constructor with means and covariance
0008 //
0009 RandomMultiGauss::RandomMultiGauss(const AlgebraicVector& aVector, const AlgebraicSymMatrix& aMatrix)
0010     : theSize(aMatrix.num_row()), theMeans(aVector), theTriangle(theSize, theSize, 0) {
0011   //
0012   // Check consistency
0013   //
0014   if (theMeans.num_row() == theSize) {
0015     initialise(aMatrix);
0016   } else {
0017     //    throw DetLogicError("RandomMultiGauss: size of vector and matrix do not match");
0018     theMeans = AlgebraicVector(theSize, 0);
0019   }
0020 }
0021 //
0022 // constructor with covariance (mean = 0)
0023 //
0024 RandomMultiGauss::RandomMultiGauss(const AlgebraicSymMatrix& aMatrix)
0025     : theSize(aMatrix.num_row()), theMeans(theSize, 0), theTriangle(theSize, theSize, 0) {
0026   //
0027   initialise(aMatrix);
0028 }
0029 //
0030 // construct triangular matrix (Cholesky decomposition)
0031 //
0032 void RandomMultiGauss::initialise(const AlgebraicSymMatrix& aMatrix) {
0033   //
0034   // Cholesky decomposition with protection against empty rows/columns
0035   //
0036   for (int i1 = 0; i1 < theSize; i1++) {
0037     if (fabs(aMatrix[i1][i1]) < FLT_MIN)
0038       continue;
0039
0040     for (int i2 = i1; i2 < theSize; i2++) {
0041       if (fabs(aMatrix[i2][i2]) < FLT_MIN)
0042         continue;
0043
0044       double sum = aMatrix[i2][i1];
0045       for (int i3 = i1 - 1; i3 >= 0; i3--) {
0046         if (fabs(aMatrix[i3][i3]) < FLT_MIN)
0047           continue;
0048         sum -= theTriangle[i1][i3] * theTriangle[i2][i3];
0049       }
0050
0051       if (i1 == i2) {
0052         //
0053         // check for positive definite input matrix, but allow for effects
0054         // due to finite precision
0055         //
0056         if (sum <= 0) {
0057           //      if ( sum<-FLT_MIN )  throw DetLogicError("RandomMultiGauss: input matrix is not positive definite");
0058           sum = FLT_MIN;
0059         }
0060         theTriangle[i1][i1] = sqrt(sum);
0061       } else {
0062         theTriangle[i1][i2] = 0.;
0063         theTriangle[i2][i1] = sum / theTriangle[i1][i1];
0064       }
0065     }
0066   }
0067 }
0068 //
0069 // generate vector of random numbers
0070 //
0071 AlgebraicVector RandomMultiGauss::fire() {
0072   AlgebraicVector vRandom(theSize, 0);
0073   for (int i = 0; i < theSize; i++) {
0074     if (theTriangle[i][i] != 0)
0075       vRandom[i] = CLHEP::RandGauss::shoot();
0076   }
0077   return theTriangle * vRandom + theMeans;
0078 }