Some useful mathematic functions. More...
Functions | |
| double | ApproxBinomial (int n, int k, double p) |
| double | ApproxLogFact (double x) |
| TH1D * | ECDF (const std::vector< double > &data) |
| double | LogApproxBinomial (int n, int k, double p) |
| double | LogBinomFactor (int n, int k) |
| 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 | LogFact (int n) |
| double | LogGaus (double x, double mean=0, double sigma=1, bool norm=false) |
| double | LogNoverK (int n, int k) |
| double | LogPoisson (double x, double par) |
| double | LogVoigtian (double x, double sigma, double gamma) |
| double | longestRunFrequency (unsigned int longestObserved, unsigned int nTrials) |
| 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) |
| int | Max (int x, int y) |
| double | Max (double x, double y) |
| int | Min (int x, int y) |
| double | Min (double x, double y) |
| int | Nint (double x) |
| void | RandomChi2 (std::vector< double > &randoms, int K) |
| double | rms (int n, const double *a) |
| double | SplitGaussian (double *x, double *par) |
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 65 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 117 of file BCMath.cxx.
{
if (x > BCMATH_NFACT_ALIMIT)
return x * log(x) - x + log(x * (1. + 4. * x * (1. + 2. * x))) / 6. + log(M_PI) / 2.;
else
return LogFact((int) x);
}
| 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 293 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 72 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 93 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 - k) < 5)
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 208 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 220 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 227 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 127 of file BCMath.cxx.
{
double ln = 0.;
for (int i = 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 25 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 139 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 49 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 245 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 |
| 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 350 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 = TMath::Max(maxRunAbove, currRun);
else
maxRunBelow = TMath::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 = TMath::Max(maxRunAbove, currRun);
else
maxRunBelow = TMath::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 403 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 = TMath::Max(maxRunAbove, currRun);
else
maxRunBelow = TMath::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 = TMath::Max(maxRunAbove, currRun);
else
maxRunBelow = TMath::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;
}
| int BCMath::Max | ( | int | x, | |
| int | y | |||
| ) | [inline] |
| double BCMath::Max | ( | double | x, | |
| double | y | |||
| ) | [inline] |
| int BCMath::Min | ( | int | x, | |
| int | y | |||
| ) | [inline] |
| double BCMath::Min | ( | double | x, | |
| double | y | |||
| ) | [inline] |
| int BCMath::Nint | ( | double | x | ) |
Returns the nearest integer of a double number.
Definition at line 167 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 278 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 188 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 256 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*TMath::Pi())/sigma * exp(- (x[0]-mean)*(x[0]-mean)/2./sigma/sigma);
}
1.7.1