LinRegBevPeb Class Reference

Online linear regression with incremental covariance updates. More...

#include <LinReg_Bevington_Pebay.h>

Public Member Functions
	LinRegBevPeb ()
	LinRegBevPeb (const LinRegBevPeb &source)
virtual	~LinRegBevPeb ()
void	solve ()
bool	failed ()
void	printSummaryP () const
void	printSummaryY () const
void	printSummaryYP () const
void	printSummaryAlphas () const
void	printSummaryMeansWithUnc () const
void	printSummaryMeansWithUncCorrected () const
void	print ()
void	init ()
void	clear ()
void	setDims (int a, int b)
Int_t	getMeanP (const int i, Double_t &mean) const
	Get mean value of a variable, returns error code.
Int_t	getMeanY (const int i, Double_t &mean) const
Int_t	getMeanYprime (const int i, Double_t &mean) const
Int_t	getSigmaP (const int i, Double_t &sigma) const
	Get mean value of a variable, returns error code.
Int_t	getSigmaY (const int i, Double_t &sigma) const
Int_t	getSigmaYprime (const int i, Double_t &sigma) const
Int_t	getCovarianceP (int i, int j, Double_t &covar) const
	Get mean value of a variable, returns error code.
Int_t	getCovariancePY (int ip, int iy, Double_t &covar) const
Int_t	getCovarianceY (int i, int j, Double_t &covar) const
double	getUsedEve () const
LinRegBevPeb &	operator+= (const std::pair< TVectorD, TVectorD > &rhs)
LinRegBevPeb &	operator+= (const LinRegBevPeb &rhs)

Private Attributes
int	nP
int	nY
Int_t	fErrorFlag
	is information valid
Long64_t	fGoodEventNumber
	accumulated so far
TMatrixD	mRPY
	correlations
TMatrixD	mRYP
TMatrixD	mRPP
TMatrixD	mRYY
TMatrixD	mRYYp
TMatrixD	mVPY
	unnormalized covariances
TMatrixD	mVYP
TMatrixD	mVPP
TMatrixD	mVYY
TMatrixD	mVYYp
TVectorD	mVP
	variances
TVectorD	mVY
TVectorD	mVYp
TMatrixD	mSPY
	normalized covariances
TMatrixD	mSYP
TMatrixD	mSPP
TMatrixD	mSYY
TMatrixD	mSYYp
TVectorD	mSP
	sigmas
TVectorD	mSY
TVectorD	mSYp
TVectorD	mMP
	mean values
TVectorD	mMY
TVectorD	mMYp
TMatrixD	Axy
	slopes
TMatrixD	Ayx
TMatrixD	dAxy
TMatrixD	dAyx

Friends
class	QwCorrelator
	Friend class with correlator for ROOT tree output.
LinRegBevPeb	operator+ (LinRegBevPeb lhs, const LinRegBevPeb &rhs)
std::ostream &	operator<< (std::ostream &stream, const LinRegBevPeb &h)
	Output stream operator.

Detailed Description

Online linear regression with incremental covariance updates.

Definition at line 30 of file LinReg_Bevington_Pebay.h.

Constructor & Destructor Documentation

LinRegBevPeb() [1/2]

LinRegBevPeb::LinRegBevPeb

(

)

Definition at line 19 of file LinReg_Bevington_Pebay.cc.

: nP(0),nY(0),
  fErrorFlag(-1),
  fGoodEventNumber(0)
{ }

References fErrorFlag, fGoodEventNumber, nP, and nY.

Referenced by LinRegBevPeb(), operator+, operator+=(), operator+=(), and operator<<.

Here is the caller graph for this function:

LinRegBevPeb() [2/2]

LinRegBevPeb::LinRegBevPeb

(

const LinRegBevPeb &

source

)

Definition at line 27 of file LinReg_Bevington_Pebay.cc.

: nP(source.nP),nY(source.nY),
  fErrorFlag(-1),
  fGoodEventNumber(0)
{
  QwMessage << fGoodEventNumber << QwLog::endl;
}

References QwLog::endl(), fErrorFlag, fGoodEventNumber, LinRegBevPeb(), nP, nY, and QwMessage.

Here is the call graph for this function:

~LinRegBevPeb()

virtual LinRegBevPeb::~LinRegBevPeb ( )

inlinevirtual

Definition at line 70 of file LinReg_Bevington_Pebay.h.

{ };

Member Function Documentation

clear()

void LinRegBevPeb::clear

