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
#	Content
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	}