RandomMultiGauss Class Reference

#include <RandomMultiGauss.h>

Public Member Functions
AlgebraicVector	fire ()
	RandomMultiGauss (const AlgebraicVector &aVector, const AlgebraicSymMatrix &aMatrix)
	RandomMultiGauss (const AlgebraicSymMatrix &aMatrix)
	~RandomMultiGauss ()

Private Member Functions
void	initialise (const AlgebraicSymMatrix &)

Private Attributes
AlgebraicVector	theMeans
int	theSize
AlgebraicMatrix	theTriangle

Detailed Description

Vector of random numbers according to covariance matrix. Generates vectors of random numbers given a vector of mean values (optional) and a covariance matrix. Will accept empty rows/columns in the input matrix. Uses CLHEP::RandGauss with default engine for generation.

Definition at line 14 of file RandomMultiGauss.h.

Constructor & Destructor Documentation

RandomMultiGauss::RandomMultiGauss	(	const AlgebraicVector &	aVector,
		const AlgebraicSymMatrix &	aMatrix
	)

constructor with explicit vector of mean values

Definition at line 9 of file RandomMultiGauss.cc.

References initialise(), theMeans, and theSize.

                                                                                                     :
  theSize(aMatrix.num_row()),
  theMeans(aVector),
  theTriangle(theSize,theSize,0) {
  //
  // Check consistency
  //
  if ( theMeans.num_row() == theSize ) {
    initialise(aMatrix);
  }
  else {
//    throw DetLogicError("RandomMultiGauss: size of vector and matrix do not match");
    theMeans = AlgebraicVector(theSize,0);
  }
}

RandomMultiGauss::RandomMultiGauss

(

const AlgebraicSymMatrix &

aMatrix

)

constructor with covariance matrix only (all means = 0)

Definition at line 27 of file RandomMultiGauss.cc.

References initialise().

                                                                     :
  theSize(aMatrix.num_row()),
  theMeans(theSize,0),
  theTriangle(theSize,theSize,0) {
  //
  initialise(aMatrix);
}

RandomMultiGauss::~RandomMultiGauss ( )

inline

Definition at line 23 of file RandomMultiGauss.h.

References fire(), and initialise().

{}

Member Function Documentation

AlgebraicVector RandomMultiGauss::fire

(

)

Generation of a vector of random numbers.

Definition at line 74 of file RandomMultiGauss.cc.

References i, theMeans, theSize, and theTriangle.

Referenced by ~RandomMultiGauss().

                                       {
  AlgebraicVector vRandom(theSize,0);
  for ( int i=0; i<theSize; i++ ) {
    if ( theTriangle[i][i]!=0 ) 
      vRandom[i] = CLHEP::RandGauss::shoot();
  }
  return theTriangle*vRandom+theMeans;
}

void RandomMultiGauss::initialise ( const AlgebraicSymMatrix &  aMatrix )

private

generation of the cholesky decomposition

Definition at line 37 of file RandomMultiGauss.cc.

References mathSSE::sqrt(), theSize, and theTriangle.

Referenced by RandomMultiGauss(), and ~RandomMultiGauss().

                                                                    {
  //
  // Cholesky decomposition with protection against empty rows/columns
  //
  for ( int i1=0; i1<theSize; i1++ ) {
    if ( fabs(aMatrix[i1][i1])<FLT_MIN )  continue;

    for ( int i2=i1; i2<theSize; i2++ ) {
      if ( fabs(aMatrix[i2][i2])<FLT_MIN )  continue;

      double sum = aMatrix[i2][i1];
      for ( int i3=i1-1; i3>=0; i3-- ) {
    if ( fabs(aMatrix[i3][i3])<FLT_MIN )  continue;
    sum -= theTriangle[i1][i3]*theTriangle[i2][i3];
      }
      
      if ( i1==i2 ) {
    //
    // check for positive definite input matrix, but allow for effects
    // due to finite precision
    //
    if ( sum<=0 ) {
//    if ( sum<-FLT_MIN )  throw DetLogicError("RandomMultiGauss: input matrix is not positive definite");
      sum = FLT_MIN;
    }
    theTriangle[i1][i1] = sqrt(sum);
      }
      else {
    theTriangle[i1][i2] = 0.;
    theTriangle[i2][i1] = sum / theTriangle[i1][i1];
      }
    }
  }
}

Member Data Documentation

AlgebraicVector RandomMultiGauss::theMeans

private

Definition at line 35 of file RandomMultiGauss.h.

Referenced by fire(), and RandomMultiGauss().

int RandomMultiGauss::theSize

private

Definition at line 34 of file RandomMultiGauss.h.

Referenced by fire(), initialise(), and RandomMultiGauss().

AlgebraicMatrix RandomMultiGauss::theTriangle

private

Definition at line 36 of file RandomMultiGauss.h.

Referenced by fire(), and initialise().