Teddy

Description

class Teddy

Use MC techniques to apply a function \(f\) to a distribution $x$ with (non-diagonal) covariance matrix \(\mathbf{V}\).

Procedure:

1. Diagonalise the covariance matrix
2. Generate a multi-dimensional Gaussian event \(\delta\) distributed according to the (sqrt of the) eigenvalues \(\epsilon_i\)

   \[ \delta_i \sim \mathcal{N}(0,\sqrt{\epsilon_i}) \]

3. Transform back to the original base (where the cov matrix is non-diagonal)
4. Add rotated vector to the input distribution an apply the function \(f\):

   \[ \mathbf{y}_k = \]

5. Go back to 2. until a satisfying statistics is reached, with Teddy::play() (the criterion of satisfaction and the loop are to be defined outside of the Teddy instance)
6. Get output histogram and covariance matrix with Teddy::buy()

#include <Teddy.h>

Public Member Functions
	Teddy (const Eigen::VectorXd &, const Eigen::MatrixXd &, std::function< Eigen::VectorXd(const Eigen::VectorXd &)>=identity, int=42)
	Teddy (const std::unique_ptr< TH1D > &, const std::unique_ptr< TH2D > &, std::function< Eigen::VectorXd(const Eigen::VectorXd &)>=identity, int=42)
void	play ()
std::pair< Eigen::VectorXd, Eigen::MatrixXd >	buy () const
template<typename... Args>
std::pair< std::unique_ptr< TH1D >, std::unique_ptr< TH2D > >	buy (TString name, Args... args) const

Static Public Attributes
static bool	verbose = true

Static Private Member Functions
static Eigen::VectorXd	H2V (const std::unique_ptr< TH1D > &)
static Eigen::MatrixXd	H2M (const std::unique_ptr< TH2D > &)
static int	GuessOutDim (int, std::function< Eigen::VectorXd(const Eigen::VectorXd &)>)
static void	SanityCheck (const Eigen::VectorXd &, Eigen::VectorXd &, double=0)

Private Attributes
const Eigen::VectorXd	vec
const Eigen::MatrixXd	cov
const int	nOldBins
const int	nNewBins
long long	N
Eigen::VectorXd	eigVals
Eigen::MatrixXd	rotation
Eigen::VectorXd	sigma
Eigen::VectorXd	sum
Eigen::MatrixXd	sum2
std::function< Eigen::VectorXd(const Eigen::VectorXd &)>	func
TRandom3	random

Constructor & Destructor Documentation

Teddy() [1/2]

Teddy	(	const Eigen::VectorXd &	,
		const Eigen::MatrixXd &	,
		std::function< Eigen::VectorXd(const Eigen::VectorXd &)>	= identity,
		int	= 42
	)

The rotation and the eigenvalues are determined from the input covariance matrix. Default function is the identity (useful for testing).

Teddy() [2/2]

Teddy	(	const std::unique_ptr< TH1D > &	,
		const std::unique_ptr< TH2D > &	,
		std::function< Eigen::VectorXd(const Eigen::VectorXd &)>	= identity,
		int	= 42
	)

The rotation and the eigenvalues are determined from the input covariance matrix. Default function is the identity (useful for testing).

Member Function Documentation

buy() [1/2]

pair< VectorXd, MatrixXd > buy

(

)

const

Get output distributions, can be called several times.

{
    VectorXd outputDist = sum*(1./N);
 
    MatrixXd avSq = sum2*(1./N);
    MatrixXd outputCov = avSq - outputDist * outputDist.transpose();
 
    return {outputDist, outputCov};
}

buy() [2/2]

std::pair<std::unique_ptr<TH1D>,std::unique_ptr<TH2D> > buy ( TString  name, Args...  args  ) const

inline

Get output distributions, can be called several times. Input histrogams are defined from scratch and filled. The statistical uncertainty of the 1D histogram correspondings to the square-root of the diagonal elements of the covariance matrix

    {
        using namespace std;
        using namespace Eigen;
 
        TAxis * axis = new TAxis(args...);
        int ndiv = axis->GetNbins();
        double * edges = new double[ndiv];
        axis->GetLowEdge(edges);
        edges[ndiv-1] = axis->GetBinUpEdge(ndiv);
        auto h   = make_unique<TH1D>(name         , name, ndiv-1, edges);
        auto cov = make_unique<TH2D>(name + "_cov", name, ndiv-1, edges, ndiv-1, edges);
 
        VectorXd x(nNewBins);
        MatrixXd err(nNewBins, nNewBins);
        tie(x,err) = buy();
    
        for (int i=0; i<nNewBins; ++i) {
            double content = x(i),
                   error = sqrt(err(i,i));
    
            h->SetBinContent(i+1, content);
            h->SetBinError  (i+1, error  );
    
            for (int j=0; j<nNewBins; ++j)
                cov->SetBinContent(i+1,j+1, err(i,j));
        }
 
        return make_pair<unique_ptr<TH1D>,unique_ptr<TH2D>>(move(h),move(cov));
    }

