Stops/BATLimits/statistics.cc

#include "statistics.hh"

#include <sstream>
#include <cassert>

#include <TGraph.h>
#include <TMath.h>
#include <TH1D.h>
#include <TF1.h>

int RandomPrior::counter_=0;

std::pair<double, double> evaluateInterval(TGraph* posterior, double alpha, double leftsidetail)
{
  double lowerCutOff = leftsidetail * alpha;
  double upperCutOff = 1. - (1.- leftsidetail) * alpha;

  double upper = 0, lower = 0;

  // normalize the interval, first
  double normalization=0.0;
  for(int i=0; i<posterior->GetN()-1; i++) {
    double firstx, firsty;
    double nextx, nexty;
    posterior->GetPoint(i, firstx, firsty);
    posterior->GetPoint(i+1, nextx, nexty);

    double intervalIntegral=(nextx-firstx)*0.5*(firsty+nexty);
    normalization+=intervalIntegral;
  }

  // now compute the intervals
  double integral=0.0;
  for(int i=0; i<posterior->GetN()-1; i++) {
    double firstx, firsty;
    double nextx, nexty;
    posterior->GetPoint(i, firstx, firsty);
    posterior->GetPoint(i+1, nextx, nexty);

    double intervalIntegral=(nextx-firstx)*0.5*(firsty+nexty)/normalization;

    // interpolate lower
    if(integral<=lowerCutOff && (integral+intervalIntegral)>=lowerCutOff) {
      lower=firstx;
    }
    if(integral<=upperCutOff && (integral+intervalIntegral)>=upperCutOff) {
      double m=(nexty-firsty)/(nextx-firstx);
      upper = firstx+(-firsty+sqrt(firsty*firsty+2*m*(upperCutOff-integral)*normalization))/m;
    }
    integral+=intervalIntegral;
  }

  std::pair<double, double> p(lower, upper);
  return p;
}

void getQuantiles(std::vector<double>& limits, double &median_, std::pair<double, double>& onesigma_, std::pair<double, double>& twosigma_) {
  unsigned int nit=limits.size();
  if(nit==0) return;

  // sort the vector with limits
  std::sort(limits.begin(), limits.end());

  // median for the expected limit
  median_ = TMath::Median(nit, &limits[0]);

  // quantiles for the expected limit bands
  double prob[4]; // array with quantile boundaries
  prob[0] = 0.021;
  prob[1] = 0.159;
  prob[2] = 0.841;
  prob[3] = 0.979;

  // array for the results
  double quantiles[4];
  TMath::Quantiles(nit, 4, &limits[0], quantiles, prob); // evaluate quantiles
  onesigma_.first=quantiles[1];
  onesigma_.second=quantiles[2];
  twosigma_.first=quantiles[0];
  twosigma_.second=quantiles[3];

  return;
}

double lognormal(double *x, double *par)
{
  if(par[0]<0.0) {
    std::cout << "par[0] = " << par[0] << std::endl;
    assert(0);
  }
  if(par[1]<0.0) {
    std::cout << "par[1] = " << par[1] << std::endl;
    assert(0);
  }

  double m0=par[0];
  double k=par[1]/par[0]+1.;
  double s=TMath::Log(k);
  return TMath::LogNormal(x[0], s, 0.0, m0);
}

double gaussian(double *x, double *par)
{
  if(par[1]<0.0) {
    std::cout << "par[1] = " << par[1] << std::endl;
    assert(0);
  }

  double m0=par[0];
  double s=par[1];
  return TMath::Gaus(x[0],m0,s,kTRUE);
}

double gamma(double *x, double *par)
{
  if(par[0]<0.0) {
    std::cout << "par[0] = " << par[0] << std::endl;
    assert(0);
  }
  if(par[1]<0.0) {
    std::cout << "par[1] = " << par[1] << std::endl;
    assert(0);
  }

  double s=par[0]*par[0]/par[1]/par[1]+1.;
  double tau = par[0]/par[1]/par[1];
  return TMath::GammaDist(x[0], s, 0.0, 1./tau);
}

