Some useful mathematic functions. More...
Functions | |
| double | LogGaus (double x, double mean=0, double sigma=1, bool norm=false) |
| double | LogPoisson (double x, double par) |
| double | ApproxBinomial (int n, int k, double p) |
| double | LogApproxBinomial (int n, int k, double p) |
| double | LogBinomFactor (int n, int k) |
| double | ApproxLogFact (double x) |
| double | LogNoverK (int n, int k) |
| double | LogFact (int n) |
| int | Max (int x, int y) |
| int | Max (unsigned int x, unsigned int y) |
| double | Max (double x, double y) |
| double | Max (float x, float y) |
| int | Min (int x, int y) |
| int | Min (unsigned int x, unsigned int y) |
| double | Min (double x, double y) |
| double | Min (float x, float y) |
| int | Nint (double x) |
| double | rms (int n, const double *a) |
| double | LogBreitWignerNonRel (double x, double mean, double Gamma, bool norm=false) |
| double | LogBreitWignerRel (double x, double mean, double Gamma) |
| double | LogChi2 (double x, int n) |
| double | LogVoigtian (double x, double sigma, double gamma) |
| void | RandomChi2 (std::vector< double > &randoms, int K) |
| TH1D * | ECDF (const std::vector< double > &data) |
| std::vector< int > | longestRuns (const std::vector< bool > &bitStream) |
| std::vector< double > | longestRunsChi2 (const std::vector< double > &yMeasured, const std::vector< double > &yExpected, const std::vector< double > &sigma) |
| double | longestRunFrequency (unsigned int longestObserved, unsigned int nTrials) |
| double | SplitGaussian (double *x, double *par) |
| void | CacheFactorial (unsigned int n) |
| double | chi2 (double *x, double *par) |
| double | longestRunFrequency (unsigned longestObserved, unsigned int nTrials) |
Variables | |
| static unsigned int | nCacheFact = 10000 |
| static double * | logfact = 0 |
Some useful mathematic functions.
| double BCMath::ApproxBinomial | ( | int | n, | |
| int | k, | |||
| double | p | |||
| ) |
Calculates Binomial probability using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function.
Definition at line 70 of file BCMath.cxx.
{
return exp(LogApproxBinomial(n, k, p));
}
| double BCMath::ApproxLogFact | ( | double | x | ) |
Calculates natural logarithm of the n-factorial (n!) using Srinivasa Ramanujan approximation log(n!) = n*log(n) - n + log(n*(1.+4.*n*(1.+2.*n)))/6. + log(PI)/2. if n > 20. If n <= 20 it uses BCMath::LogFact to calculate it exactly.
Definition at line 122 of file BCMath.cxx.
{
if (x > BCMath::nCacheFact)
return x * log(x) - x + log(x * (1. + 4. * x * (1. + 2. * x))) / 6. + log(M_PI) / 2.;
else
return LogFact((int) x);
}
| void BCMath::CacheFactorial | ( | unsigned int | n | ) |
Cache factorials for first
n integers. The cache is filled upon first call of LogFact(). Definition at line 166 of file BCMath.cxx.
{
nCacheFact = n;
}
| double BCMath::chi2 | ( | double * | x, | |
| double * | par | |||
| ) |
Definition at line 306 of file BCMath.cxx.
{
return ROOT::Math::chisquared_pdf(x[0], par[0]);
}
| TH1D * BCMath::ECDF | ( | const std::vector< double > & | data | ) |
Calculate the empirical cumulative distribution function for one dimensional data vector. For consistency, the ECDF of value smaller than the minimum observed (underflow bin) is zero, and for larger than maximum (overflow bin) it is one.
| data | the observations |
Definition at line 327 of file BCMath.cxx.
{
int N = data.size();
std::set<double> uniqueObservations;
// sort and filter out multiple instances
for (int i = 0; i < N; ++i)
uniqueObservations.insert(data[i]);
// extract lower edges for CDF histogram
int nUnique = uniqueObservations.size();
double lowerEdges[nUnique];
// traverse the set
std::set<double>::iterator iter;
int counter = 0;
for (iter = uniqueObservations.begin(); iter != uniqueObservations.end(); ++iter) {
lowerEdges[counter] = *iter;
counter++;
}
// create histogram where
// lower edge of first bin = min. data
// upper edge of last bin = max. data
TH1D * ECDF = new TH1D("ECDF", "Empirical cumulative distribution function", nUnique - 1, lowerEdges);
// fill the data in to find multiplicities
for (int i = 0; i < N; ++i)
ECDF -> Fill(data[i]);
// just in case, empty the underflow
ECDF -> SetBinContent(0, 0.0);
// construct the ecdf
for (int nBin = 1; nBin <= ECDF->GetNbinsX(); nBin++) {
double previousBin = ECDF -> GetBinContent(nBin - 1);
// BCLog::OutDebug(Form("n_%d = %.2f", nBin, ECDF -> GetBinContent(nBin) ));
// BCLog::OutDebug(Form("previous_%d = %.2f", nBin, previousBin));
double thisBin = ECDF -> GetBinContent(nBin) / double(N);
ECDF -> SetBinContent(nBin, thisBin + previousBin);
// the uncertainty is only correctly estimated in the model
ECDF -> SetBinError(nBin, 0.0);
}
// set the endpoint to 1, so all larger values are at CDF=1
ECDF -> SetBinContent(ECDF->GetNbinsX() + 1, 1.);
// adjust for nice plotting
ECDF -> SetMinimum(0.);
ECDF -> SetMaximum(1.);
return ECDF;
}
| double BCMath::LogApproxBinomial | ( | int | n, | |
| int | k, | |||
| double | p | |||
| ) |
Calculates natural logarithm of the Binomial probability using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function.
Definition at line 77 of file BCMath.cxx.
{
// check p
if (p == 0)
return -1e99;
else if (p == 1)
return 0;
// switch parameters if n < k
if (n < k) {
int a = n;
n = k;
k = a;
}
return LogBinomFactor(n, k) + (double) k * log(p) + (double) (n - k) * log(1. - p);
}
| double BCMath::LogBinomFactor | ( | int | n, | |
| int | k | |||
| ) |
Calculates natural logarithm of the Binomial factor "n over k" using approximations for factorial calculations if calculation for number greater than 20 required using the BCMath::ApproxLogFact function. Even for large numbers the calculation is performed precisely, if n-k < 5
Definition at line 98 of file BCMath.cxx.
{
// switch parameters if n < k
if (n < k) {
int a = n;
n = k;
k = a;
}
if (n == k || k == 0)
return 0.;
if (k == 1 || k == n - 1)
return log((double) n);
// if no approximation needed
if ( n < BCMATH_NFACT_ALIMIT || (n < (int) BCMath::nCacheFact && (n - k) < 10) )
return LogNoverK(n, k);
// calculate final log(n over k) using approximations if necessary
return ApproxLogFact((double)n) - ApproxLogFact((double)k) - ApproxLogFact((double)(n - k));
}
| double BCMath::LogBreitWignerNonRel | ( | double | x, | |
| double | mean, | |||
| double | Gamma, | |||
| bool | norm = false | |||
| ) |
Calculates the logarithm of the non-relativistic Breit-Wigner distribution.
Definition at line 242 of file BCMath.cxx.
{
double bw = log(Gamma) - log((x - mean) * (x - mean) + Gamma*Gamma / 4.);
if (norm)
bw -= log(2. * M_PI);
return bw;
}
| double BCMath::LogBreitWignerRel | ( | double | x, | |
| double | mean, | |||
| double | Gamma | |||
| ) |
Definition at line 254 of file BCMath.cxx.
{
return -log((x*x - mean*mean) * (x*x - mean*mean) + mean*mean * Gamma*Gamma);
}
| double BCMath::LogChi2 | ( | double | x, | |
| int | n | |||
| ) |
Calculates the logarithm of chi square function: chi2(double x; size_t n)
Definition at line 261 of file BCMath.cxx.
{
if (x < 0) {
BCLog::OutWarning("BCMath::LogChi2 : parameter cannot be negative!");
return -1e99;
}
if (x == 0 && n == 1) {
BCLog::OutWarning("BCMath::LogChi2 : returned value is infinity!");
return 1e99;
}
double nOver2 = ((double) n) / 2.;
return (nOver2 - 1.) * log(x) - x / 2. - nOver2 * log(2) - log(TMath::Gamma(nOver2));
}
| double BCMath::LogFact | ( | int | n | ) |
Calculates natural logarithm of the n-factorial (n!)
Definition at line 132 of file BCMath.cxx.
{
// return NaN for negative argument
if (n<0)
return std::numeric_limits<double>::quiet_NaN();
// cache the factorials on first call
if ( !BCMath::logfact && BCMath::nCacheFact>=0 ) {
BCMath::logfact = new double[BCMath::nCacheFact+1];
double tmplogfact = 0;
BCMath::logfact[0] = tmplogfact;
for (unsigned int i=1; i<=BCMath::nCacheFact; i++) {
tmplogfact += log((double) i);
BCMath::logfact[i] = tmplogfact;
}
}
// return cached value if available
if (n <= (int) BCMath::nCacheFact)
return BCMath::logfact[n];
// calculate factorial starting from the highest cached value
double ln(0.);
if (BCMath::logfact)
ln = BCMath::logfact[nCacheFact];
for (int i = BCMath::nCacheFact+1; i <= n; i++)
ln += log((double) i);
return ln;
}
| double BCMath::LogGaus | ( | double | x, | |
| double | mean = 0, |
|||
| double | sigma = 1, |
|||
| bool | norm = false | |||
| ) |
Calculate the natural logarithm of a gaussian function with mean and sigma. If norm=true (default is false) the result is multiplied by the normalization constant, i.e. divided by sqrt(2*Pi)*sigma.
Definition at line 30 of file BCMath.cxx.
{
// if we have a delta function, return fixed value
if (sigma == 0.)
return 0;
// if sigma is negative use absolute value
if (sigma < 0.)
sigma *= -1.;
double arg = (x - mean) / sigma;
double result = -.5 * arg * arg;
// check if we should add the normalization constant
if (!norm)
return result;
// subtract the log of the denominator of the normalization constant
// and return
return result - log(sqrt(2. * M_PI) * sigma);
}
| double BCMath::LogNoverK | ( | int | n, | |
| int | k | |||
| ) |
Calculates natural logarithm of the Binomial factor "n over k".
Definition at line 173 of file BCMath.cxx.
{
// switch parameters if n < k
if (n < k) {
int a = n;
n = k;
k = a;
}
if (n == k || k == 0)
return 0.;
if (k == 1 || k == n - 1)
return log((double) n);
int lmax = Max(k, n - k);
int lmin = Min(k, n - k);
double ln = 0.;
for (int i = n; i > lmax; i--)
ln += log((double) i);
ln -= LogFact(lmin);
return ln;
}
| double BCMath::LogPoisson | ( | double | x, | |
| double | par | |||
| ) |
Calculate the natural logarithm of a poisson distribution.
Definition at line 54 of file BCMath.cxx.
{
if (par > 899)
return LogGaus(x, par, sqrt(par), true);
if (x < 0)
return 0;
if (x == 0.)
return -par;
return x * log(par) - par - ApproxLogFact(x);
}
| double BCMath::LogVoigtian | ( | double | x, | |
| double | sigma, | |||
| double | gamma | |||
| ) |
Calculates the logarithm of normalized voigtian function: voigtian(double x, double sigma, double gamma)
voigtian is a convolution of the following two functions: gaussian(x) = 1/(sqrt(2*pi)*sigma) * exp(x*x/(2*sigma*sigma) and lorentz(x) = (1/pi)*(gamma/2) / (x*x + (gamma/2)*(gamma/2))
it is singly peaked at x=0. The width of the peak is decided by sigma and gamma, so they should be positive.
Definition at line 279 of file BCMath.cxx.
{
if (sigma <= 0 || gamma <= 0) {
BCLog::OutWarning("BCMath::LogVoigtian : widths are negative or zero!");
return -1e99;
}
return log(TMath::Voigt(x, sigma, gamma));
}
| double BCMath::longestRunFrequency | ( | unsigned int | longestObserved, | |
| unsigned int | nTrials | |||
| ) |
Find the sampling probability that, given n independent Bernoulli trials with success rate = failure rate = 1/2, the longest run of consecutive successes is greater than the longest observed run. Key idea from Burr, E.J. & Cane, G. Longest Run of Consecutive Observations Having a Specified Attribute. Biometrika 48, 461-465 (1961).
| longestObserved | actual longest run | |
| nTrials | number of independent trials |
| double BCMath::longestRunFrequency | ( | unsigned | longestObserved, | |
| unsigned int | nTrials | |||
| ) |
Definition at line 513 of file BCMath.cxx.
{
// can't observe run that's longer than the whole sequence
if (longestObserved >= nTrials)
return 0.;
// return value
double prob = 0.;
// short cuts
typedef unsigned int uint;
uint Lobs = longestObserved;
uint n = nTrials;
// the result of the inner loop is the cond. P given r successes
double conditionalProb;
// first method: use the gamma function for the factorials: bit slower and more inaccurate
// in fact may return NaN for n >= 1000
// double Gup, Gdown1, Gdown2, Gdown3;
//
// for (uint r = 0; r <= n; r++) {
// conditionalProb = 0.0;
// for (uint i = 1; ( i <= n-r+1) && (i <= uint(r / double(Lobs + 1)) ); i++) {
//
// Gup = TMath::Gamma(1 - i * (Lobs + 1) + n);
// Gdown1 = TMath::Gamma(1 + i);
// Gdown2 = TMath::Gamma(2 - i + n - r);
// Gdown3 = TMath::Gamma(1 - i * (Lobs + 1) + r);
//
// //consider the sign of contribution
// Gup = i%2 ? Gup : - Gup;
//
// conditionalProb += Gup/(Gdown1 * Gdown2 * Gdown3);
//
//
// }
// prob += (1 + n - r)*conditionalProb;
// }
// alternative using log factorial approximations, is faster and more accurate
double tempLog = 0;
for (uint r = 0; r <= n; r++) {
conditionalProb = 0.0;
for (uint i = 1; (i <= n - r + 1) && (i <= uint(r / double(Lobs + 1))); i++) {
tempLog = ApproxLogFact(n - i * (Lobs + 1)) - ApproxLogFact(i)
- ApproxLogFact(n - r - i + 1)
- ApproxLogFact(r - i * (Lobs + 1));
if (i % 2)
conditionalProb += exp(tempLog);
else
conditionalProb -= exp(tempLog);
// printf("tempLog inside = %.2f",prob);
}
// printf("tempLog outside = %.2f",prob);
prob += (1 + n - r) * conditionalProb;
}
// Bernoulli probability of each permutation
prob *= pow(2., -double(n));
return prob;
}
| std::vector< int > BCMath::longestRuns | ( | const std::vector< bool > & | bitStream | ) |
Find the longest runs of zeros and ones in the bit stream
| bitStream | input sequence of boolean values |
Definition at line 384 of file BCMath.cxx.
{
// initialize counter variables
unsigned int maxRunAbove, maxRunBelow, currRun;
maxRunAbove = 0;
maxRunBelow = 0;
currRun = 1;
// set both entries to zero
std::vector<int> runs(2, 0);
if (bitStream.empty())
return runs;
// flag about kind of the currently considered run
bool aboveRun = bitStream.at(0);
// start at second variable
std::vector<bool>::const_iterator iter = bitStream.begin();
++iter;
while (iter != bitStream.end()) {
// increase counter if run continues
if (*(iter - 1) == *iter)
currRun++;
else {
// compare terminated run to maximum
if (aboveRun)
maxRunAbove = Max(maxRunAbove, currRun);
else
maxRunBelow = Max(maxRunBelow, currRun);
// set flag to run of opposite kind
aboveRun = !aboveRun;
// restart at length one
currRun = 1;
}
// move to next bit
++iter;
}
// check last run
if (aboveRun)
maxRunAbove = Max(maxRunAbove, currRun);
else
maxRunBelow = Max(maxRunBelow, currRun);
// save the longest runs
runs.at(0) = maxRunBelow;
runs.at(1) = maxRunAbove;
return runs;
}
| std::vector< double > BCMath::longestRunsChi2 | ( | const std::vector< double > & | yMeasured, | |
| const std::vector< double > & | yExpected, | |||
| const std::vector< double > & | sigma | |||
| ) |
Find the longest success/failure run in set of norm. distributed variables. Success = observation >= expectation. Runs are weighted by the total chi^2 of all elements in the run
| yMeasured | the observations | |
| yExpected | the expected values | |
| sigma | the theoretical uncertainties on the expectations |
Definition at line 437 of file BCMath.cxx.
{
//initialize counter variables
double maxRunAbove, maxRunBelow, currRun;
maxRunAbove = 0;
maxRunBelow = 0;
currRun = 0;
//set both entries to zero
std::vector<double> runs(2, 0);
//check input size
if (yMeasured.size() != yExpected.size() || yMeasured.size() != sigma.size()
|| yExpected.size() != sigma.size()) {
//should throw exception
return runs;
}
//exclude zero uncertainty
//...
int N = yMeasured.size();
if ( N<=0)
return runs;
//BCLog::OutDebug(Form("N = %d", N));
//flag about kind of the currently considered run
double residue = (yMeasured.at(0) - yExpected.at(0)) / sigma.at(0);
bool aboveRun = residue >= 0 ? true : false;
currRun = residue * residue;
//start at second variable
for (int i = 1; i < N; i++) {
residue = (yMeasured.at(i) - yExpected.at(i)) / sigma.at(i);
//run continues
if ((residue >= 0) == aboveRun) {
currRun += residue * residue;
} else {
//compare terminated run to maximum
if (aboveRun)
maxRunAbove = Max(maxRunAbove, currRun);
else
maxRunBelow = Max(maxRunBelow, currRun);
//set flag to run of opposite kind
aboveRun = !aboveRun;
//restart at current residual
currRun = residue * residue;
}
//BCLog::OutDebug(Form("maxRunBelow = %g", maxRunBelow));
//BCLog::OutDebug(Form("maxRunAbove = %g", maxRunAbove));
//BCLog::OutDebug(Form("currRun = %g", currRun));
}
//BCLog::OutDebug(Form("maxRunBelow = %g", maxRunBelow));
//BCLog::OutDebug(Form("maxRunAbove = %g", maxRunAbove));
//BCLog::OutDebug(Form("currRun = %g", currRun));
//check last run
if (aboveRun)
maxRunAbove = Max(maxRunAbove, currRun);
else
maxRunBelow = Max(maxRunBelow, currRun);
//BCLog::OutDebug(Form("maxRunBelow = %g", maxRunBelow));
//BCLog::OutDebug(Form("maxRunAbove = %g", maxRunAbove));
//save the longest runs
runs.at(0) = maxRunBelow;
runs.at(1) = maxRunAbove;
return runs;
}
| double BCMath::Max | ( | double | x, | |
| double | y | |||
| ) | [inline] |
| int BCMath::Max | ( | int | x, | |
| int | y | |||
| ) | [inline] |
| int BCMath::Max | ( | unsigned int | x, | |
| unsigned int | y | |||
| ) | [inline] |
| double BCMath::Max | ( | float | x, | |
| float | y | |||
| ) | [inline] |
| double BCMath::Min | ( | float | x, | |
| float | y | |||
| ) | [inline] |
| int BCMath::Min | ( | int | x, | |
| int | y | |||
| ) | [inline] |
| double BCMath::Min | ( | double | x, | |
| double | y | |||
| ) | [inline] |
| int BCMath::Min | ( | unsigned int | x, | |
| unsigned int | y | |||
| ) | [inline] |
| int BCMath::Nint | ( | double | x | ) |
Returns the nearest integer of a double number.
Definition at line 201 of file BCMath.cxx.
{
// round to integer
int i;
if (x >= 0) {
i = (int) (x + .5);
if (x + .5 == (double) i && i&1)
i--;
}
else {
i = int(x - 0.5);
if (x - 0.5 == double(i) && i&1)
i++;
}
return i;
}
| void BCMath::RandomChi2 | ( | std::vector< double > & | randoms, | |
| int | K | |||
| ) |
Get N random numbers distributed according to chi square function with K degrees of freedom
Definition at line 312 of file BCMath.cxx.
{
// fixed upper cutoff to 1000, might be too small
TF1 *f = new TF1("chi2", chi2, 0.0, 1000, 1);
f->SetParameter(0, K);
f->SetNpx(500);
// uses inverse-transform method
// fortunately CDF only built once
for (unsigned int i = 0; i < randoms.size(); i++)
randoms.at(i) = f->GetRandom();
delete f;
}
| double BCMath::rms | ( | int | n, | |
| const double * | a | |||
| ) |
Returns the rms of an array.
Definition at line 222 of file BCMath.cxx.
{
if (n <= 0 || !a)
return 0;
double sum = 0., sum2 = 0.;
for (int i = 0; i < n; i++) {
sum += a[i];
sum2 += a[i] * a[i];
}
double n1 = 1. / (double) n;
double mean = sum * n1;
return sqrt(fabs(sum2 * n1 - mean * mean));
}
| double BCMath::SplitGaussian | ( | double * | x, | |
| double * | par | |||
| ) |
Definition at line 290 of file BCMath.cxx.
{
double mean = par[0];
double sigmadown = par[1];
double sigmaup = par[2];
double sigma = sigmadown;
if (x[0] > mean)
sigma = sigmaup;
return 1.0/sqrt(2.0*M_PI)/sigma * exp(- (x[0]-mean)*(x[0]-mean)/2./sigma/sigma);
}
double* BCMath::logfact = 0 [static] |
Definition at line 26 of file BCMath.cxx.
unsigned int BCMath::nCacheFact = 10000 [static] |
Definition at line 25 of file BCMath.cxx.
1.7.1