/*
Fisher and Mahalanobis Discriminants (Linear Discriminant Analysis)
<p>
In the method of Fisher discriminants event selection is performed
in a transformed variable space with zero linear correlations, by
distinguishing the mean values of the signal and background
distributions.<br></p>
<p>
The linear discriminant analysis determines an axis in the (correlated)
hyperspace of the input variables
such that, when projecting the output classes (signal and background)
upon this axis, they are pushed as far as possible away from each other,
while events of a same class are confined in a close vicinity.
The linearity property of this method is reflected in the metric with
which "far apart" and "close vicinity" are determined: the covariance
matrix of the discriminant variable space.
</p>
<p>
The classification of the events in signal and background classes
relies on the following characteristics (only): overall sample means,
<i><my:o>x</my:o><sub>i</sub></i>, for each input variable, <i>i</i>,
class-specific sample means, <i><my:o>x</my:o><sub>S(B),i</sub></i>,
and total covariance matrix <i>T<sub>ij</sub></i>. The covariance matrix
can be decomposed into the sum of a <i>within-</i> (<i>W<sub>ij</sub></i>)
and a <i>between-class</i> (<i>B<sub>ij</sub></i>) class matrix. They describe
the dispersion of events relative to the means of their own class (within-class
matrix), and relative to the overall sample means (between-class matrix).
The Fisher coefficients, <i>F<sub>i</sub></i>, are then given by <br>
<center>
<img vspace=6 src="gif/tmva_fisherC.gif" align="bottom" >
</center>
where in TMVA is set <i>N<sub>S</sub>=N<sub>B</sub></i>, so that the factor
in front of the sum simplifies to ½.
The Fisher discriminant then reads<br>
<center>
<img vspace=6 src="gif/tmva_fisherD.gif" align="bottom" >
</center>
The offset <i>F</i><sub>0</sub> centers the sample mean of <i>x</i><sub>Fi</sub>
at zero. Instead of using the within-class matrix, the Mahalanobis variant
determines the Fisher coefficients as follows:<br>
<center>
<img vspace=6 src="gif/tmva_mahaC.gif" align="bottom" >
</center>
with resulting <i>x</i><sub>Ma</sub> that are very similar to the
<i>x</i><sub>Fi</sub>. <br></p>
TMVA provides two outputs for the ranking of the input variables:<br><p></p>
<ul>
<li> <u>Fisher test:</u> the Fisher analysis aims at simultaneously maximising
the between-class separation, while minimising the within-class dispersion.
A useful measure of the discrimination power of a variable is hence given
by the diagonal quantity: <i>B<sub>ii</sub>/W<sub>ii</sub></i>.
</li>
<li> <u>Discrimination power:</u> the value of the Fisher coefficient is a
measure of the discriminating power of a variable. The discrimination power
of set of input variables can therefore be measured by the scalar
<center>
<img vspace=6 src="gif/tmva_discpower.gif" align="bottom" >
</center>
</li>
</ul>
The corresponding numbers are printed on standard output.
*/
// End_Html
#include <assert.h>
#include "TMath.h"
#include "TMVA/MethodFisher.h"
#include "TMVA/Tools.h"
#include "TMatrix.h"
#include "TMVA/Ranking.h"
ClassImp(TMVA::MethodFisher)
TMVA::MethodFisher::MethodFisher( TString jobName, TString methodTitle, DataSet& theData,
TString theOption, TDirectory* theTargetDir )
: TMVA::MethodBase( jobName, methodTitle, theData, theOption, theTargetDir )
{
InitFisher();
DeclareOptions();
ParseOptions();
ProcessOptions();
if (HasTrainingTree()) InitMatrices();
}
TMVA::MethodFisher::MethodFisher( DataSet& theData,
TString theWeightFile,
TDirectory* theTargetDir )
: TMVA::MethodBase( theData, theWeightFile, theTargetDir )
{
InitFisher();
DeclareOptions();
}
void TMVA::MethodFisher::InitFisher( void )
{
SetMethodName( "Fisher" );
SetMethodType( TMVA::Types::kFisher );
SetTestvarName();
SetNormalised( kTRUE );
fMeanMatx = 0;
fBetw = 0;
fWith = 0;
fCov = 0;
fSumOfWeightsS = fSumOfWeightsB = 0;
fF0 = 0;
fFisherCoeff = new vector<Double_t>( GetNvar() );
SetSignalReferenceCut( 0.0 );
}
void TMVA::MethodFisher::DeclareOptions()
{
DeclareOptionRef( fTheMethod = "Fisher", "Method", "Discrimination method" );
AddPreDefVal(TString("Fisher"));
AddPreDefVal(TString("Mahalanobis"));
}
void TMVA::MethodFisher::ProcessOptions()
{
MethodBase::ProcessOptions();
if (fTheMethod == "Fisher" ) fFisherMethod = kFisher;
else fFisherMethod = kMahalanobis;
PrintOptions();
CheckForUnusedOptions();
}
TMVA::MethodFisher::~MethodFisher( void )
{
delete fBetw;
delete fWith;
delete fCov;
delete fDiscrimPow;
delete fFisherCoeff;
}
void TMVA::MethodFisher::Train( void )
{
if (!CheckSanity()) fLogger << kFATAL << "<Train> sanity check failed" << Endl;
GetMean();
GetCov_WithinClass();
GetCov_BetweenClass();
GetCov_Full();
GetFisherCoeff();
GetDiscrimPower();
PrintCoefficients();
}
Double_t TMVA::MethodFisher::GetMvaValue()
{
Double_t result = fF0;
for (Int_t ivar=0; ivar<GetNvar(); ivar++)
result += (*fFisherCoeff)[ivar]*GetEventVal(ivar);
return result;
}
void TMVA::MethodFisher::InitMatrices( void )
{
if (!HasTrainingTree()) {
fLogger << kFATAL << "<InitMatrices> fatal error: Data().TrainingTree() is zero pointer" << Endl;
}
fMeanMatx = new TMatrixD( GetNvar(), 3 );
fBetw = new TMatrixD( GetNvar(), GetNvar() );
fWith = new TMatrixD( GetNvar(), GetNvar() );
fCov = new TMatrixD( GetNvar(), GetNvar() );
fDiscrimPow = new vector<Double_t>( GetNvar() );
}
void TMVA::MethodFisher::GetMean( void )
{
fSumOfWeightsS = 0;
fSumOfWeightsB = 0;
for (Int_t ievt=0; ievt<Data().GetNEvtTrain(); ievt++) {
ReadTrainingEvent(ievt);
if (IsSignalEvent()) fSumOfWeightsS += GetEventWeight();
else fSumOfWeightsB += GetEventWeight();
}
Double_t *sumS = new Double_t[(const Int_t)GetNvar()];
Double_t *sumB = new Double_t[(const Int_t)GetNvar()];
for (Int_t ivar=0; ivar<GetNvar(); ivar++) { sumS[ivar] = sumB[ivar] = 0; }
for (Int_t ievt=0; ievt<Data().GetNEvtTrain(); ievt++) {
ReadTrainingEvent(ievt);
Double_t weight = GetEventWeight();
if (IsSignalEvent()) fSumOfWeightsS += weight;
else fSumOfWeightsB += weight;
Double_t* sum = IsSignalEvent() ? sumS : sumB;
for (Int_t ivar=0; ivar<GetNvar(); ivar++) sum[ivar] += GetEventVal( ivar )*weight;
}
for (Int_t ivar=0; ivar<GetNvar(); ivar++) {
(*fMeanMatx)( ivar, 2 ) = sumS[ivar];
(*fMeanMatx)( ivar, 0 ) = sumS[ivar]/fSumOfWeightsS;
(*fMeanMatx)( ivar, 2 ) += sumB[ivar];
(*fMeanMatx)( ivar, 1 ) = sumB[ivar]/fSumOfWeightsB;
(*fMeanMatx)( ivar, 2 ) /= (fSumOfWeightsS + fSumOfWeightsB);
}
delete [] sumS;
delete [] sumB;
}
void TMVA::MethodFisher::GetCov_WithinClass( void )
{
assert( fSumOfWeightsS > 0 && fSumOfWeightsB > 0 );
const Int_t nvar = GetNvar();
const Int_t nvar2 = nvar*nvar;
Double_t *sumSig = new Double_t[nvar2];
Double_t *sumBgd = new Double_t[nvar2];
Double_t *xval = new Double_t[nvar];
memset(sumSig,0,nvar2*sizeof(Double_t));
memset(sumBgd,0,nvar2*sizeof(Double_t));
Int_t k=0;
for (Int_t ievt=0; ievt<Data().GetNEvtTrain(); ievt++) {
ReadTrainingEvent(ievt);
Double_t weight = GetEventWeight();
for (Int_t x=0; x<nvar; x++) xval[x] = GetEventVal( x );
Int_t k=0;
for (Int_t x=0; x<nvar; x++) {
for (Int_t y=0; y<nvar; y++) {
Double_t v = ( (xval[x] - (*fMeanMatx)(x, 0))*(xval[y] - (*fMeanMatx)(y, 0)) )*weight;
if (IsSignalEvent()) sumSig[k] += v;
else sumBgd[k] += v;
k++;
}
}
}
k=0;
for (Int_t x=0; x<nvar; x++) {
for (Int_t y=0; y<nvar; y++) {
(*fWith)(x, y) = (sumSig[k] + sumBgd[k])/(fSumOfWeightsS + fSumOfWeightsB);
k++;
}
}
delete [] sumSig;
delete [] sumBgd;
delete [] xval;
}
void TMVA::MethodFisher::GetCov_BetweenClass( void )
{
assert( fSumOfWeightsS > 0 && fSumOfWeightsB > 0);
Double_t prodSig, prodBgd;
for (Int_t x=0; x<GetNvar(); x++) {
for (Int_t y=0; y<GetNvar(); y++) {
prodSig = ( ((*fMeanMatx)(x, 0) - (*fMeanMatx)(x, 2))*
((*fMeanMatx)(y, 0) - (*fMeanMatx)(y, 2)) );
prodBgd = ( ((*fMeanMatx)(x, 1) - (*fMeanMatx)(x, 2))*
((*fMeanMatx)(y, 1) - (*fMeanMatx)(y, 2)) );
(*fBetw)(x, y) = (fSumOfWeightsS*prodSig + fSumOfWeightsB*prodBgd) / (fSumOfWeightsS + fSumOfWeightsB);
}
}
}
void TMVA::MethodFisher::GetCov_Full( void )
{
for (Int_t x=0; x<GetNvar(); x++)
for (Int_t y=0; y<GetNvar(); y++)
(*fCov)(x, y) = (*fWith)(x, y) + (*fBetw)(x, y);
}
void TMVA::MethodFisher::GetFisherCoeff( void )
{
assert( fSumOfWeightsS > 0 && fSumOfWeightsB > 0);
TMatrixD* theMat = 0;
switch (GetFisherMethod()) {
case kFisher:
theMat = fWith;
break;
case kMahalanobis:
theMat = fCov;
break;
default:
fLogger << kFATAL << "<GetFisherCoeff> undefined method" << GetFisherMethod() << Endl;
}
TMatrixD invCov( *theMat );
if ( TMath::Abs(invCov.Determinant()) < 10E-24 ) {
fLogger << kWARNING << "<GetFisherCoeff> matrix is almost singular with deterninant="
<< TMath::Abs(invCov.Determinant())
<< " did you use the variables that are linear combinations or highly correlated?"
<< Endl;
}
if ( TMath::Abs(invCov.Determinant()) < 10E-120 ) {
fLogger << kFATAL << "<GetFisherCoeff> matrix is singular with determinant="
<< TMath::Abs(invCov.Determinant())
<< " did you use the variables that are linear combinations?"
<< Endl;
}
invCov.Invert();
Double_t xfact = TMath::Sqrt( fSumOfWeightsS*fSumOfWeightsB ) / (fSumOfWeightsS + fSumOfWeightsB);
vector<Double_t> diffMeans( GetNvar() );
Int_t ivar, jvar;
for (ivar=0; ivar<GetNvar(); ivar++) {
(*fFisherCoeff)[ivar] = 0;
for (jvar=0; jvar<GetNvar(); jvar++) {
Double_t d = (*fMeanMatx)(jvar, 0) - (*fMeanMatx)(jvar, 1);
(*fFisherCoeff)[ivar] += invCov(ivar, jvar)*d;
}
(*fFisherCoeff)[ivar] *= xfact;
}
fF0 = 0.0;
for (ivar=0; ivar<GetNvar(); ivar++){
fF0 += (*fFisherCoeff)[ivar]*((*fMeanMatx)(ivar, 0) + (*fMeanMatx)(ivar, 1));
}
fF0 /= -2.0;
}
void TMVA::MethodFisher::GetDiscrimPower( void )
{
for (Int_t ivar=0; ivar<GetNvar(); ivar++) {
if ((*fCov)(ivar, ivar) != 0)
(*fDiscrimPow)[ivar] = (*fBetw)(ivar, ivar)/(*fCov)(ivar, ivar);
else
(*fDiscrimPow)[ivar] = 0;
}
}
const TMVA::Ranking* TMVA::MethodFisher::CreateRanking()
{
fRanking = new Ranking( GetName(), "Discr. power" );
for (Int_t ivar=0; ivar<GetNvar(); ivar++) {
fRanking->AddRank( *new Rank( GetInputExp(ivar), (*fDiscrimPow)[ivar] ) );
}
return fRanking;
}
void TMVA::MethodFisher::PrintCoefficients( void )
{
fLogger << kINFO << "Results for Fisher coefficients:" << Endl;
vector<TString> vars;
vector<Double_t> coeffs;
for (Int_t ivar=0; ivar<GetNvar(); ivar++) {
vars .push_back( GetInputExp(ivar) );
coeffs.push_back( (*fFisherCoeff)[ivar] );
}
vars .push_back( "(offset)" );
coeffs.push_back( fF0 );
TMVA::Tools::FormattedOutput( coeffs, vars, "Variable" , "Coefficient", fLogger );
}
void TMVA::MethodFisher::WriteWeightsToStream( ostream& o ) const
{
o << fF0 << endl;
for (Int_t ivar=0; ivar<GetNvar(); ivar++) o << setprecision(12) << (*fFisherCoeff)[ivar] << endl;
}
void TMVA::MethodFisher::ReadWeightsFromStream( istream& istr )
{
istr >> fF0;
for (Int_t ivar=0; ivar<GetNvar(); ivar++) istr >> (*fFisherCoeff)[ivar];
}
void TMVA::MethodFisher::MakeClassSpecific( std::ostream& fout, const TString& className ) const
{
fout << " double fFisher0;" << endl;
fout << " std::vector<double> fFisherCoefficients;" << endl;
fout << "};" << endl;
fout << "" << endl;
fout << "inline void " << className << "::Initialize() " << endl;
fout << "{" << endl;
fout << " fFisher0 = " << setprecision(12) << fF0 << ";" << endl;
for (Int_t ivar=0; ivar<GetNvar(); ivar++) {
fout << " fFisherCoefficients.push_back( " << setprecision(12) << (*fFisherCoeff)[ivar] << " );" << endl;
}
fout << endl;
fout << " // sanity check" << endl;
fout << " if (fFisherCoefficients.size() != fNvars) {" << endl;
fout << " std::cout << \"Problem in class \\\"\" << fClassName << \"\\\"::Initialize: mismatch in number of input values\"" << endl;
fout << " << fFisherCoefficients.size() << \" != \" << fNvars << std::endl;" << endl;
fout << " fStatusIsClean = false;" << endl;
fout << " } " << endl;
fout << "}" << endl;
fout << endl;
fout << "inline double " << className << "::GetMvaValue__( const std::vector<double>& inputValues ) const" << endl;
fout << "{" << endl;
fout << " double retval = fFisher0;" << endl;
fout << " for (size_t ivar = 0; ivar < fNvars; ivar++) {" << endl;
fout << " retval += fFisherCoefficients[ivar]*inputValues[ivar];" << endl;
fout << " }" << endl;
fout << endl;
fout << " return retval;" << endl;
fout << "}" << endl;
fout << endl;
fout << "// Clean up" << endl;
fout << "inline void " << className << "::Clear() " << endl;
fout << "{" << endl;
fout << " // clear coefficients" << endl;
fout << " fFisherCoefficients.clear(); " << endl;
fout << "}" << endl;
}
void TMVA::MethodFisher::GetHelpMessage() const
{
fLogger << Endl;
fLogger << Tools::Color("bold") << "--- Short description:" << Tools::Color("reset") << Endl;
fLogger << Endl;
fLogger << "Fisher discriminants select events by distinguishing the mean " << Endl;
fLogger << "values of the signal and background distributions in a trans- " << Endl;
fLogger << "formed variable space where linear correlations are removed." << Endl;
fLogger << Endl;
fLogger << " (More precisely: the \"linear discriminator\" determines" << Endl;
fLogger << " an axis in the (correlated) hyperspace of the input " << Endl;
fLogger << " variables such that, when projecting the output classes " << Endl;
fLogger << " (signal and background) upon this axis, they are pushed " << Endl;
fLogger << " as far as possible away from each other, while events" << Endl;
fLogger << " of a same class are confined in a close vicinity. The " << Endl;
fLogger << " linearity property of this classifier is reflected in the " << Endl;
fLogger << " metric with which \"far apart\" and \"close vicinity\" are " << Endl;
fLogger << " determined: the covariance matrix of the discriminating" << Endl;
fLogger << " variable space.)" << Endl;
fLogger << Endl;
fLogger << Tools::Color("bold") << "--- Performance optimisation:" << Tools::Color("reset") << Endl;
fLogger << Endl;
fLogger << "Optimal performance for Fisher discriminants is obtained for " << Endl;
fLogger << "linearly correlated Gaussian-distributed variables. Any deviation" << Endl;
fLogger << "from this ideal reduces the achievable separation power. In " << Endl;
fLogger << "particular, no discrimination at all is achieved for a variable" << Endl;
fLogger << "that has the same sample mean for signal and background, even if " << Endl;
fLogger << "the shapes of the distributions are very different. Thus, Fisher " << Endl;
fLogger << "discriminants often benefit from suitable transformations of the " << Endl;
fLogger << "input variables. For example, if a variable x in [-1,1] has a " << Endl;
fLogger << "a parabolic signal distributions, and a uniform background" << Endl;
fLogger << "distributions, their mean value is zero in both cases, leading " << Endl;
fLogger << "to no separation. The simple transformation x -> |x| renders this " << Endl;
fLogger << "variable powerful for the use in a Fisher discriminant." << Endl;
fLogger << Endl;
fLogger << Tools::Color("bold") << "--- Performance tuning via configuration options:" << Tools::Color("reset") << Endl;
fLogger << Endl;
fLogger << "<None>" << Endl;
}
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.