(

)

Definition at line 76 of file LinReg_Bevington_Pebay.cc.

{
  mMP.Zero();
  mMY.Zero();
  mMYp.Zero();
 
  mVPP.Zero();
  mVPY.Zero();
  mVYP.Zero();
  mVYY.Zero();
  mVYYp.Zero();
  mVP.Zero();
  mVY.Zero();
  mVYp.Zero();
 
  mSPP.Zero();
  mSPY.Zero();
  mSYP.Zero();
  mSYY.Zero();
  mSYYp.Zero();
 
  Axy.Zero();
  Ayx.Zero();
  dAxy.Zero();
  dAyx.Zero();
 
  mSP.Zero();
  mSY.Zero();
  mSYp.Zero();
 
  mRPP.Zero();
  mRPY.Zero();
  mRYP.Zero();
  mRYY.Zero();
  mRYYp.Zero();
 
  fErrorFlag = -1;
  fGoodEventNumber = 0;
}

References Axy, Ayx, dAxy, dAyx, fErrorFlag, fGoodEventNumber, mMP, mMY, mMYp, mRPP, mRPY, mRYP, mRYY, mRYYp, mSP, mSPP, mSPY, mSY, mSYP, mSYp, mSYY, mSYYp, mVP, mVPP, mVPY, mVY, mVYP, mVYp, mVYY, and mVYYp.

failed()

bool LinRegBevPeb::failed ( )

inline

Definition at line 73 of file LinReg_Bevington_Pebay.h.

{ return fGoodEventNumber < nP + 1; }

References fGoodEventNumber, and nP.

getCovarianceP()

Int_t LinRegBevPeb::getCovarianceP	(	int	i,
		int	j,
		Double_t &	covar ) const

Get mean value of a variable, returns error code.

Definition at line 280 of file LinReg_Bevington_Pebay.cc.

{
    covar=-1e50;
    if( i>j) { int k=i; i=j; j=k; }//swap i & j
    //... now we need only upper right triangle
    if(i<0 || i >= nP ) return -11;
    if( fGoodEventNumber<2) return -14;
    covar=mVPP(i,j)/(fGoodEventNumber-1.);
    return 0;
}

References fGoodEventNumber, mVPP, and nP.

Referenced by printSummaryP().

Here is the caller graph for this function:

getCovariancePY()

Int_t LinRegBevPeb::getCovariancePY	(	int	ip,
		int	iy,
		Double_t &	covar ) const

Definition at line 293 of file LinReg_Bevington_Pebay.cc.

{
    covar=-1e50;
    //... now we need only upper right triangle
    if(ip<0 || ip >= nP ) return -11;
    if(iy<0 || iy >= nY ) return -12;
    if( fGoodEventNumber<2) return -14;
    covar=mVPY(ip,iy)/(fGoodEventNumber-1.);
    return 0;
}

References fGoodEventNumber, mVPY, nP, and nY.

Referenced by printSummaryYP().

Here is the caller graph for this function:

getCovarianceY()

Int_t LinRegBevPeb::getCovarianceY	(	int	i,
		int	j,
		Double_t &	covar ) const

Definition at line 306 of file LinReg_Bevington_Pebay.cc.

{
    covar=-1e50;
    if( i>j) { int k=i; i=j; j=k; }//swap i & j
    //... now we need only upper right triangle
    if(i<0 || i >= nY ) return -11;
    if( fGoodEventNumber<2) return -14;
    covar=mVYY(i,j)/(fGoodEventNumber-1.);
    return 0;
}

References fGoodEventNumber, mVYY, and nY.

getMeanP()

Int_t LinRegBevPeb::getMeanP	(	const int	i,
		Double_t &	mean ) const

Get mean value of a variable, returns error code.

Definition at line 213 of file LinReg_Bevington_Pebay.cc.

{
   mean=-1e50;
   if(i<0 || i >= nP ) return -1;
   if( fGoodEventNumber<1) return -3;
   mean = mMP(i);    return 0;
}

References fGoodEventNumber, mMP, and nP.

Referenced by printSummaryP().

Here is the caller graph for this function:

getMeanY()

Int_t LinRegBevPeb::getMeanY	(	const int	i,
		Double_t &	mean ) const

Definition at line 224 of file LinReg_Bevington_Pebay.cc.

{
  mean=-1e50;
  if(i<0 || i >= nY ) return -1;
  if( fGoodEventNumber<1) return -3;
  mean = mMY(i);    return 0;
}

