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

testVectorDists.cc
Go to the documentation of this file.
1 // $Id: testVectorDists.cc,v 1.2 2003/08/13 20:00:13 garren Exp $
2 // -*- C++ -*-
3 // ----------------------------------------------------------------------
4 
5 // ----------------------------------------------------------------------
6 //
7 // testVectorDists -- tests of the correctness of vector random distributions
8 //
9 // Currently tested:
10 // RandMultiGauss
11 //
12 // ----------------------------------------------------------------------
13 
15 #include "CLHEP/Random/Randomize.h"
17 #include "CLHEP/Matrix/SymMatrix.h"
18 #include "CLHEP/Matrix/Matrix.h"
19 #include "CLHEP/Matrix/Vector.h"
20 #include <iostream>
21 
22 using namespace std;
23 using namespace CLHEP;
24 
25 static const int MultiGaussBAD = 1 << 0;
26 
27 
28 static HepMatrix modifiedOutput(const HepMatrix& D) {
29  HepMatrix DD (D);
30  int n = DD.num_row();
31  int m = DD.num_col();
32  int i;
33  int j;
34  for ( i = 1; i <= n; ++i ) {
35  for ( j = 1; j <= m; ++j ) {
36  if ( DD(i,j)*DD(i,j) < 1.0e-24 * DD(i,i) * DD(j,j) ) DD (i,j) = 0;
37  }
38  }
39  return DD;
40 }
41 
42 
43 // --------------
44 // RandMultiGauss
45 // --------------
46 
48 
49  cout << "\n--------------------------------------------\n";
50  cout << "Test of RandMultiGauss distribution \n\n";
51 
52  long seed;
53  cout << "Please enter an integer seed: ";
54  cin >> seed;
55 
56  int nvectors;
57  cout << "How many vectors should we generate: ";
58  cin >> nvectors;
59  double rootn = sqrt((double)nvectors);
60 
61  int nMu;
62  int nS;
63  cout << "Enter the dimensions of mu, then S (normally the same): ";
64  cin >> nMu >> nS;
65  if ( nMu != nS ) {
66  cout << "Usually mu and S will be of equal dimensions.\n";
67  cout << "You may be testing the behavior when that is not the case.\n";
68  cout << "Please verify by re-entering the correct dimensions: ";
69  cin >> nMu >> nS;
70  }
71  int dim = (nMu >= nS) ? nMu : nS;
72 
73  HepVector mu(nMu);
74  HepSymMatrix S(nS);
75 
76  cout << "Enter mu, one component at a time: \n";
77  int imu;
78  double muElement;
79  for (imu = 1; imu <= nMu; imu++) {
80  cout << imu << ": ";
81  cin >> muElement;
82  mu(imu) = muElement;
83  }
84 
85  cout << "Enter S, one white-space-separated row at a time. \n";
86  cout << "The length needed for each row is given in {braces}.\n";
87  cout <<
88  "The diagonal elements of S will be the first numbers on each line:\n";
89  int row, col;
90  double sij;
91  for (row = 1; row <= nS; row++) {
92  cout << row << " {" << nS - row + 1 << " numbers}: ";
93  for (col = row; col <= nS; col++) {
94  cin >> sij;
95  S(row, col) = sij;
96  }
97  }
98 
99  cout << "mu is: \n";
100  cout << mu;
101  cout << endl;
102 
103  cout << "S is: \n";
104  cout << S << endl;
105 
106  HepSymMatrix tempS ( S ); // Since diagonalize does not take a const s
107  // we have to copy S.
108  HepMatrix U = diagonalize ( &tempS );
109  HepSymMatrix D = S.similarityT(U);
110  cout << "D = Diagonalized S is " << modifiedOutput(D) << endl;
111  bool pd = true;
112  for ( int npd = 1; npd <= dim; npd++) {
113  if ( D(npd,npd) < 0 ) {
114  pd = false;
115  }
116  }
117  if (!pd) {
118  cout << "S is not positive definite.\n" <<
119  "The negative elements of D will have been raised to zero.\n" <<
120  "The second moment matrix at the end will not match S.\n";
121  }
122 
123  cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
124  TripleRand eng (seed);
125  RandMultiGauss dist (eng, mu, S);
126 
127  HepVector x(dim);
128 
129  cout << "\n Sample fire(): \n";
130 
131  x = dist.fire();
132  cout << x;
133 
134  int ntrials;
135  cout << "Normal operation will try a single group of " << nvectors
136  << " random vectors.\n"
137  << "Enter 1 for a single trial with " << nvectors
138  << " random vectors.\n"
139  << "Alternatively some number of groups of " << nvectors
140  << " vectors can be produced to accumulate deviation statistics.\n"
141  << "Enter " << 5000/(dim*(dim+1))+1
142  << " or some other number of trials to do this: ";
143  cin >> ntrials;
144  if (ntrials < 1) return 0;
145 
146  cout << "\n Testing fire() ... \n";
147 
148 // I'm going to accumulate correlation matrix by equation (28.9) RPP
149 // and compare to the specified matrix S. That is, E(x-<x>)(y-<y>) should
150 // be Sxy.
151 //
152 // For off-diagonal terms, Sxy = <xy> - <x><y>.
153 //
154 // For diagonal terms, Sxx = <x^2> - <x>^2.
155 
156  HepSymMatrix Sumxy(nS);
157  HepSymMatrix Dprime(dim);
158  HepSymMatrix VarD(dim);
159  HepSymMatrix Delta(dim);
160 
161  int ipr = nvectors / 10; if (ipr < 1) ipr = 1;
162  int in1 = 0;
163  int in2 = 0;
164  int in3 = 0;
165  int nentries = 0;
166  float binno;
167  int nbin;
168  int bins[30];
169  int ix, iy;
170 // double root2 = sqrt(2.0);
171  double sumDelta = 0.0;
172  double sumDelta2 = 0.0;
173  int nunder = 0;
174  int nover = 0;
175  double worstDeviation=0;
176 
177  int k;
178  for(k=0; k<30; ++k) {
179  bins[k] = 0;
180  }
181  for(k=0; k<ntrials; ++k ) {
182  HepVector sumx(dim,0);
183  HepMatrix sumxy(dim,dim,0);
184  for ( int ifire = 1; ifire <= nvectors; ifire++) {
185  x = dist.fire();
186  for (ix = 1; ix <= dim; ix++) {
187  sumx(ix) += x(ix);
188  for (iy = 1; iy <= dim; iy++) {
189  sumxy(ix,iy) += x(ix)*x(iy);
190  }
191  }
192  if ( (ifire % ipr) == 0 && k == 0 ) {
193  cout << ifire << ", ";
194  }
195  }
196  if( k == 0 )
197  cout << "Statistics for the first (or only) trial of " << nvectors
198  << " vectors:\n\n";
199 
200  sumx = sumx / nvectors;
201  if( k == 0 ) cout << "\nAverage (should match mu): \n" << sumx << endl;
202  for (ix = 1; ix <= dim; ix++) {
203  for (iy = 1; iy <= dim; iy++) {
204  sumxy(ix,iy) = sumxy(ix,iy) / nvectors - sumx(ix)*sumx(iy);
205  }
206  }
207  if (pd) {
208  if( k == 0 ) cout << "Second Moments (should match S)\n" << sumxy << endl;
209  } else {
210  if( k == 0 ) cout << "Second Moments \n" << sumxy << endl;
211  }
212 
213 // Now transform sumxy with the same matrix that diagonalized S. Call the
214 // result Dprime. There is a bit of fooling around here because sumxy is a
215 // general matrix and similarityT() acts on a symmetric matrix.
216 
217  Sumxy.assign(sumxy);
218  Dprime = Sumxy.similarityT(U);
219 
220  if( k == 0 ) cout << "\nDprime = Second moment matrix transformed by the same matrix that diagonalized S\n" << Dprime << endl;
221 
222  for( ix=1; ix<=dim; ++ix ) {
223  for( iy=ix; iy<=dim; ++iy ) {
224  if( ix == iy ) {
225  VarD(ix,iy) = 2.0*Dprime(ix,iy)*Dprime(ix,iy)/rootn;
226  } else {
227  VarD(ix,iy) = Dprime(ix,ix)*Dprime(iy,iy)/rootn;
228  }
229  }
230  }
231  if( k == 0 ) cout << "\nThe expected variance for Dprime elements is \n"
232  << VarD << endl;
233 
234  for( ix=1; ix<=dim; ++ix ) {
235  for( iy=ix; iy<=dim; ++iy ) {
236  Delta(ix,iy) = sqrt(rootn)*(D(ix,iy)-Dprime(ix,iy))/sqrt(VarD(ix,iy));
237  if (k==0) {
238  if (abs(Delta(ix,iy)) > abs(worstDeviation)) {
239  worstDeviation = Delta(ix,iy);
240  }
241  }
242  }
243  }
244 
245  if( k == 0 ) {
246  cout
247  << "\nDifferences between each element and its normative value,\n"
248  << "scaled by the expected deviation (sqrt(variance)) are: \n"
249  << Delta << endl;
250  }
251 
252  if( k == 0 ) {
253  cout <<
254  "About 5% of the above values should have absolute value more than 2.\n"
255  << "Deviations of more than 4 sigma would probably indicate a problem.\n";
256  }
257 
258 // Do a little counting
259 
260  for( ix=1; ix<=dim; ++ix ) {
261  for( iy=ix; iy<=dim; ++iy ) {
262  if( Delta(ix,iy) >= -1.0 && Delta(ix,iy) <= 1.0 ) in1++;
263  if( Delta(ix,iy) >= -2.0 && Delta(ix,iy) <= 2.0 ) in2++;
264  if( Delta(ix,iy) >= -3.0 && Delta(ix,iy) <= 3.0 ) in3++;
265  sumDelta += Delta(ix,iy);
266  sumDelta2 += Delta(ix,iy)*Delta(ix,iy);
267  binno = 5.0*(Delta(ix,iy)+3.0);
268  if( binno < 0.0 ) ++nunder;
269  if( binno > 30.0 ) ++nover;
270  if( binno >= 0.0 && binno < 30.0 ) {
271  nbin = (int)binno;
272  ++nentries;
273  ++bins[nbin];
274  }
275  }
276  }
277  }
278 
279  if (ntrials == 1) {
280  cout << "The worst deviation of any element of D in this trial was "
281  << worstDeviation << "\n";
282  if (abs(worstDeviation) > 4) {
283  cout << "\nREJECT this distribution: \n"
284  << "This value indicates a PROBLEM!!!!\n\n";
285  return MultiGaussBAD;
286  } else {
287  return 0;
288  }
289  }
290 
291  float ndf = ntrials*dim*(dim+1)/2.0;
292  cout << "\nOut of a total of " << ndf << " entries" << endl;
293  cout << "There are " << in1 << " within 1 sigma or "
294  << 100.0*(float)in1/ndf << "%" << endl;
295  cout << "There are " << in2 << " within 2 sigma or "
296  << 100.0*(float)in2/ndf << "%" << endl;
297  cout << "There are " << in3 << " within 3 sigma or "
298  << 100.0*(float)in3/ndf << "%" << endl;
299  double aveDelta = sumDelta/(double)ndf;
300  double rmsDelta = sumDelta2/(double)ndf - aveDelta*aveDelta;
301  cout << "\nFor dim = " << dim << " Average(Delta) = " << aveDelta << " RMS(Delta) = " << rmsDelta << endl;
302  cout << "\nPoor man's histogram of deviations in 30 bins from -3.0 to 3.0" << endl;
303  cout << "This should be a standard unit Gaussian.\n" << endl;
304  for(k=0; k<30; ++k ) {
305  cout << setw(3) << k+1 << " " << setw(4) << bins[k] << endl;
306  }
307  cout << "\nTotal number of entries in this histogram is " << nentries << endl;
308  cout << "\twith " << nunder << " underflow(s) and " << nover << " overflow(s)" << endl;
309 
310  int status;
311 
312  cout << "The mean squared delta should be between .85 and 1.15; it is "
313  << rmsDelta << "\n";
314 
315  if( abs(rmsDelta-1.0) <= 0.15 ) {
316  status = false;
317  } else {
318  cout << "REJECT this distribution based on improper spread of delta!\n";
319  status = MultiGaussBAD;
320  }
321  if (abs(worstDeviation)>4) {
322  cout << "REJECT this distribution: Bad correlation in first trial!\n";
323  status = MultiGaussBAD;
324  }
325 
326  return status;
327 
328 } // testRandMultiGauss()
329 
330 int main() {
331 
332  int mask = 0;
333 
334  mask |= testRandMultiGauss();
335 
336  return mask;
337 
338 }
339 
340 
double
#define double(obj)
Definition: excDblThrow.cc:32
CLHEP::TripleRand
Definition: Matrix/CLHEP/Random/TripleRand.h:52
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
CLHEP::HepSymMatrix::assign
void assign(const HepMatrix &hm2)
Definition: SymMatrix.cc:718
D
Definition: excDblThrow.cc:17
CLHEP::detail::n
n
Definition: Ranlux64Engine.cc:85
CLHEP::HepMatrix
Definition: Matrix/CLHEP/Matrix/Matrix.h:209
testRandMultiGauss
int testRandMultiGauss()
Definition: testVectorDists.cc:47
CLHEP::diagonalize
HepMatrix diagonalize(HepSymMatrix *s)
Definition: MatrixLinear.cc:315
CLHEP
Definition: ClhepVersion.h:13
CLHEP::RandMultiGauss
Definition: CLHEP/RandomObjects/RandMultiGauss.h:45
main
int main()
Definition: testVectorDists.cc:330
CLHEP::HepSymMatrix
Definition: Matrix/CLHEP/Matrix/SymMatrix.h:89
j
long j
Definition: JamesRandomSeeding.txt:28
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::HepVector
Definition: Matrix/CLHEP/Matrix/Vector.h:39
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
k
long k
Definition: JamesRandomSeeding.txt:29
n_constructors::m
int m
Definition: testSharedPtr.cc:370
defs.h