JAPAn
Just Another Parity Analyzer
Loading...
Searching...
No Matches
LinReg_Bevington_Pebay.cc
Go to the documentation of this file.
1/*!
2 * \file LinReg_Bevington_Pebay.cc
3 * \brief Implementation of linear regression utility using Bevington and Pebay algorithms
4 * \author Jan Balewski, MIT
5 * \date 2010
6 */
7
8#include <assert.h>
9#include <math.h>
10
11#include "TString.h"
12
13#include "LinReg_Bevington_Pebay.h"
14
15#include "QwLog.h"
16
17//=================================================
18//=================================================
24
25//=================================================
26//=================================================
27LinRegBevPeb::LinRegBevPeb(const LinRegBevPeb& source)
28: nP(source.nP),nY(source.nY),
29 fErrorFlag(-1),
30 fGoodEventNumber(0)
31{
32 QwMessage << fGoodEventNumber << QwLog::endl;
33}
34
35//=================================================
36//=================================================
37void LinRegBevPeb::init()
38{
39 mMP.ResizeTo(nP);
40 mMY.ResizeTo(nY);
41 mMYp.ResizeTo(nY);
42
43 mVPP.ResizeTo(nP,nP);
44 mVPY.ResizeTo(nP,nY);
45 mVYP.ResizeTo(nY,nP);
46 mVYY.ResizeTo(nY,nY);
47 mVYYp.ResizeTo(nY,nY);
48 mVP.ResizeTo(nP);
49 mVY.ResizeTo(nY);
50 mVYp.ResizeTo(nY);
51
52 mSPP.ResizeTo(mVPP);
53 mSPY.ResizeTo(mVPY);
54 mSYP.ResizeTo(mVYP);
55 mSYY.ResizeTo(mVYY);
56 mSYYp.ResizeTo(mVYYp);
57
58 Axy.ResizeTo(nP,nY);
59 Ayx.ResizeTo(nY,nP);
60 dAxy.ResizeTo(Axy);
61 dAyx.ResizeTo(Ayx);
62
63 mSP.ResizeTo(nP);
64 mSY.ResizeTo(nY);
65 mSYp.ResizeTo(nY);
66
67 mRPP.ResizeTo(mVPP);
68 mRPY.ResizeTo(mVPY);
69 mRYP.ResizeTo(mVYP);
70 mRYY.ResizeTo(mVYY);
71 mRYYp.ResizeTo(mVYYp);
72
73 fGoodEventNumber = 0;
74}
75
76void LinRegBevPeb::clear()
77{
78 mMP.Zero();
79 mMY.Zero();
80 mMYp.Zero();
81
82 mVPP.Zero();
83 mVPY.Zero();
84 mVYP.Zero();
85 mVYY.Zero();
86 mVYYp.Zero();
87 mVP.Zero();
88 mVY.Zero();
89 mVYp.Zero();
90
91 mSPP.Zero();
92 mSPY.Zero();
93 mSYP.Zero();
94 mSYY.Zero();
95 mSYYp.Zero();
96
97 Axy.Zero();
98 Ayx.Zero();
99 dAxy.Zero();
100 dAyx.Zero();
101
102 mSP.Zero();
103 mSY.Zero();
104 mSYp.Zero();
105
106 mRPP.Zero();
107 mRPY.Zero();
108 mRYP.Zero();
109 mRYY.Zero();
110 mRYYp.Zero();
111
112 fErrorFlag = -1;
113 fGoodEventNumber = 0;
114}
115
116//=================================================
117//=================================================
118void LinRegBevPeb::print()
119{
120 QwMessage << "LinReg dims: nP=" << nP << " nY=" << nY << QwLog::endl;
121
122 QwMessage << "MP:"; mMP.Print();
123 QwMessage << "MY:"; mMY.Print();
124 QwMessage << "VPP:"; mVPP.Print();
125 QwMessage << "VPY:"; mVPY.Print();
126 QwMessage << "VYY:"; mVYY.Print();
127 QwMessage << "VYYprime:"; mVYYp.Print();
128}
129
130
131//==========================================================
132//==========================================================
133LinRegBevPeb& LinRegBevPeb::operator+=(const std::pair<TVectorD,TVectorD>& rhs)
134{
135 // Get independent and dependent components
136 const TVectorD& P = rhs.first;
137 const TVectorD& Y = rhs.second;
138
139 // Update number of events
140 fGoodEventNumber++;
141
142 if (fGoodEventNumber <= 1) {
143 // First event, set covariances to zero and means to first value
144 mVPP.Zero();
145 mVPY.Zero();
146 mVYY.Zero();
147 mMP = P;
148 mMY = Y;
149 } else {
150 // Deviations from mean
151 TVectorD delta_y(Y - mMY);
152 TVectorD delta_p(P - mMP);
153
154 // Update covariances
155 Double_t alpha = (fGoodEventNumber - 1.0) / fGoodEventNumber;
156 mVPP.Rank1Update(delta_p, alpha);
157 mVPY.Rank1Update(delta_p, delta_y, alpha);
158 mVYY.Rank1Update(delta_y, alpha);
159
160 // Update means
161 Double_t beta = 1.0 / fGoodEventNumber;
162 mMP += delta_p * beta;
163 mMY += delta_y * beta;
164 }
165
166 return *this;
167}
168
169
170//==========================================================
171//==========================================================
172LinRegBevPeb& LinRegBevPeb::operator+=(const LinRegBevPeb& rhs)
173{
174 // If set X = A + B, then
175 // Cov[X] = Cov[A] + Cov[B]
176 // + (E[x_A] - E[x_B]) * (E[y_A] - E[y_B]) * n_A * n_B / n_X
177 // Ref: E. Schubert, M. Gertz (9 July 2018).
178 // "Numerically stable parallel computation of (co-)variance".
179 // SSDBM '18 Proceedings of the 30th International Conference
180 // on Scientific and Statistical Database Management.
181 // https://doi.org/10.1145/3221269.3223036
182
183 if(fGoodEventNumber + rhs.fGoodEventNumber == 0)
184 return *this;
185
186 // Deviations from mean
187 TVectorD delta_y(mMY - rhs.mMY);
188 TVectorD delta_p(mMP - rhs.mMP);
189
190 // Update covariances
191 Double_t alpha = fGoodEventNumber * rhs.fGoodEventNumber
192 / (fGoodEventNumber + rhs.fGoodEventNumber);
193 mVYY += rhs.mVYY;
194 mVYY.Rank1Update(delta_y, alpha);
195 mVPY += rhs.mVPY;
196 mVPY.Rank1Update(delta_p, delta_y, alpha);
197 mVPP += rhs.mVPP;
198 mVPP.Rank1Update(delta_p, alpha);
199
200 // Update means
201 Double_t beta = rhs.fGoodEventNumber / (fGoodEventNumber + rhs.fGoodEventNumber);
202 mMY += delta_y * beta;
203 mMP += delta_p * beta;
204
205 fGoodEventNumber += rhs.fGoodEventNumber;
206
207 return *this;
208}
209
210
211//==========================================================
212//==========================================================
213Int_t LinRegBevPeb::getMeanP(const int i, Double_t &mean) const
214{
215 mean=-1e50;
216 if(i<0 || i >= nP ) return -1;
217 if( fGoodEventNumber<1) return -3;
218 mean = mMP(i); return 0;
219}
220
221
222//==========================================================
223//==========================================================
224Int_t LinRegBevPeb::getMeanY(const int i, Double_t &mean) const
225{
226 mean=-1e50;
227 if(i<0 || i >= nY ) return -1;
228 if( fGoodEventNumber<1) return -3;
229 mean = mMY(i); return 0;
230}
231
232
233//==========================================================
234//==========================================================
235Int_t LinRegBevPeb::getMeanYprime(const int i, Double_t &mean) const
236{
237 mean=-1e50;
238 if(i<0 || i >= nY ) return -1;
239 if( fGoodEventNumber<1) return -3;
240 mean = mMYp(i); return 0;
241}
242
243
244//==========================================================
245//==========================================================
246Int_t LinRegBevPeb::getSigmaP(const int i, Double_t &sigma) const
247{
248 sigma=-1e50;
249 if(i<0 || i >= nP ) return -1;
250 if( fGoodEventNumber<2) return -3;
251 sigma=sqrt(mVPP(i,i)/(fGoodEventNumber-1.));
252 return 0;
253}
254
255
256//==========================================================
257//==========================================================
258Int_t LinRegBevPeb::getSigmaY(const int i, Double_t &sigma) const
259{
260 sigma=-1e50;
261 if(i<0 || i >= nY ) return -1;
262 if( fGoodEventNumber<2) return -3;
263 sigma=sqrt(mVYY(i,i)/(fGoodEventNumber-1.));
264 return 0;
265}
266
267//==========================================================
268//==========================================================
269Int_t LinRegBevPeb::getSigmaYprime(const int i, Double_t &sigma) const
270{
271 sigma=-1e50;
272 if(i<0 || i >= nY ) return -1;
273 if( fGoodEventNumber<2) return -3;
274 sigma=sqrt(mVYYp(i,i)/(fGoodEventNumber-1.));
275 return 0;
276}
277
278//==========================================================
279//==========================================================
280Int_t LinRegBevPeb::getCovarianceP( int i, int j, Double_t &covar) const
281{
282 covar=-1e50;
283 if( i>j) { int k=i; i=j; j=k; }//swap i & j
284 //... now we need only upper right triangle
285 if(i<0 || i >= nP ) return -11;
286 if( fGoodEventNumber<2) return -14;
287 covar=mVPP(i,j)/(fGoodEventNumber-1.);
288 return 0;
289}
290
291//==========================================================
292//==========================================================
293Int_t LinRegBevPeb::getCovariancePY( int ip, int iy, Double_t &covar) const
294{
295 covar=-1e50;
296 //... now we need only upper right triangle
297 if(ip<0 || ip >= nP ) return -11;
298 if(iy<0 || iy >= nY ) return -12;
299 if( fGoodEventNumber<2) return -14;
300 covar=mVPY(ip,iy)/(fGoodEventNumber-1.);
301 return 0;
302}
303
304//==========================================================
305//==========================================================
306Int_t LinRegBevPeb::getCovarianceY( int i, int j, Double_t &covar) const
307{
308 covar=-1e50;
309 if( i>j) { int k=i; i=j; j=k; }//swap i & j
310 //... now we need only upper right triangle
311 if(i<0 || i >= nY ) return -11;
312 if( fGoodEventNumber<2) return -14;
313 covar=mVYY(i,j)/(fGoodEventNumber-1.);
314 return 0;
315}
316
317//==========================================================
318//==========================================================
319void LinRegBevPeb::printSummaryP() const
320{
321 QwMessage << Form("\nLinRegBevPeb::printSummaryP seen good eve=%lld",fGoodEventNumber)<<QwLog::endl;
322
323 size_t dim=nP;
324 if(fGoodEventNumber>2) { // print full matrix
325 QwMessage << Form("\nname: ");
326 for (size_t i = 1; i <dim; i++) {
327 QwMessage << Form("P%d%11s",(int)i," ");
328 }
329 QwMessage << Form("\n mean sig(distrib) nSig(mean) correlation-matrix ....\n");
330 for (size_t i = 0; i <dim; i++) {
331 double meanI,sigI;
332 if (getMeanP(i,meanI) < 0) QwWarning << "LRB::getMeanP failed" << QwLog::endl;
333 if (getSigmaP(i,sigI) < 0) QwWarning << "LRB::getSigmaP failed" << QwLog::endl;
334 double nSig=-1;
335 double err=sigI/sqrt(fGoodEventNumber);
336 if(sigI>0.) nSig=meanI/err;
337
338 QwMessage << Form("P%d: %+12.4g %12.3g %.1f ",(int)i,meanI,sigI,nSig);
339 for (size_t j = 1; j <dim; j++) {
340 if( j<=i) { QwMessage << Form(" %12s","._._._."); continue;}
341 double sigJ,cov;
342 if (getSigmaP(j,sigJ) < 0) QwWarning << "LRB::getSigmaP failed" << QwLog::endl;
343 if (getCovarianceP(i,j,cov) < 0) QwWarning << "LRB::getCovarianceP failed" << QwLog::endl;
344 double corel=cov / sigI / sigJ;
345
346 QwMessage << Form(" %12.3g",corel);
347 }
348 QwMessage << Form("\n");
349 }
350 }
351}
352
353
354//==========================================================
355//==========================================================
356void LinRegBevPeb::printSummaryY() const
357{
358 QwMessage << Form("\nLinRegBevPeb::printSummaryY seen good eve=%lld (CSV-format)",fGoodEventNumber)<<QwLog::endl;
359 QwMessage << Form(" j, mean, sig(mean), nSig(mean), sig(distribution) \n");
360
361 for (int i = 0; i <nY; i++) {
362 double meanI,sigI;
363 if (getMeanY(i,meanI) < 0) QwWarning << "LRB::getMeanY failed" << QwLog::endl;
364 if (getSigmaY(i,sigI) < 0) QwWarning << "LRB::getSigmaY failed" << QwLog::endl;
365 double err = sigI / sqrt(fGoodEventNumber);
366 double nSigErr = meanI / err;
367 QwMessage << Form("Y%02d, %+11.4g, %12.4g, %8.1f, %12.4g "" ",i,meanI,err,nSigErr,sigI)<<QwLog::endl;
368
369 }
370}
371
372
373
374//==========================================================
375//==========================================================
376void LinRegBevPeb::printSummaryAlphas() const
377{
378 QwMessage << Form("\nLinRegBevPeb::printSummaryAlphas seen good eve=%lld",fGoodEventNumber)<<QwLog::endl;
379 QwMessage << Form("\n j slope sigma mean/sigma\n");
380 for (int iy = 0; iy <nY; iy++) {
381 QwMessage << Form("dv=Y%d: ",iy)<<QwLog::endl;
382 for (int j = 0; j < nP; j++) {
383 double val=Axy(j,iy);
384 double err=dAxy(j,iy);
385 double nSig=val/err;
386 char x=' ';
387 if(fabs(nSig)>3.) x='*';
388 QwMessage << Form(" slope_%d = %11.3g +/-%11.3g (nSig=%.2f) %c\n",j,val, err,nSig,x);
389 }
390 }
391}
392
393
394//==========================================================
395//==========================================================
396void LinRegBevPeb::printSummaryYP() const
397{
398 QwMessage << Form("\nLinRegBevPeb::printSummaryYP seen good eve=%lld",fGoodEventNumber)<<QwLog::endl;
399
400 if(fGoodEventNumber<2) { QwMessage<<" too few events, skip"<<QwLog::endl; return;}
401 QwMessage << Form("\n name: ");
402 for (int i = 0; i <nP; i++) {
403 QwMessage << Form(" %10sP%d "," ",i);
404 }
405 QwMessage << Form("\n j meanY sigY correlation with Ps ....\n");
406 for (int iy = 0; iy <nY; iy++) {
407 double meanI,sigI;
408 if (getMeanY(iy,meanI) < 0) QwWarning << "LRB::getMeanY failed" << QwLog::endl;
409 if (getSigmaY(iy,sigI) < 0) QwWarning << "LRB::getSigmaY failed" << QwLog::endl;
410
411 QwMessage << Form(" %3d %6sY%d: %+12.4g %12.4g ",iy," ",iy,meanI,sigI);
412 for (int ip = 0; ip <nP; ip++) {
413 double sigJ,cov;
414 if (getSigmaP(ip,sigJ) < 0) QwWarning << "LRB::getSigmaP failed" << QwLog::endl;
415 if (getCovariancePY(ip,iy,cov) < 0) QwWarning << "LRB::getCovariancePY failed" << QwLog::endl;
416 double corel = cov / sigI / sigJ;
417 QwMessage << Form(" %12.3g",corel);
418 }
419 QwMessage << Form("\n");
420 }
421}
422
423
424//==========================================================
425//==========================================================
426void LinRegBevPeb::printSummaryMeansWithUnc() const
427{
428 QwMessage << "Uncorrected Y values:" << QwLog::endl;
429 QwMessage << " mean sig" << QwLog::endl;
430 for (int i = 0; i < nY; i++){
431 QwMessage << "Y" << i << ": " << mMY(i) << " +- " << mSY(i) << QwLog::endl;
432 }
433 QwMessage << QwLog::endl;
434}
435
436
437//==========================================================
438//==========================================================
439void LinRegBevPeb::printSummaryMeansWithUncCorrected() const
440{
441 QwMessage << "Corrected Y values:" << QwLog::endl;
442 QwMessage << " mean sig" << QwLog::endl;
443 for (int i = 0; i < nY; i++){
444 QwMessage << "Y" << i << ": " << mMYp(i) << " +- " << mSYp(i) << QwLog::endl;
445 }
446 QwMessage << QwLog::endl;
447}
448
449
450//==========================================================
451//==========================================================
452void LinRegBevPeb::solve()
453{
454 // off-diagonal raw covariance
455 mVYP.Transpose(mVPY);
456
457 // diagonal variances
458 mVP = TMatrixDDiag(mVPP); mVP.Sqrt();
459 mVY = TMatrixDDiag(mVYY); mVY.Sqrt();
460
461 // correlation matrices
462 mRPP = mVPP; mRPP.NormByColumn(mVP); mRPP.NormByRow(mVP);
463 mRYY = mVYY; mRYY.NormByColumn(mVY); mRYY.NormByRow(mVY);
464 mRPY = mVPY; mRPY.NormByColumn(mVP); mRPY.NormByRow(mVY);
465
466 /// normalized covariances
467 mSYY = mVYY * (1.0 / (fGoodEventNumber - 1.));
468 mSPP = mVPP * (1.0 / (fGoodEventNumber - 1.));
469 mSPY = mVPY * (1.0 / (fGoodEventNumber - 1.));
470 mSYP.Transpose(mSPY);
471
472 // uncertainties on the means
473 mSP = TMatrixDDiag(mSPP); mSP.Sqrt();
474 mSY = TMatrixDDiag(mSYY); mSY.Sqrt();
475
476 // Warn if correlation matrix determinant close to zero (heuristic)
477 if (mRPP.Determinant() < std::pow(10,-(2*nP))) {
478 QwWarning << "LRB: correlation matrix nearly singular, "
479 << "determinant = " << mRPP.Determinant()
480 << " (set includes highly correlated variable pairs)"
481 << QwLog::endl;
482 if (fGoodEventNumber > 10*nP) {
483 QwMessage << fGoodEventNumber << " events" << QwLog::endl;
484 QwMessage << "Covariance matrix: " << QwLog::endl; mVPP.Print();
485 QwMessage << "Correlation matrix: " << QwLog::endl; mRPP.Print();
486 }
487 QwWarning << "LRB: solving failed (this happens when only few events)."
488 << QwLog::endl;
489 return;
490 }
491 // slopes
492 TMatrixD invRPP(TMatrixD::kInverted, mRPP);
493 Axy = TMatrixD(invRPP, TMatrixD::kMult, mRPY);
494 Axy.NormByColumn(mSP); // divide
495 Axy.NormByRow(mSY, ""); // mult
496 Ayx.Transpose(Axy);
497
498 // new means
499 mMYp = mMY - Ayx * mMP;
500
501 // new raw covariance
502 mVYYp = mVYY + Ayx * mVPP * Axy - (Ayx * mVPY + mVYP * Axy);
503 // new variances
504 mVYp = TMatrixDDiag(mVYYp); mVYp.Sqrt();
505
506 // new normalized covariance
507 mSYYp = mSYY + Ayx * mSPP * Axy - (Ayx * mSPY + mSYP * Axy);
508 // uncertainties on the new means
509 mSYp = TMatrixDDiag(mSYYp); mSYp.Sqrt();
510
511 // new correlation matrix
512 mRYYp = mVYYp; mRYYp.NormByColumn(mVYp); mRYYp.NormByRow(mVYp);
513
514 // slope uncertainties
515 double norm = 1. / (fGoodEventNumber - nP - 1);
516 dAxy.Zero();
517 dAxy.Rank1Update(TMatrixDDiag(invRPP), TMatrixDDiag(mRYYp), norm); // diag mRYYp = row of ones
518 dAxy.Sqrt();
519 dAxy.NormByColumn(mSP); // divide
520 dAxy.NormByRow(mSYp, ""); // mult
521 dAyx.Transpose(dAxy);
522
523 fErrorFlag = 0;
524}
525
QwLog.h
A logfile class, based on an identical class in the Hermes analyzer.
QwWarning
#define QwWarning
Predefined log drain for warnings.
Definition QwLog.h:44
QwMessage
#define QwMessage
Predefined log drain for regular messages.
Definition QwLog.h:49
LinReg_Bevington_Pebay.h
Linear regression utility class based on Bevington and Pebay algorithms.
QwLog::endl
static std::ostream & endl(std::ostream &)
End of the line.
Definition QwLog.cc:297
LinRegBevPeb::getCovarianceP
Int_t getCovarianceP(int i, int j, Double_t &covar) const
Get mean value of a variable, returns error code.
Definition LinReg_Bevington_Pebay.cc:280
LinRegBevPeb::printSummaryYP
void printSummaryYP() const
Definition LinReg_Bevington_Pebay.cc:396
LinRegBevPeb::getCovarianceY
Int_t getCovarianceY(int i, int j, Double_t &covar) const
Definition LinReg_Bevington_Pebay.cc:306
LinRegBevPeb::getSigmaY
Int_t getSigmaY(const int i, Double_t &sigma) const
Definition LinReg_Bevington_Pebay.cc:258
LinRegBevPeb::mVPY
TMatrixD mVPY
unnormalized covariances
Definition LinReg_Bevington_Pebay.h:44
LinRegBevPeb::printSummaryMeansWithUncCorrected
void printSummaryMeansWithUncCorrected() const
Definition LinReg_Bevington_Pebay.cc:439
LinRegBevPeb::printSummaryAlphas
void printSummaryAlphas() const
Definition LinReg_Bevington_Pebay.cc:376
LinRegBevPeb::mRPY
TMatrixD mRPY
correlations
Definition LinReg_Bevington_Pebay.h:39
LinRegBevPeb::getMeanY
Int_t getMeanY(const int i, Double_t &mean) const
Definition LinReg_Bevington_Pebay.cc:224
LinRegBevPeb::getMeanYprime
Int_t getMeanYprime(const int i, Double_t &mean) const
Definition LinReg_Bevington_Pebay.cc:235
LinRegBevPeb::getMeanP
Int_t getMeanP(const int i, Double_t &mean) const
Get mean value of a variable, returns error code.
Definition LinReg_Bevington_Pebay.cc:213
LinRegBevPeb::getCovariancePY
Int_t getCovariancePY(int ip, int iy, Double_t &covar) const
Definition LinReg_Bevington_Pebay.cc:293
LinRegBevPeb::fGoodEventNumber
Long64_t fGoodEventNumber
accumulated so far
Definition LinReg_Bevington_Pebay.h:36
LinRegBevPeb::printSummaryMeansWithUnc
void printSummaryMeansWithUnc() const
Definition LinReg_Bevington_Pebay.cc:426
LinRegBevPeb::getSigmaYprime
Int_t getSigmaYprime(const int i, Double_t &sigma) const
Definition LinReg_Bevington_Pebay.cc:269
LinRegBevPeb::mSPY
TMatrixD mSPY
normalized covariances
Definition LinReg_Bevington_Pebay.h:52
LinRegBevPeb::fErrorFlag
Int_t fErrorFlag
is information valid
Definition LinReg_Bevington_Pebay.h:35
LinRegBevPeb::getSigmaP
Int_t getSigmaP(const int i, Double_t &sigma) const
Get mean value of a variable, returns error code.
Definition LinReg_Bevington_Pebay.cc:246
LinRegBevPeb::mVP
TVectorD mVP
variances
Definition LinReg_Bevington_Pebay.h:48
LinRegBevPeb::operator+=
LinRegBevPeb & operator+=(const std::pair< TVectorD, TVectorD > &rhs)
Definition LinReg_Bevington_Pebay.cc:133
LinRegBevPeb::mMP
TVectorD mMP
mean values
Definition LinReg_Bevington_Pebay.h:60