References fGoodEventNumber, mMY, and nY.

Referenced by printSummaryY(), and printSummaryYP().

Here is the caller graph for this function:

getMeanYprime()

Int_t LinRegBevPeb::getMeanYprime	(	const int	i,
		Double_t &	mean ) const

Definition at line 235 of file LinReg_Bevington_Pebay.cc.

{
  mean=-1e50;
  if(i<0 || i >= nY ) return -1;
  if( fGoodEventNumber<1) return -3;
  mean = mMYp(i);    return 0;
}

References fGoodEventNumber, mMYp, and nY.

getSigmaP()

Int_t LinRegBevPeb::getSigmaP	(	const int	i,
		Double_t &	sigma ) const

Get mean value of a variable, returns error code.

Definition at line 246 of file LinReg_Bevington_Pebay.cc.

{
  sigma=-1e50;
  if(i<0 || i >= nP ) return -1;
  if( fGoodEventNumber<2) return -3;
  sigma=sqrt(mVPP(i,i)/(fGoodEventNumber-1.));
  return 0;
}

References fGoodEventNumber, mVPP, and nP.

Referenced by printSummaryP(), and printSummaryYP().

Here is the caller graph for this function:

getSigmaY()

Int_t LinRegBevPeb::getSigmaY	(	const int	i,
		Double_t &	sigma ) const

Definition at line 258 of file LinReg_Bevington_Pebay.cc.

{
  sigma=-1e50;
  if(i<0 || i >= nY ) return -1;
  if( fGoodEventNumber<2) return -3;
  sigma=sqrt(mVYY(i,i)/(fGoodEventNumber-1.));
  return 0;
}

References fGoodEventNumber, mVYY, and nY.

Referenced by printSummaryY(), and printSummaryYP().

Here is the caller graph for this function:

getSigmaYprime()

Int_t LinRegBevPeb::getSigmaYprime	(	const int	i,
		Double_t &	sigma ) const

Definition at line 269 of file LinReg_Bevington_Pebay.cc.

{
  sigma=-1e50;
  if(i<0 || i >= nY ) return -1;
  if( fGoodEventNumber<2) return -3;
  sigma=sqrt(mVYYp(i,i)/(fGoodEventNumber-1.));
  return 0;
}

References fGoodEventNumber, mVYYp, and nY.

getUsedEve()

double LinRegBevPeb::getUsedEve ( ) const

inline

Definition at line 103 of file LinReg_Bevington_Pebay.h.

{ return fGoodEventNumber; };

References fGoodEventNumber.

init()

void LinRegBevPeb::init

(

)

Definition at line 37 of file LinReg_Bevington_Pebay.cc.

{
  mMP.ResizeTo(nP);
  mMY.ResizeTo(nY);
  mMYp.ResizeTo(nY);
 
  mVPP.ResizeTo(nP,nP);
  mVPY.ResizeTo(nP,nY);
  mVYP.ResizeTo(nY,nP);
  mVYY.ResizeTo(nY,nY);
  mVYYp.ResizeTo(nY,nY);
  mVP.ResizeTo(nP);
  mVY.ResizeTo(nY);
  mVYp.ResizeTo(nY);
 
  mSPP.ResizeTo(mVPP);
  mSPY.ResizeTo(mVPY);
  mSYP.ResizeTo(mVYP);
  mSYY.ResizeTo(mVYY);
  mSYYp.ResizeTo(mVYYp);
 
  Axy.ResizeTo(nP,nY);
  Ayx.ResizeTo(nY,nP);
  dAxy.ResizeTo(Axy);
  dAyx.ResizeTo(Ayx);
 
  mSP.ResizeTo(nP);
  mSY.ResizeTo(nY);
  mSYp.ResizeTo(nY);
 
  mRPP.ResizeTo(mVPP);
  mRPY.ResizeTo(mVPY);
  mRYP.ResizeTo(mVYP);
  mRYY.ResizeTo(mVYY);
  mRYYp.ResizeTo(mVYYp);
 
  fGoodEventNumber = 0;
}

References Axy, Ayx, dAxy, dAyx, fGoodEventNumber, mMP, mMY, mMYp, mRPP, mRPY, mRYP, mRYY, mRYYp, mSP, mSPP, mSPY, mSY, mSYP, mSYp, mSYY, mSYYp, mVP, mVPP, mVPY, mVY, mVYP, mVYp, mVYY, mVYYp, nP, and nY.

operator+=() [1/2]

