CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandMultiGauss.cc,v 1.3 2003/08/13 20:00:13 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandMultiGauss --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // Mark Fischler - Created: 17th September 1998 00012 // ======================================================================= 00013 00014 // Some theory about how to get the Multivariate Gaussian from a bunch 00015 // of Gaussian deviates. For the purpose of this discussion, take mu = 0. 00016 // 00017 // We want an n-vector x with distribution (See table 28.1 of Review of PP) 00018 // 00019 // exp ( .5 * x.T() * S.inverse() * x ) 00020 // f(x;S) = ------------------------------------ 00021 // sqrt ( (2*pi)^n * S.det() ) 00022 // 00023 // Suppose S = U * D * U.T() with U orthogonal ( U*U.T() = 1 ) and D diagonal. 00024 // Consider a random n-vector y such that each element y(i)is distributed as a 00025 // Gaussian with sigma = sqrt(D(i,i)). Then the distribution of y is the 00026 // product of n Gaussains which can be written as 00027 // 00028 // exp ( .5 * y.T() * D.inverse() * y ) 00029 // f(y;D) = ------------------------------------ 00030 // sqrt ( (2*pi)^n * D.det() ) 00031 // 00032 // Now take an n-vector x = U * y (or y = U.T() * x ). Then substituting, 00033 // 00034 // exp ( .5 * x * U * D.inverse() U.T() * x ) 00035 // f(x;D,U) = ------------------------------------------ 00036 // sqrt ( (2*pi)^n * D.det() ) 00037 // 00038 // and this simplifies to the desired f(x;S) when we note that 00039 // D.det() = S.det() and U * D.inverse() * U.T() = S.inverse() 00040 // 00041 // So the strategy is to diagonalize S (finding U and D), form y with each 00042 // element a Gaussian random with sigma of sqrt(D(i,i)), and form x = U*y. 00043 // (It turns out that the D information needed is the sigmas.) 00044 // Since for moderate or large n the work needed to diagonalize S can be much 00045 // greater than the work to generate n Gaussian variates, we save the U and 00046 // sigmas for both the default S and the latest S value provided. 00047 00048 #include "CLHEP/RandomObjects/RandMultiGauss.h" 00049 #include "CLHEP/RandomObjects/defs.h" 00050 #include <cmath> // for log() 00051 00052 namespace CLHEP { 00053 00054 // ------------ 00055 // Constructors 00056 // ------------ 00057 00058 RandMultiGauss::RandMultiGauss( HepRandomEngine & anEngine, 00059 const HepVector & mu, 00060 const HepSymMatrix & S ) 00061 : localEngine(&anEngine), 00062 deleteEngine(false), 00063 set(false), 00064 nextGaussian(0.0) 00065 { 00066 if (S.num_row() != mu.num_row()) { 00067 std::cerr << "In constructor of RandMultiGauss distribution: \n" << 00068 " Dimension of mu (" << mu.num_row() << 00069 ") does not match dimension of S (" << S.num_row() << ")\n"; 00070 std::cerr << "---Exiting to System\n"; 00071 exit(1); 00072 } 00073 defaultMu = mu; 00074 defaultSigmas = HepVector(S.num_row()); 00075 prepareUsigmas (S, defaultU, defaultSigmas); 00076 } 00077 00078 RandMultiGauss::RandMultiGauss( HepRandomEngine * anEngine, 00079 const HepVector & mu, 00080 const HepSymMatrix & S ) 00081 : localEngine(anEngine), 00082 deleteEngine(true), 00083 set(false), 00084 nextGaussian(0.0) 00085 { 00086 if (S.num_row() != mu.num_row()) { 00087 std::cerr << "In constructor of RandMultiGauss distribution: \n" << 00088 " Dimension of mu (" << mu.num_row() << 00089 ") does not match dimension of S (" << S.num_row() << ")\n"; 00090 std::cerr << "---Exiting to System\n"; 00091 exit(1); 00092 } 00093 defaultMu = mu; 00094 defaultSigmas = HepVector(S.num_row()); 00095 prepareUsigmas (S, defaultU, defaultSigmas); 00096 } 00097 00098 RandMultiGauss::RandMultiGauss( HepRandomEngine & anEngine ) 00099 : localEngine(&anEngine), 00100 deleteEngine(false), 00101 set(false), 00102 nextGaussian(0.0) 00103 { 00104 defaultMu = HepVector(2,0); 00105 defaultU = HepMatrix(2,1); 00106 defaultSigmas = HepVector(2); 00107 defaultSigmas(1) = 1.; 00108 defaultSigmas(2) = 1.; 00109 } 00110 00111 RandMultiGauss::RandMultiGauss( HepRandomEngine * anEngine ) 00112 : localEngine(anEngine), 00113 deleteEngine(true), 00114 set(false), 00115 nextGaussian(0.0) 00116 { 00117 defaultMu = HepVector(2,0); 00118 defaultU = HepMatrix(2,1); 00119 defaultSigmas = HepVector(2); 00120 defaultSigmas(1) = 1.; 00121 defaultSigmas(2) = 1.; 00122 } 00123 00124 RandMultiGauss::~RandMultiGauss() { 00125 if ( deleteEngine ) delete localEngine; 00126 } 00127 00128 // ---------------------------- 00129 // prepareUsigmas() 00130 // ---------------------------- 00131 00132 void RandMultiGauss::prepareUsigmas( const HepSymMatrix & S, 00133 HepMatrix & U, 00134 HepVector & sigmas ) { 00135 00136 HepSymMatrix tempS ( S ); // Since diagonalize does not take a const s 00137 // we have to copy S. 00138 00139 U = diagonalize ( &tempS ); // S = U Sdiag U.T() 00140 HepSymMatrix D = S.similarityT(U); // D = U.T() S U = Sdiag 00141 for (int i = 1; i <= S.num_row(); i++) { 00142 double s2 = D(i,i); 00143 if ( s2 > 0 ) { 00144 sigmas(i) = sqrt ( s2 ); 00145 } else { 00146 std::cerr << "In RandMultiGauss distribution: \n" << 00147 " Matrix S is not positive definite. Eigenvalues are:\n"; 00148 for (int ixx = 1; ixx <= S.num_row(); ixx++) { 00149 std::cerr << " " << D(ixx,ixx) << std::endl; 00150 } 00151 std::cerr << "---Exiting to System\n"; 00152 exit(1); 00153 } 00154 } 00155 } // prepareUsigmas 00156 00157 // ----------- 00158 // deviates() 00159 // ----------- 00160 00161 HepVector RandMultiGauss::deviates ( const HepMatrix & U, 00162 const HepVector & sigmas, 00163 HepRandomEngine * engine, 00164 bool& available, 00165 double& next) 00166 { 00167 // Returns vector of gaussian randoms based on sigmas, rotated by U, 00168 // with means of 0. 00169 00170 int n = sigmas.num_row(); 00171 HepVector v(n); // The vector to be returned 00172 00173 double r,v1,v2,fac; 00174 00175 int i = 1; 00176 if (available) { 00177 v(1) = next; 00178 i = 2; 00179 available = false; 00180 } 00181 00182 while ( i <= n ) { 00183 do { 00184 v1 = 2.0 * engine->flat() - 1.0; 00185 v2 = 2.0 * engine->flat() - 1.0; 00186 r = v1*v1 + v2*v2; 00187 } while ( r > 1.0 ); 00188 fac = sqrt(-2.0*log(r)/r); 00189 v(i++) = v1*fac; 00190 if ( i <= n ) { 00191 v(i++) = v2*fac; 00192 } else { 00193 next = v2*fac; 00194 available = true; 00195 } 00196 } 00197 00198 for ( i = 1; i <= n; i++ ) { 00199 v(i) *= sigmas(i); 00200 } 00201 00202 return U*v; 00203 00204 } // deviates() 00205 00206 // --------------- 00207 // fire signatures 00208 // --------------- 00209 00210 HepVector RandMultiGauss::fire() { 00211 // Returns a pair of unit normals, using the S and mu set in constructor, 00212 // utilizing the engine belonging to this instance of RandMultiGauss. 00213 00214 return defaultMu + deviates ( defaultU, defaultSigmas, 00215 localEngine, set, nextGaussian ); 00216 00217 } // fire(); 00218 00219 00220 HepVector RandMultiGauss::fire( const HepVector& mu, const HepSymMatrix& S ) { 00221 00222 HepMatrix U; 00223 HepVector sigmas; 00224 00225 if (mu.num_row() == S.num_row()) { 00226 prepareUsigmas ( S, U, sigmas ); 00227 return mu + deviates ( U, sigmas, localEngine, set, nextGaussian ); 00228 } else { 00229 std::cerr << "In firing RandMultiGauss distribution with explicit mu and S: \n" 00230 << " Dimension of mu (" << mu.num_row() << 00231 ") does not match dimension of S (" << S.num_row() << ")\n"; 00232 std::cerr << "---Exiting to System\n"; 00233 exit(1); 00234 } 00235 return mu; // This line cannot be reached. But without returning 00236 // some HepVector here, KCC 3.3 complains. 00237 00238 } // fire(mu, S); 00239 00240 00241 // -------------------- 00242 // fireArray signatures 00243 // -------------------- 00244 00245 void RandMultiGauss::fireArray( const int size, HepVector* array ) { 00246 00247 int i; 00248 for (i = 0; i < size; ++i) { 00249 array[i] = defaultMu + deviates ( defaultU, defaultSigmas, 00250 localEngine, set, nextGaussian ); 00251 } 00252 00253 } // fireArray ( size, vect ) 00254 00255 00256 void RandMultiGauss::fireArray( const int size, HepVector* array, 00257 const HepVector& mu, const HepSymMatrix& S ) { 00258 00259 // For efficiency, we diagonalize S once and generate all the vectors based 00260 // on that U and sigmas. 00261 00262 HepMatrix U; 00263 HepVector sigmas; 00264 HepVector mu_ (mu); 00265 00266 if (mu.num_row() == S.num_row()) { 00267 prepareUsigmas ( S, U, sigmas ); 00268 } else { 00269 std::cerr << 00270 "In fireArray for RandMultiGauss distribution with explicit mu and S: \n" 00271 << " Dimension of mu (" << mu.num_row() << 00272 ") does not match dimension of S (" << S.num_row() << ")\n"; 00273 std::cerr << "---Exiting to System\n"; 00274 exit(1); 00275 } 00276 00277 int i; 00278 for (i=0; i<size; ++i) { 00279 array[i] = mu_ + deviates(U, sigmas, localEngine, set, nextGaussian); 00280 } 00281 00282 } // fireArray ( size, vect, mu, S ) 00283 00284 // ---------- 00285 // operator() 00286 // ---------- 00287 00288 HepVector RandMultiGauss::operator()() { 00289 return fire(); 00290 } 00291 00292 HepVector RandMultiGauss::operator() 00293 ( const HepVector& mu, const HepSymMatrix& S ) { 00294 return fire(mu,S); 00295 } 00296 00297 00298 } // namespace CLHEP