RandomPrior::RandomPrior(int priorType, double median, double variance, double min, double max)
{
  // set the prior type
  priorType_=priorType;

  // create function
  std::ostringstream oss;
  if(priorType==1) // Lognormal
  {
    oss << "_Random_Lognormal__priorfcn_" << (counter_++);
    priorfcn_ = new TF1(oss.str().c_str(), lognormal, std::max(0., (min<(median-5*variance) ? (median-5*variance) : min)), (max>(median+5*variance) ? (median+5*variance) : max), 2);
    priorfcn_->SetParameter(0, median);
    priorfcn_->SetParameter(1, variance);
  }
  else if(priorType==2) // Gaussian
  {
    oss << "_Random_Gaussian__priorfcn_" << (counter_++);
    priorfcn_ = new TF1(oss.str().c_str(), gaussian, (min<(median-5*variance) ? (median-5*variance) : min), (max>(median+5*variance) ? (median+5*variance) : max), 2);
    priorfcn_->SetParameter(0, median);
    priorfcn_->SetParameter(1, variance);
  }
  else if(priorType==3) // Gamma
  {
    oss << "_Random_Gamma__priorfcn_" << (counter_++);
    priorfcn_ = new TF1(oss.str().c_str(), gamma, std::max(0., (min<(median-5*variance) ? (median-5*variance) : min)), (max>(median+5*variance) ? (median+5*variance) : max), 2);
    priorfcn_->SetParameter(0, median);
    priorfcn_->SetParameter(1, variance);
  }
  else // Uniform
  {
    oss << "_Random_Uniform__priorfcn_" << (counter_++);
    priorfcn_ = new TF1(oss.str().c_str(), "pol0", min, max);
    priorfcn_->SetParameter(0, 1./(max-min));
  }
}

RandomPrior::~RandomPrior()
{
  delete priorfcn_;
}

double RandomPrior::getRandom(void) const
{
  return priorfcn_->GetRandom();
}

double RandomPrior::getXmin(void) const
{
  return priorfcn_->GetXmin();
}

double RandomPrior::getXmax(void) const
{
  return priorfcn_->GetXmax();
}

double RandomPrior::eval(double x) const
{
  return priorfcn_->Eval(x);
}