LinRegBevPeb & LinRegBevPeb::operator+=

(

const LinRegBevPeb &

rhs

)

Definition at line 172 of file LinReg_Bevington_Pebay.cc.

{
  // If set X = A + B, then
  //   Cov[X] = Cov[A] + Cov[B]
  //          + (E[x_A] - E[x_B]) * (E[y_A] - E[y_B]) * n_A * n_B / n_X
  // Ref: E. Schubert, M. Gertz (9 July 2018).
  // "Numerically stable parallel computation of (co-)variance".
  // SSDBM '18 Proceedings of the 30th International Conference
  // on Scientific and Statistical Database Management.
  // https://doi.org/10.1145/3221269.3223036
 
  if(fGoodEventNumber + rhs.fGoodEventNumber == 0)
    return *this;
 
  // Deviations from mean
  TVectorD delta_y(mMY - rhs.mMY);
  TVectorD delta_p(mMP - rhs.mMP);
 
  // Update covariances
  Double_t alpha = fGoodEventNumber * rhs.fGoodEventNumber
                / (fGoodEventNumber + rhs.fGoodEventNumber);
  mVYY += rhs.mVYY;
  mVYY.Rank1Update(delta_y, alpha);
  mVPY += rhs.mVPY;
  mVPY.Rank1Update(delta_p, delta_y, alpha);
  mVPP += rhs.mVPP;
  mVPP.Rank1Update(delta_p, alpha);
 
  // Update means
  Double_t beta = rhs.fGoodEventNumber / (fGoodEventNumber + rhs.fGoodEventNumber);
  mMY += delta_y * beta;
  mMP += delta_p * beta;
 
  fGoodEventNumber += rhs.fGoodEventNumber;
 
  return *this;
}

References fGoodEventNumber, LinRegBevPeb(), mMP, mMY, mVPP, mVPY, and mVYY.

Here is the call graph for this function:

operator+=() [2/2]

LinRegBevPeb & LinRegBevPeb::operator+=

(

const std::pair< TVectorD, TVectorD > &

rhs

)

Definition at line 133 of file LinReg_Bevington_Pebay.cc.

{
  // Get independent and dependent components
  const TVectorD& P = rhs.first;
  const TVectorD& Y = rhs.second;
 
  // Update number of events
  fGoodEventNumber++;
 
  if (fGoodEventNumber <= 1) {
    // First event, set covariances to zero and means to first value
    mVPP.Zero();
    mVPY.Zero();
    mVYY.Zero();
    mMP = P;
    mMY = Y;
  } else {
    // Deviations from mean
    TVectorD delta_y(Y - mMY);
    TVectorD delta_p(P - mMP);
 
    // Update covariances
    Double_t alpha = (fGoodEventNumber - 1.0) / fGoodEventNumber;
    mVPP.Rank1Update(delta_p, alpha);
    mVPY.Rank1Update(delta_p, delta_y, alpha);
    mVYY.Rank1Update(delta_y, alpha);
 
    // Update means
    Double_t beta = 1.0 / fGoodEventNumber;
    mMP += delta_p * beta;
    mMY += delta_y * beta;
  }
 
  return *this;
}

References fGoodEventNumber, LinRegBevPeb(), mMP, mMY, mVPP, mVPY, and mVYY.

Here is the call graph for this function:

print()

void LinRegBevPeb::print

(

)

Definition at line 118 of file LinReg_Bevington_Pebay.cc.

{
  QwMessage << "LinReg dims:  nP=" << nP << " nY=" << nY << QwLog::endl;
 
  QwMessage << "MP:"; mMP.Print();
  QwMessage << "MY:"; mMY.Print();
  QwMessage << "VPP:"; mVPP.Print();
  QwMessage << "VPY:"; mVPY.Print();
  QwMessage << "VYY:"; mVYY.Print();
  QwMessage << "VYYprime:"; mVYYp.Print();
}

References QwLog::endl(), mMP, mMY, mVPP, mVPY, mVYY, mVYYp, nP, nY, and QwMessage.

Here is the call graph for this function:

printSummaryAlphas()

void LinRegBevPeb::printSummaryAlphas

(

)

const

Definition at line 376 of file LinReg_Bevington_Pebay.cc.

