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
52namespace 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
132void 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
161HepVector 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
245void 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
256void 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
292HepVector RandMultiGauss::operator()
293 ( const HepVector& mu, const HepSymMatrix& S ) {
294 return fire(mu,S);
295}
296
297
298} // namespace CLHEP
int num_row() const
HepSymMatrix similarityT(const HepMatrix &hm1) const
Definition: SymMatrix.cc:816
virtual int num_row() const
Definition: Vector.cc:117
void fireArray(const int size, HepVector *array)
RandMultiGauss(HepRandomEngine &anEngine, const HepVector &mu, const HepSymMatrix &S)
#define exit(x)
HepMatrix diagonalize(HepSymMatrix *s)
Definition: excDblThrow.cc:17