GuessOutDim()

int GuessOutDim ( int  , std::function< Eigen::VectorXd(const Eigen::VectorXd &)>    )

staticprivate

Guess size of the output distribution by calling once the function with dummy input.

{
    if (verbose) cout << __FILE__ << ':' << __LINE__ << '\t' << __func__ << endl;
    VectorXd In(nOldBins);
    VectorXd Out = f(In); // guessing the size from the output of the function
    return Out.size();
}

H2M()

MatrixXd H2M ( const std::unique_ptr< TH2D > &  h )

staticprivate

conversion from (2D) histogram to matrix

{
    if (verbose) cout << __FILE__ << ':' << __LINE__ << '\t' << __func__ << endl;
    int sizx = h->GetNbinsX(),
        sizy = h->GetNbinsY();
    MatrixXd mat(sizx,sizy);
    for(int j = 1; j <= sizx; ++j)
    for(int i = 1; i <= sizy; ++i)
        mat(i-1,j-1) = h->GetBinContent(j,i);
    return mat;
}

H2V()

VectorXd H2V ( const std::unique_ptr< TH1D > &  h )

staticprivate

conversion from (1D) histogram to vector

{
    if (verbose) cout << __FILE__ << ':' << __LINE__ << '\t' << __func__ << endl;
    VectorXd v(h->GetNbinsX());
    for(int i = 1; i <= h->GetNbinsX(); ++i)
        v(i-1) = h->GetBinContent(i);
    return v;
}

play()

void play

(

)

Needs to be called from a loop, but requires no argument.

{
    // 1) get vector in rotated basis
    VectorXd deltaP(nOldBins); // P for prime
    for (int k = 0; k < nOldBins; ++k) 
        deltaP(k) = random.Gaus(0, sigma(k));
 
    // 2) rotate the shift back
    VectorXd delta = rotation*deltaP;
    VectorXd x = vec + delta;
 
    // 3) compute whatever transformation you like on x
    VectorXd y = func(x); // dimension may change here!
    sum+= y;
 
    // 4) get its uncertainty
    MatrixXd square = y * y.transpose(); // tensor product
    sum2 += square;
 
    // 5) increment count for normalisation in Teddy::buy()
    ++N;
}

SanityCheck()

void SanityCheck ( const Eigen::VectorXd &  , Eigen::VectorXd &  , double  = 0  )

staticprivate

Check the number of negative eigenvalues and fail if any is found.

{
    int nEigv = 0, nNegEigVals = 0, nEmptyBins = 0;
    for (long k = 0; k < eigVals.size(); ++k) {
        if (eigVals(k) > -tolerance) {
            eigVals(k) = max(eigVals(k), 0.);
            ++nEigv;
        }
        else {
            ++nNegEigVals;
            cout << __FILE__ << ':' << __LINE__ << '\t' << k << ' ' << eigVals(k) << '\n';
        }
 
        if (x(k) == 0)
            ++nEmptyBins;
    }
    cout << flush;
 
    if (nNegEigVals > 0)
        BOOST_THROW_EXCEPTION( invalid_argument(Form("%d negative eigenvalues (%d empty bins in input dist)",
                                                    nNegEigVals,                nEmptyBins)) );
}

Member Data Documentation

cov

const Eigen::MatrixXd cov

private

covariance matrix

eigVals

Eigen::VectorXd eigVals

private

eigenvalues of input covariance matrix

func

std::function<Eigen::VectorXd(const Eigen::VectorXd &)> func

private

transformation defining output as a function of the input

N

long long N

private

number of calls

nNewBins

const int nNewBins

private

number of bins of output distributions

nOldBins

const int nOldBins

private

number of bins of input distributions

random

TRandom3 random

private

random generator

rotation

Eigen::MatrixXd rotation

private

rotation from diagonal basis to original basis

sigma

Eigen::VectorXd sigma

private

Gaussian widths in diagonal base.

sum

Eigen::VectorXd sum

private

sum of events, increased at each call and from which output distribution is estimated

sum2

Eigen::MatrixXd sum2

private

sum of tensor products of each event, increased at each call and from which output covariance is estimated

vec

const Eigen::VectorXd vec

private

distribution

verbose

bool verbose = true

static

---

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

• /builds/cms-analysis/general/DasAnalysisSystem/Core/Installer/Core/Uncertainties/interface/Teddy.h
• /builds/cms-analysis/general/DasAnalysisSystem/Core/Installer/Core/Uncertainties/src/Teddy.cc