{
  QwMessage << Form("\nLinRegBevPeb::printSummaryAlphas seen good eve=%lld",fGoodEventNumber)<<QwLog::endl;
  QwMessage << Form("\n  j                slope         sigma     mean/sigma\n");
  for (int iy = 0; iy <nY; iy++) {
    QwMessage << Form("dv=Y%d: ",iy)<<QwLog::endl;
    for (int j = 0; j < nP; j++) {
      double val=Axy(j,iy);
      double err=dAxy(j,iy);
      double nSig=val/err;
      char x=' ';
      if(fabs(nSig)>3.) x='*';
      QwMessage << Form("  slope_%d = %11.3g  +/-%11.3g  (nSig=%.2f) %c\n",j,val, err,nSig,x);
    }
  }
}

References Axy, dAxy, QwLog::endl(), fGoodEventNumber, nP, nY, and QwMessage.

Here is the call graph for this function:

printSummaryMeansWithUnc()

void LinRegBevPeb::printSummaryMeansWithUnc

(

)

const

Definition at line 426 of file LinReg_Bevington_Pebay.cc.

{
  QwMessage << "Uncorrected Y values:" << QwLog::endl;
  QwMessage << "     mean          sig" << QwLog::endl;
  for (int i = 0; i < nY; i++){
    QwMessage << "Y" << i << ":  " << mMY(i) << " +- " << mSY(i) << QwLog::endl;
  }
  QwMessage << QwLog::endl;
}

References QwLog::endl(), mMY, mSY, nY, and QwMessage.

Here is the call graph for this function:

printSummaryMeansWithUncCorrected()

void LinRegBevPeb::printSummaryMeansWithUncCorrected

(

)

const

Definition at line 439 of file LinReg_Bevington_Pebay.cc.

{
  QwMessage << "Corrected Y values:" << QwLog::endl;
  QwMessage << "     mean          sig" << QwLog::endl;
  for (int i = 0; i < nY; i++){
    QwMessage << "Y" << i << ":  " << mMYp(i) << " +- " << mSYp(i) << QwLog::endl;
  }
  QwMessage << QwLog::endl;
}

References QwLog::endl(), mMYp, mSYp, nY, and QwMessage.

Here is the call graph for this function:

printSummaryP()

void LinRegBevPeb::printSummaryP

(

)

const

Definition at line 319 of file LinReg_Bevington_Pebay.cc.

{
  QwMessage << Form("\nLinRegBevPeb::printSummaryP seen good eve=%lld",fGoodEventNumber)<<QwLog::endl;
 
  size_t dim=nP;
  if(fGoodEventNumber>2) { // print full matrix
    QwMessage << Form("\nname:                                               ");
    for (size_t i = 1; i <dim; i++) {
      QwMessage << Form("P%d%11s",(int)i," ");
    }
    QwMessage << Form("\n           mean     sig(distrib)   nSig(mean)       correlation-matrix ....\n");
    for (size_t i = 0; i <dim; i++) {
      double meanI,sigI;
      if (getMeanP(i,meanI) < 0) QwWarning << "LRB::getMeanP failed" << QwLog::endl;
      if (getSigmaP(i,sigI) < 0) QwWarning << "LRB::getSigmaP failed" << QwLog::endl;
      double nSig=-1;
      double err=sigI/sqrt(fGoodEventNumber);
      if(sigI>0.) nSig=meanI/err;
 
      QwMessage << Form("P%d:  %+12.4g  %12.3g      %.1f    ",(int)i,meanI,sigI,nSig);
      for (size_t j = 1; j <dim; j++) {
        if( j<=i) { QwMessage << Form("  %12s","._._._."); continue;}
        double sigJ,cov;
        if (getSigmaP(j,sigJ) < 0) QwWarning << "LRB::getSigmaP failed" << QwLog::endl;
        if (getCovarianceP(i,j,cov) < 0) QwWarning << "LRB::getCovarianceP failed" << QwLog::endl;
        double corel=cov / sigI / sigJ;
        
        QwMessage << Form("  %12.3g",corel);
      }
      QwMessage << Form("\n");
    }
  }
}

References QwLog::endl(), fGoodEventNumber, getCovarianceP(), getMeanP(), getSigmaP(), nP, QwMessage, and QwWarning.

Here is the call graph for this function:

printSummaryY()

void LinRegBevPeb::printSummaryY

(

)

const

Definition at line 356 of file LinReg_Bevington_Pebay.cc.

