47 {
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
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
107
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";
126
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
149
150
151
152
153
154
155
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
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 ) {
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
214
215
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
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}
HepMatrix diagonalize(HepSymMatrix *s)
1.9.6