int RandomPrior::getPriorType(void) const
{
  return priorType_;
}
Revision:	1.1
Committed:	Thu Sep 5 03:58:17 2013 UTC (11 years, 8 months ago) by algomez
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	First version
#	User	Rev	Content
1	algomez	1.1	#include "statistics.hh"
2
3			#include <sstream>
4			#include <cassert>
5
6			#include <TGraph.h>
7			#include <TMath.h>
8			#include <TH1D.h>
9			#include <TF1.h>
10
11			int RandomPrior::counter_=0;
12
13			std::pair<double, double> evaluateInterval(TGraph* posterior, double alpha, double leftsidetail)
14			{
15			double lowerCutOff = leftsidetail * alpha;
16			double upperCutOff = 1. - (1.- leftsidetail) * alpha;
17
18			double upper = 0, lower = 0;
19
20			// normalize the interval, first
21			double normalization=0.0;
22			for(int i=0; i<posterior->GetN()-1; i++) {
23			double firstx, firsty;
24			double nextx, nexty;
25			posterior->GetPoint(i, firstx, firsty);
26			posterior->GetPoint(i+1, nextx, nexty);
27
28			double intervalIntegral=(nextx-firstx)0.5(firsty+nexty);
29			normalization+=intervalIntegral;
30			}
31
32			// now compute the intervals
33			double integral=0.0;
34			for(int i=0; i<posterior->GetN()-1; i++) {
35			double firstx, firsty;
36			double nextx, nexty;
37			posterior->GetPoint(i, firstx, firsty);
38			posterior->GetPoint(i+1, nextx, nexty);
39
40			double intervalIntegral=(nextx-firstx)0.5(firsty+nexty)/normalization;
41
42			// interpolate lower
43			if(integral<=lowerCutOff && (integral+intervalIntegral)>=lowerCutOff) {
44			lower=firstx;
45			}
46			if(integral<=upperCutOff && (integral+intervalIntegral)>=upperCutOff) {
47			double m=(nexty-firsty)/(nextx-firstx);
48			upper = firstx+(-firsty+sqrt(firstyfirsty+2m(upperCutOff-integral)normalization))/m;
49			}
50			integral+=intervalIntegral;
51			}
52
53			std::pair<double, double> p(lower, upper);
54			return p;
55			}
56
57			void getQuantiles(std::vector<double>& limits, double &median_, std::pair<double, double>& onesigma_, std::pair<double, double>& twosigma_) {
58			unsigned int nit=limits.size();
59			if(nit==0) return;
60
61			// sort the vector with limits
62			std::sort(limits.begin(), limits.end());
63
64			// median for the expected limit
65			median_ = TMath::Median(nit, &limits[0]);
66
67			// quantiles for the expected limit bands
68			double prob[4]; // array with quantile boundaries
69			prob[0] = 0.021;
70			prob[1] = 0.159;
71			prob[2] = 0.841;
72			prob[3] = 0.979;
73
74			// array for the results
75			double quantiles[4];
76			TMath::Quantiles(nit, 4, &limits[0], quantiles, prob); // evaluate quantiles
77			onesigma_.first=quantiles[1];
78			onesigma_.second=quantiles[2];
79			twosigma_.first=quantiles[0];
80			twosigma_.second=quantiles[3];
81
82			return;
83			}
84
85			double lognormal(double x, double par)
86			{
87			if(par[0]<0.0) {
88			std::cout << "par[0] = " << par[0] << std::endl;
89			assert(0);
90			}
91			if(par[1]<0.0) {
92			std::cout << "par[1] = " << par[1] << std::endl;
93			assert(0);
94			}
95
96			double m0=par[0];
97			double k=par[1]/par[0]+1.;
98			double s=TMath::Log(k);
99			return TMath::LogNormal(x[0], s, 0.0, m0);
100			}
101
102			double gaussian(double x, double par)
103			{
104			if(par[1]<0.0) {
105			std::cout << "par[1] = " << par[1] << std::endl;
106			assert(0);
107			}
108
109			double m0=par[0];
110			double s=par[1];
111			return TMath::Gaus(x[0],m0,s,kTRUE);
112			}
113
114			double gamma(double x, double par)
115			{
116			if(par[0]<0.0) {
117			std::cout << "par[0] = " << par[0] << std::endl;
118			assert(0);
119			}
120			if(par[1]<0.0) {
121			std::cout << "par[1] = " << par[1] << std::endl;
122			assert(0);
123			}
124
125			double s=par[0]*par[0]/par[1]/par[1]+1.;
126			double tau = par[0]/par[1]/par[1];
127			return TMath::GammaDist(x[0], s, 0.0, 1./tau);
128			}
129
130			RandomPrior::RandomPrior(int priorType, double median, double variance, double min, double max)
131			{
132			// set the prior type
133			priorType_=priorType;
134
135			// create function
136			std::ostringstream oss;
137			if(priorType==1) // Lognormal
138			{
139			oss << "_Random_Lognormal__priorfcn_" << (counter_++);
140			priorfcn_ = new TF1(oss.str().c_str(), lognormal, std::max(0., (min<(median-5variance) ? (median-5variance) : min)), (max>(median+5variance) ? (median+5variance) : max), 2);
141			priorfcn_->SetParameter(0, median);
142			priorfcn_->SetParameter(1, variance);
143			}
144			else if(priorType==2) // Gaussian
145			{
146			oss << "_Random_Gaussian__priorfcn_" << (counter_++);
147			priorfcn_ = new TF1(oss.str().c_str(), gaussian, (min<(median-5variance) ? (median-5variance) : min), (max>(median+5variance) ? (median+5variance) : max), 2);
148			priorfcn_->SetParameter(0, median);
149			priorfcn_->SetParameter(1, variance);
150			}
151			else if(priorType==3) // Gamma
152			{
153			oss << "_Random_Gamma__priorfcn_" << (counter_++);
154			priorfcn_ = new TF1(oss.str().c_str(), gamma, std::max(0., (min<(median-5variance) ? (median-5variance) : min)), (max>(median+5variance) ? (median+5variance) : max), 2);
155			priorfcn_->SetParameter(0, median);
156			priorfcn_->SetParameter(1, variance);
157			}
158			else // Uniform
159			{
160			oss << "_Random_Uniform__priorfcn_" << (counter_++);
161			priorfcn_ = new TF1(oss.str().c_str(), "pol0", min, max);
162			priorfcn_->SetParameter(0, 1./(max-min));
163			}
164			}
165
166			RandomPrior::~RandomPrior()
167			{
168			delete priorfcn_;
169			}
170
171			double RandomPrior::getRandom(void) const
172			{
173			return priorfcn_->GetRandom();
174			}
175
176			double RandomPrior::getXmin(void) const
177			{
178			return priorfcn_->GetXmin();
179			}
180
181			double RandomPrior::getXmax(void) const
182			{
183			return priorfcn_->GetXmax();
184			}
185
186			double RandomPrior::eval(double x) const
187			{
188			return priorfcn_->Eval(x);
189			}
190
191			int RandomPrior::getPriorType(void) const
192			{
193			return priorType_;
194			}