{
  QwMessage << Form("\nLinRegBevPeb::printSummaryY seen good eve=%lld  (CSV-format)",fGoodEventNumber)<<QwLog::endl;
  QwMessage << Form("  j,       mean,     sig(mean),   nSig(mean),  sig(distribution)    \n");
  
  for (int i = 0; i <nY; i++) {
    double meanI,sigI;
    if (getMeanY(i,meanI) < 0) QwWarning << "LRB::getMeanY failed" << QwLog::endl;
    if (getSigmaY(i,sigI) < 0) QwWarning << "LRB::getSigmaY failed" << QwLog::endl;
    double err = sigI / sqrt(fGoodEventNumber);
    double nSigErr = meanI / err;
    QwMessage << Form("Y%02d, %+11.4g, %12.4g, %8.1f, %12.4g "" ",i,meanI,err,nSigErr,sigI)<<QwLog::endl;
 
  }
}

References QwLog::endl(), fGoodEventNumber, getMeanY(), getSigmaY(), nY, QwMessage, and QwWarning.

Here is the call graph for this function:

printSummaryYP()

void LinRegBevPeb::printSummaryYP

(

)

const

Definition at line 396 of file LinReg_Bevington_Pebay.cc.

{
  QwMessage << Form("\nLinRegBevPeb::printSummaryYP seen good eve=%lld",fGoodEventNumber)<<QwLog::endl;
 
  if(fGoodEventNumber<2) { QwMessage<<"  too few events, skip"<<QwLog::endl; return;}
  QwMessage << Form("\n         name:             ");
  for (int i = 0; i <nP; i++) {
    QwMessage << Form(" %10sP%d "," ",i);
  }
  QwMessage << Form("\n  j                   meanY         sigY      correlation with Ps ....\n");
  for (int iy = 0; iy <nY; iy++) {
    double meanI,sigI;
    if (getMeanY(iy,meanI) < 0) QwWarning << "LRB::getMeanY failed" << QwLog::endl;
    if (getSigmaY(iy,sigI) < 0) QwWarning << "LRB::getSigmaY failed" << QwLog::endl;
    
    QwMessage << Form(" %3d %6sY%d:  %+12.4g  %12.4g ",iy," ",iy,meanI,sigI);
    for (int ip = 0; ip <nP; ip++) {
      double sigJ,cov; 
      if (getSigmaP(ip,sigJ) < 0) QwWarning << "LRB::getSigmaP failed" << QwLog::endl;
      if (getCovariancePY(ip,iy,cov) < 0) QwWarning << "LRB::getCovariancePY failed" << QwLog::endl;
      double corel = cov / sigI / sigJ;
      QwMessage << Form("  %12.3g",corel);
    }
    QwMessage << Form("\n");
  }
}

References QwLog::endl(), fGoodEventNumber, getCovariancePY(), getMeanY(), getSigmaP(), getSigmaY(), nP, nY, QwMessage, and QwWarning.

Here is the call graph for this function:

setDims()

void LinRegBevPeb::setDims ( int a, int b )

inline

Definition at line 86 of file LinReg_Bevington_Pebay.h.

{ nP = a; nY = b;}

References nP, and nY.

solve()

void LinRegBevPeb::solve

(

)

normalized covariances

Definition at line 452 of file LinReg_Bevington_Pebay.cc.

