CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandMultiGauss.cc
Go to the documentation of this file.
1 // $Id: RandMultiGauss.cc,v 1.3 2003/08/13 20:00:13 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandMultiGauss ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 
10 // =======================================================================
11 // Mark Fischler - Created: 17th September 1998
12 // =======================================================================
13 
14 // Some theory about how to get the Multivariate Gaussian from a bunch
15 // of Gaussian deviates. For the purpose of this discussion, take mu = 0.
16 //
17 // We want an n-vector x with distribution (See table 28.1 of Review of PP)
18 //
19 // exp ( .5 * x.T() * S.inverse() * x )
20 // f(x;S) = ------------------------------------
21 // sqrt ( (2*pi)^n * S.det() )
22 //
23 // Suppose S = U * D * U.T() with U orthogonal ( U*U.T() = 1 ) and D diagonal.
24 // Consider a random n-vector y such that each element y(i)is distributed as a
25 // Gaussian with sigma = sqrt(D(i,i)). Then the distribution of y is the
26 // product of n Gaussains which can be written as
27 //
28 // exp ( .5 * y.T() * D.inverse() * y )
29 // f(y;D) = ------------------------------------
30 // sqrt ( (2*pi)^n * D.det() )
31 //
32 // Now take an n-vector x = U * y (or y = U.T() * x ). Then substituting,
33 //
34 // exp ( .5 * x * U * D.inverse() U.T() * x )
35 // f(x;D,U) = ------------------------------------------
36 // sqrt ( (2*pi)^n * D.det() )
37 //
38 // and this simplifies to the desired f(x;S) when we note that
39 // D.det() = S.det() and U * D.inverse() * U.T() = S.inverse()
40 //
41 // So the strategy is to diagonalize S (finding U and D), form y with each
42 // element a Gaussian random with sigma of sqrt(D(i,i)), and form x = U*y.
43 // (It turns out that the D information needed is the sigmas.)
44 // Since for moderate or large n the work needed to diagonalize S can be much
45 // greater than the work to generate n Gaussian variates, we save the U and
46 // sigmas for both the default S and the latest S value provided.
47 
50 #include <cmath> // for log()
51 
52 namespace CLHEP {
53 
54 // ------------
55 // Constructors
56 // ------------
57 
59  const HepVector & mu,
60  const HepSymMatrix & S )
61  : localEngine(&anEngine),
62  deleteEngine(false),
63  set(false),
64  nextGaussian(0.0)
65 {
66  if (S.num_row() != mu.num_row()) {
67  std::cerr << "In constructor of RandMultiGauss distribution: \n" <<
68  " Dimension of mu (" << mu.num_row() <<
69  ") does not match dimension of S (" << S.num_row() << ")\n";
70  std::cerr << "---Exiting to System\n";
71  exit(1);
72  }
73  defaultMu = mu;
74  defaultSigmas = HepVector(S.num_row());
75  prepareUsigmas (S, defaultU, defaultSigmas);
76 }
77 
79  const HepVector & mu,
80  const HepSymMatrix & S )
81  : localEngine(anEngine),
82  deleteEngine(true),
83  set(false),
84  nextGaussian(0.0)
85 {
86  if (S.num_row() != mu.num_row()) {
87  std::cerr << "In constructor of RandMultiGauss distribution: \n" <<
88  " Dimension of mu (" << mu.num_row() <<
89  ") does not match dimension of S (" << S.num_row() << ")\n";
90  std::cerr << "---Exiting to System\n";
91  exit(1);
92  }
93  defaultMu = mu;
94  defaultSigmas = HepVector(S.num_row());
95  prepareUsigmas (S, defaultU, defaultSigmas);
96 }
97 
99  : localEngine(&anEngine),
100  deleteEngine(false),
101  set(false),
102  nextGaussian(0.0)
103 {
104  defaultMu = HepVector(2,0);
105  defaultU = HepMatrix(2,1);
106  defaultSigmas = HepVector(2);
107  defaultSigmas(1) = 1.;
108  defaultSigmas(2) = 1.;
109 }
110 
112  : localEngine(anEngine),
113  deleteEngine(true),
114  set(false),
115  nextGaussian(0.0)
116 {
117  defaultMu = HepVector(2,0);
118  defaultU = HepMatrix(2,1);
119  defaultSigmas = HepVector(2);
120  defaultSigmas(1) = 1.;
121  defaultSigmas(2) = 1.;
122 }
123 
125  if ( deleteEngine ) delete localEngine;
126 }
127 
128 // ----------------------------
129 // prepareUsigmas()
130 // ----------------------------
131 
132 void RandMultiGauss::prepareUsigmas( const HepSymMatrix & S,
133  HepMatrix & U,
134  HepVector & sigmas ) {
135 
136  HepSymMatrix tempS ( S ); // Since diagonalize does not take a const s
137  // we have to copy S.
138 
139  U = diagonalize ( &tempS ); // S = U Sdiag U.T()
140  HepSymMatrix D = S.similarityT(U); // D = U.T() S U = Sdiag
141  for (int i = 1; i <= S.num_row(); i++) {
142  double s2 = D(i,i);
143  if ( s2 > 0 ) {
144  sigmas(i) = sqrt ( s2 );
145  } else {
146  std::cerr << "In RandMultiGauss distribution: \n" <<
147  " Matrix S is not positive definite. Eigenvalues are:\n";
148  for (int ixx = 1; ixx <= S.num_row(); ixx++) {
149  std::cerr << " " << D(ixx,ixx) << std::endl;
150  }
151  std::cerr << "---Exiting to System\n";
152  exit(1);
153  }
154  }
155 } // prepareUsigmas
156 
157 // -----------
158 // deviates()
159 // -----------
160 
161 HepVector RandMultiGauss::deviates ( const HepMatrix & U,
162  const HepVector & sigmas,
163  HepRandomEngine * engine,
164  bool& available,
165  double& next)
166 {
167  // Returns vector of gaussian randoms based on sigmas, rotated by U,
168  // with means of 0.
169 
170  int n = sigmas.num_row();
171  HepVector v(n); // The vector to be returned
172 
173  double r,v1,v2,fac;
174 
175  int i = 1;
176  if (available) {
177  v(1) = next;
178  i = 2;
179  available = false;
180  }
181 
182  while ( i <= n ) {
183  do {
184  v1 = 2.0 * engine->flat() - 1.0;
185  v2 = 2.0 * engine->flat() - 1.0;
186  r = v1*v1 + v2*v2;
187  } while ( r > 1.0 );
188  fac = sqrt(-2.0*log(r)/r);
189  v(i++) = v1*fac;
190  if ( i <= n ) {
191  v(i++) = v2*fac;
192  } else {
193  next = v2*fac;
194  available = true;
195  }
196  }
197 
198  for ( i = 1; i <= n; i++ ) {
199  v(i) *= sigmas(i);
200  }
201 
202  return U*v;
203 
204 } // deviates()
205 
206 // ---------------
207 // fire signatures
208 // ---------------
209 
211  // Returns a pair of unit normals, using the S and mu set in constructor,
212  // utilizing the engine belonging to this instance of RandMultiGauss.
213 
214  return defaultMu + deviates ( defaultU, defaultSigmas,
215  localEngine, set, nextGaussian );
216 
217 } // fire();
218 
219 
221 
222  HepMatrix U;
223  HepVector sigmas;
224 
225  if (mu.num_row() == S.num_row()) {
226  prepareUsigmas ( S, U, sigmas );
227  return mu + deviates ( U, sigmas, localEngine, set, nextGaussian );
228  } else {
229  std::cerr << "In firing RandMultiGauss distribution with explicit mu and S: \n"
230  << " Dimension of mu (" << mu.num_row() <<
231  ") does not match dimension of S (" << S.num_row() << ")\n";
232  std::cerr << "---Exiting to System\n";
233  exit(1);
234  }
235  return mu; // This line cannot be reached. But without returning
236  // some HepVector here, KCC 3.3 complains.
237 
238 } // fire(mu, S);
239 
240 
241 // --------------------
242 // fireArray signatures
243 // --------------------
244 
245 void RandMultiGauss::fireArray( const int size, HepVector* array ) {
246 
247  int i;
248  for (i = 0; i < size; ++i) {
249  array[i] = defaultMu + deviates ( defaultU, defaultSigmas,
250  localEngine, set, nextGaussian );
251  }
252 
253 } // fireArray ( size, vect )
254 
255 
256 void RandMultiGauss::fireArray( const int size, HepVector* array,
257  const HepVector& mu, const HepSymMatrix& S ) {
258 
259  // For efficiency, we diagonalize S once and generate all the vectors based
260  // on that U and sigmas.
261 
262  HepMatrix U;
263  HepVector sigmas;
264  HepVector mu_ (mu);
265 
266  if (mu.num_row() == S.num_row()) {
267  prepareUsigmas ( S, U, sigmas );
268  } else {
269  std::cerr <<
270  "In fireArray for RandMultiGauss distribution with explicit mu and S: \n"
271  << " Dimension of mu (" << mu.num_row() <<
272  ") does not match dimension of S (" << S.num_row() << ")\n";
273  std::cerr << "---Exiting to System\n";
274  exit(1);
275  }
276 
277  int i;
278  for (i=0; i<size; ++i) {
279  array[i] = mu_ + deviates(U, sigmas, localEngine, set, nextGaussian);
280  }
281 
282 } // fireArray ( size, vect, mu, S )
283 
284 // ----------
285 // operator()
286 // ----------
287 
289  return fire();
290 }
291 
292 HepVector RandMultiGauss::operator()
293  ( const HepVector& mu, const HepSymMatrix& S ) {
294  return fire(mu,S);
295 }
296 
297 
298 } // namespace CLHEP
CLHEP::RandMultiGauss::operator()
HepVector operator()()
Definition: RandMultiGauss.cc:288
CLHEP::RandMultiGauss::RandMultiGauss
RandMultiGauss(HepRandomEngine &anEngine, const HepVector &mu, const HepSymMatrix &S)
Definition: RandMultiGauss.cc:58
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
mu
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss RandGaussT RandGaussQ bins stepwise bins RandPoisson RandPoissonT mu< 100 N=50, 000, 000 ------- mu > RandPoissonQ mu< 100 N=50, 000, 000 -------(same as RandPoissonT) mu > RandGauss shoot() etc 2.5 units Validation tests applied and a very accurate series is substituted for the table method shoot() etc 1.7 units Validation tests applied and below about **a series approximation is but we have applied this test with N up to and never get anywhere near the rejection point Analytical considerations indicate that the validation test would not reject until O(10 **24) samples were inspected. ---------------------------------------------------------- 2. RandGeneral Method since we wish to have good resolution power even if just one of the mu values is it would be unwise to dial this down any further Validation for several values o mu)
Definition: validation.doc:263
D
Definition: excDblThrow.cc:17
CLHEP::detail::n
n
Definition: Ranlux64Engine.cc:85
CLHEP::HepMatrix
Definition: Matrix/CLHEP/Matrix/Matrix.h:209
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::diagonalize
HepMatrix diagonalize(HepSymMatrix *s)
Definition: MatrixLinear.cc:315
CLHEP
Definition: ClhepVersion.h:13
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::RandMultiGauss::~RandMultiGauss
virtual ~RandMultiGauss()
Definition: RandMultiGauss.cc:124
fire
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss RandGaussT RandGaussQ bins stepwise bins RandPoisson RandPoissonT mu< 100 N=50, 000, 000 ------- mu > RandPoissonQ mu< 100 N=50, 000, 000 -------(same as RandPoissonT) mu > RandGauss shoot() etc 2.5 units Validation tests applied and a very accurate series is substituted for the table method shoot() etc 1.7 units Validation tests applied and below about **a series approximation is but we have applied this test with N up to and never get anywhere near the rejection point Analytical considerations indicate that the validation test would not reject until O(10 **24) samples were inspected. ---------------------------------------------------------- 2. RandGeneral Method since we wish to have good resolution power even if just one of the mu values is it would be unwise to dial this down any further Validation for several values of but we have applied this with much higher N We validated the mai fire)() method for a variety of mu values between 0 and 100 at a level of 10
Definition: validation.doc:266
CLHEP::HepSymMatrix
Definition: Matrix/CLHEP/Matrix/SymMatrix.h:89
set
set(pkginclude_HEADERS itos.h) INSTALL(FILES $
Definition: Cast/Cast/CMakeLists.txt:2
CLHEP::HepSymMatrix::num_row
int num_row() const
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::HepSymMatrix::similarityT
HepSymMatrix similarityT(const HepMatrix &hm1) const
Definition: SymMatrix.cc:816
CLHEP::RandMultiGauss::fire
HepVector fire()
Definition: RandMultiGauss.cc:210
RandMultiGauss.h
CLHEP::RandMultiGauss::fireArray
void fireArray(const int size, HepVector *array)
Definition: RandMultiGauss.cc:245
CLHEP::HepVector
Definition: Matrix/CLHEP/Matrix/Vector.h:39
exit
it has advantages For I leave the ZMthrows but substitute I replaced ZMthrow with ZMthrowA in this package when will issue its message and exit(-1). Also
defs.h