{
  // off-diagonal raw covariance
  mVYP.Transpose(mVPY);
 
  // diagonal variances
  mVP = TMatrixDDiag(mVPP); mVP.Sqrt();
  mVY = TMatrixDDiag(mVYY); mVY.Sqrt();
 
  // correlation matrices
  mRPP = mVPP; mRPP.NormByColumn(mVP); mRPP.NormByRow(mVP);
  mRYY = mVYY; mRYY.NormByColumn(mVY); mRYY.NormByRow(mVY);
  mRPY = mVPY; mRPY.NormByColumn(mVP); mRPY.NormByRow(mVY);
 
  /// normalized covariances
  mSYY = mVYY * (1.0 / (fGoodEventNumber - 1.));
  mSPP = mVPP * (1.0 / (fGoodEventNumber - 1.));
  mSPY = mVPY * (1.0 / (fGoodEventNumber - 1.));
  mSYP.Transpose(mSPY);
 
  // uncertainties on the means
  mSP = TMatrixDDiag(mSPP); mSP.Sqrt();
  mSY = TMatrixDDiag(mSYY); mSY.Sqrt();
 
  // Warn if correlation matrix determinant close to zero (heuristic)
  if (mRPP.Determinant() < std::pow(10,-(2*nP))) {
    QwWarning << "LRB: correlation matrix nearly singular, "
              << "determinant = " << mRPP.Determinant()
              << " (set includes highly correlated variable pairs)"
              << QwLog::endl;
    if (fGoodEventNumber > 10*nP) {
      QwMessage << fGoodEventNumber << " events" << QwLog::endl;
      QwMessage << "Covariance matrix: " << QwLog::endl; mVPP.Print();
      QwMessage << "Correlation matrix: " << QwLog::endl; mRPP.Print();
    }
    QwWarning << "LRB: solving failed (this happens when only few events)."
              << QwLog::endl;
    return;
  }
  // slopes
  TMatrixD invRPP(TMatrixD::kInverted, mRPP);
  Axy = TMatrixD(invRPP, TMatrixD::kMult, mRPY);
  Axy.NormByColumn(mSP); // divide
  Axy.NormByRow(mSY, ""); // mult
  Ayx.Transpose(Axy);
 
  // new means
  mMYp = mMY - Ayx * mMP;
 
  // new raw covariance
  mVYYp = mVYY + Ayx * mVPP * Axy - (Ayx * mVPY + mVYP * Axy);
  // new variances
  mVYp = TMatrixDDiag(mVYYp); mVYp.Sqrt();
 
  // new normalized covariance
  mSYYp = mSYY + Ayx * mSPP * Axy - (Ayx * mSPY + mSYP * Axy);
  // uncertainties on the new means
  mSYp = TMatrixDDiag(mSYYp); mSYp.Sqrt();
 
  // new correlation matrix
  mRYYp = mVYYp; mRYYp.NormByColumn(mVYp); mRYYp.NormByRow(mVYp);
 
  // slope uncertainties
  double norm = 1. / (fGoodEventNumber - nP - 1);
  dAxy.Zero();
  dAxy.Rank1Update(TMatrixDDiag(invRPP), TMatrixDDiag(mRYYp), norm); // diag mRYYp = row of ones
  dAxy.Sqrt();
  dAxy.NormByColumn(mSP); // divide
  dAxy.NormByRow(mSYp, ""); // mult
  dAyx.Transpose(dAxy);
 
  fErrorFlag = 0;
}

References Axy, Ayx, dAxy, dAyx, QwLog::endl(), fErrorFlag, fGoodEventNumber, mMP, mMY, mMYp, mRPP, mRPY, mRYY, mRYYp, mSP, mSPP, mSPY, mSY, mSYP, mSYp, mSYY, mSYYp, mVP, mVPP, mVPY, mVY, mVYP, mVYp, mVYY, mVYYp, nP, QwMessage, and QwWarning.

Here is the call graph for this function:

Friends And Related Symbol Documentation

operator+

LinRegBevPeb operator+ ( LinRegBevPeb lhs, const LinRegBevPeb & rhs )

friend

Definition at line 110 of file LinReg_Bevington_Pebay.h.

  {
    lhs += rhs; // reuse compound assignment
    return lhs; // return the result by value (uses move constructor)
  }

References LinRegBevPeb().

operator<<

std::ostream & operator<< ( std::ostream & stream, const LinRegBevPeb & h )

friend

Output stream operator.

Definition at line 125 of file LinReg_Bevington_Pebay.h.

{
  stream << "LRB: " << h.fGoodEventNumber << " events";
  return stream;
}

References fGoodEventNumber, and LinRegBevPeb().

QwCorrelator

friend class QwCorrelator

friend

Friend class with correlator for ROOT tree output.

Definition at line 121 of file LinReg_Bevington_Pebay.h.

References QwCorrelator.

Referenced by QwCorrelator.

Field Documentation

Axy

TMatrixD LinRegBevPeb::Axy

private

slopes

Definition at line 64 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), printSummaryAlphas(), and solve().

Ayx

TMatrixD LinRegBevPeb::Ayx

private

Definition at line 64 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

dAxy

TMatrixD LinRegBevPeb::dAxy

private

Definition at line 64 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), printSummaryAlphas(), and solve().

dAyx

TMatrixD LinRegBevPeb::dAyx

private

Definition at line 64 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

fErrorFlag

Int_t LinRegBevPeb::fErrorFlag

private

is information valid

Definition at line 35 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), LinRegBevPeb(), LinRegBevPeb(), and solve().

fGoodEventNumber

Long64_t LinRegBevPeb::fGoodEventNumber

private

accumulated so far

Definition at line 36 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), failed(), getCovarianceP(), getCovariancePY(), getCovarianceY(), getMeanP(), getMeanY(), getMeanYprime(), getSigmaP(), getSigmaY(), getSigmaYprime(), getUsedEve(), init(), LinRegBevPeb(), LinRegBevPeb(), operator+=(), operator+=(), operator<<, printSummaryAlphas(), printSummaryP(), printSummaryY(), printSummaryYP(), and solve().

mMP

TVectorD LinRegBevPeb::mMP

private

mean values

Definition at line 60 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), getMeanP(), init(), operator+=(), operator+=(), print(), and solve().

mMY

TVectorD LinRegBevPeb::mMY

private

Definition at line 60 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), getMeanY(), init(), operator+=(), operator+=(), print(), printSummaryMeansWithUnc(), and solve().

mMYp

TVectorD LinRegBevPeb::mMYp

private

Definition at line 60 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), getMeanYprime(), init(), printSummaryMeansWithUncCorrected(), and solve().

mRPP

TMatrixD LinRegBevPeb::mRPP

private

Definition at line 40 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mRPY

TMatrixD LinRegBevPeb::mRPY

private

correlations

Definition at line 39 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mRYP

TMatrixD LinRegBevPeb::mRYP

private

Definition at line 39 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), and init().

mRYY

TMatrixD LinRegBevPeb::mRYY

private

Definition at line 40 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mRYYp

TMatrixD LinRegBevPeb::mRYYp

private

Definition at line 41 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mSP

TVectorD LinRegBevPeb::mSP

private

sigmas

Definition at line 56 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mSPP

TMatrixD LinRegBevPeb::mSPP

private

Definition at line 53 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mSPY

TMatrixD LinRegBevPeb::mSPY

private

normalized covariances

Definition at line 52 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mSY

TVectorD LinRegBevPeb::mSY

private

Definition at line 56 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), printSummaryMeansWithUnc(), and solve().

mSYP

TMatrixD LinRegBevPeb::mSYP

private

Definition at line 52 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mSYp

TVectorD LinRegBevPeb::mSYp

private

Definition at line 57 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), printSummaryMeansWithUncCorrected(), and solve().

mSYY

TMatrixD LinRegBevPeb::mSYY

private

Definition at line 53 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mSYYp

TMatrixD LinRegBevPeb::mSYYp

private

Definition at line 54 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mVP

TVectorD LinRegBevPeb::mVP

private

variances

Definition at line 48 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mVPP

TMatrixD LinRegBevPeb::mVPP

private

Definition at line 45 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), getCovarianceP(), getSigmaP(), init(), operator+=(), operator+=(), print(), and solve().

mVPY

TMatrixD LinRegBevPeb::mVPY

private

unnormalized covariances

Definition at line 44 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), getCovariancePY(), init(), operator+=(), operator+=(), print(), and solve().

mVY

TVectorD LinRegBevPeb::mVY

private

Definition at line 48 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mVYP

TMatrixD LinRegBevPeb::mVYP

private

Definition at line 44 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mVYp

TVectorD LinRegBevPeb::mVYp

private

Definition at line 49 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), init(), and solve().

mVYY

TMatrixD LinRegBevPeb::mVYY

private

Definition at line 45 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), getCovarianceY(), getSigmaY(), init(), operator+=(), operator+=(), print(), and solve().

mVYYp

TMatrixD LinRegBevPeb::mVYYp

private

Definition at line 46 of file LinReg_Bevington_Pebay.h.

Referenced by clear(), getSigmaYprime(), init(), print(), and solve().

nP

int LinRegBevPeb::nP

private

Definition at line 31 of file LinReg_Bevington_Pebay.h.

Referenced by failed(), getCovarianceP(), getCovariancePY(), getMeanP(), getSigmaP(), init(), LinRegBevPeb(), LinRegBevPeb(), print(), printSummaryAlphas(), printSummaryP(), printSummaryYP(), setDims(), and solve().

nY

int LinRegBevPeb::nY

private

Definition at line 32 of file LinReg_Bevington_Pebay.h.

Referenced by getCovariancePY(), getCovarianceY(), getMeanY(), getMeanYprime(), getSigmaY(), getSigmaYprime(), init(), LinRegBevPeb(), LinRegBevPeb(), print(), printSummaryAlphas(), printSummaryMeansWithUnc(), printSummaryMeansWithUncCorrected(), printSummaryY(), printSummaryYP(), and setDims().

---

The documentation for this class was generated from the following files:

• LinReg_Bevington_Pebay.h
• LinReg_Bevington_Pebay.cc