package jhpro.stat;
/*
This is a proprietary part of DataMelt package.
License: the Creative Commons Attribution-Share Alike License; either version 3.0 of the License,
or any later version.
http://creativecommons.org/licenses/by-sa/3.0/
Permission to use, copy and modify this software and its documentation for
NON-COMMERCIAL purposes is granted,
without fee, provided that an acknowledgment to the author, Dr Sergei Chekanov
appears in all copies and associated documentation or publications.
Redistributions of the source code, or parts of the source codes,
must retain the copyright notice, this list of conditions and the following disclaimer
(all at the top of each source code) and requires the written permission of Sergei Chekanov
Redistribution in binary form of all or parts of a class must reproduce the
copyright notice, this list of conditions and the following disclaimer in the documentation and/or other
materials provided with the distribution and requires the written permission of Sergei Chekanov
Details about license agreement or commercial use are: http://jwork.org/scavis
*/
import jhplot.*;
import jhplot.math.*;
import java.util.ArrayList;
import java.util.Random;
import jhplot.gui.*;
import org.apache.commons.math3.util.FastMath;
/**
* Algorithm to compute 95% CL limits using the Likelihood ratio semi-Bayesian
* method. It takes signal, background and data histograms wrapped in a
* DataSource as input and runs a set of Monte Carlo experiments in order
* to compute the limits.
*
* When excluding new physics models, 95% CL exclusion is usually the criterion
* chosen, and one only has to compute CLs with enough precision to tell that an
* observed outcome is less probable than about 5% of the time assuming a signal
* is present. But forming discovery p-values, we must compute 1-CLb values of
* the order of 1 10^{-7}. This computation involves generating of the order of 1
* 10^{8} pseudoexperiments, just to be on the safe side.
*
* If needed, inputs are fluctuated according to systematics.
*
* Reference: HEP-EX/9902006. see: Tom Junk,NIM A434, p. 435-443, 1999 see
* http://root.cern.ch/root/doc/TomJunk.pdf
*
* @author S.Chekanov
*
*/
public class CLimits {
private int nmc;
private boolean stat;
private DataSource data;
private double[] fgTable;
private Random generator;
/**
* Main constructor. Takes the signal, background and data histograms as
* well as different systematics sources to form the limit input. For this,
* nmc Monte Carlo experiments are performed. As usual, the larger this
* number, the longer the compute time, but the better the result.
* Statistical errors are included in the calculations.
*
* @param data
* input data (signal, background, data)
* @param nmc
* Number of MC experiments. Usually 50000.
*/
public CLimits(DataSource data, int nmc) {
this.data = data;
this.nmc = nmc;
this.stat=true;
}
/**
* Main constructor. Takes the signal, background and data histograms as
* well as different systematics sources to form the limit input. For this,
* nmc Monte Carlo experiments are performed. As usual, the larger this
* number, the longer the compute time, but the better the result.
*
* @param data
* input data (signal, background, data)
* @param nmc
* Number of MC experiments. Usually 50000.
* @param stat
* set true if statistical errors are included in the estimate
*/
public CLimits(DataSource data, int nmc, boolean stat) {
this.data = data;
this.nmc = nmc;
this.stat=stat;
}
/**
* Get limits including statistical errors in the calculations.
* Statistical errors are included.
* @param s
* signal
* @param b
* background
* @param d
* data
* @param nmc
* number of MC points
*/
public CLimits(H1D s, H1D b, H1D d, int nmc) {
this(new DataSource(s, b, d), nmc, true);
}
/**
* Initialize limit calculations including statistical errors in the calculations.
*
* @param s
* signal
* @param b
* background
* @param d
* data
* @param nmc
* number of MC points
*/
public CLimits(P1D s, P1D b, P1D d, int nmc) {
this(new DataSource(s, b, d), nmc, true);
}
/**
* Initialize limit calculations for one channel with systematic errors.
* Statistical errors are included.
* One can define systematic error. Example, the number 0.05
* means scaling by 1+0.05 (5% systematic).
*
* @param s
* Signal histogram
* @param b
* Background histogram
* @param d
* Data histogram
* @param se
* Systematic error on signal
* @param be
* Systematic error on background
* @param de
* Systematic errors on data
* @param l
* Name of this channel
* @param nmc
* Number of MC histograms
*/
public CLimits(H1D s, H1D b, H1D d, double se,
double be, double de, String l, int nmc) {
this(new DataSource(s, b, d, se, be, de, l), nmc, true);
}
/**
* Build a Confidence level.
*
* @return
*/
public ConfidenceLevel getLimit() {
return getLimit(new Random());
}
/**
* Build a Confidence level with a custom Random() to build Poisson random
* numbers.
*
* @param stat
* usually false
* @param generator
* . Usually Random(). Used to build Poisson random numbers.
* @return
*/
public ConfidenceLevel getLimit(Random generator) {
this.generator = generator;
if (generator == null)
this.generator = new Random();
Poisson myrandom = new Poisson(this.generator, 0.0);
ConfidenceLevel result = new ConfidenceLevel(nmc);
int maxbins = 0;
double nsig = 0;
double nbg = 0;
int ncand = 0;
int i;
for (i = 0; i < data.getSignal().size(); i++) {
maxbins = ((data.getSignal().get(i)).getSize()) > maxbins ? ((data
.getSignal().get(i)).getSize()) : maxbins;
nsig += (data.getSignal().get(i)).getIntegral();
nbg += (data.getBackground().get(i)).getIntegral();
ncand += (int) (data.getCandidates().get(i)).getIntegral();
}
result.setBtot(nbg);
result.setStot(nsig);
result.setDtot(ncand);
double buffer = 0;
fgTable = new double[maxbins * (data.getSignal().size() + 1)];
for (int channel = 0; channel < data.getSignal().size(); channel++)
for (int bin = 0; bin < (data.getSignal().get(channel)).getSize(); bin++) {
double s = data.getSignal().get(channel).getValue(bin);
double b = data.getBackground().get(channel).getValue(bin);
double d = data.getCandidates().get(channel).getValue(bin);
// Compute the value of the "-2lnQ" for the actual data
if ((b == 0) && (s > 0)) {
System.out.println("WARNING: Ignoring bin "
+ Integer.toString(bin) + " which has s="
+ Double.toString(s) + " and b="
+ Double.toString(b));
System.out
.println(" -> Maybe the MC statistic has to be improved...");
}
if ((s > 0) && (b > 0))
buffer += LogLikelihood(s, b, b, d);
// System.out.println("bin="+Integer.toString(bin)+
// " s="+Double.toString(s)+
// " b="+Double.toString(b)+" d="+Double.toString(d) +
// " logLike="+Double.toString(LogLikelihood(s,b,b,d)));
// precompute the log(1+s/b)'s in an array to speed up
// computation
// background-free bins are set to have a maximum t.s. value
// for protection (corresponding to s/b of about 5E8)
if ((s > 0) && (b > 0)) {
fgTable[(channel * maxbins) + bin] = LogLikelihood(s, b, b,
1);
} else if ((s > 0) && (b == 0)) {
fgTable[(channel * maxbins) + bin] = 20;
}
}
result.setTSD(buffer);
// System.out.println("Accumulated Buffer="+Double.toString(buffer));
// accumulate MC experiments. Hold the test statistic function fixed,
// but fluctuate s and b within syst. errors for computing probabilities
// of
// having that outcome. (Alex Read's prescription -- errors are on the
// ensemble, not on the observed test statistic. This technique does not
// split
// outcomes.)
// keep the tstats as sum log(1+s/b). convert to -2lnQ when preparing
// the results
// (reason -- like to keep the < signs right)
double[] tss = new double[nmc];
double[] tsb = new double[nmc];
double[] lrs = new double[nmc];
double[] lrb = new double[nmc];
for (i = 0; i < nmc; i++) {
tss[i] = 0;
tsb[i] = 0;
lrs[i] = 0;
lrb[i] = 0;
// System.out.println(i);
DataSource fluctuated = fluctuate(data, generator,stat);
// make an independent set of fluctuations for reweighting
// pseudoexperiments it turns out using the same fluctuated ones for
// numerator and denominator is biased.
DataSource fluctuated2 = fluctuate(data, generator,stat);
for (int channel = 0; channel < fluctuated.getSignal().size(); channel++) {
for (int bin = 0; bin < (fluctuated.getSignal().get(channel))
.getSize(); bin++) {
if ((double) (fluctuated.getSignal().get(channel))
.getValue(bin) != 0) {
// s+b hypothesis
double rate = (fluctuated.getSignal().get(channel))
.getValue(bin)
+ (fluctuated.getBackground().get(channel))
.getValue(bin);
double rand = myrandom.next(rate);
tss[i] += rand * fgTable[(channel * maxbins) + bin];
double s = fluctuated.getSignal().get(channel)
.getValue(bin);
double s2 = fluctuated2.getSignal().get(channel)
.getValue(bin);
double b = fluctuated.getBackground().get(channel)
.getValue(bin);
double b2 = fluctuated2.getBackground().get(channel)
.getValue(bin);
if ((s > 0) && (b2 > 0))
lrs[i] += LogLikelihood(s, b, b2, rand) - s - b
+ b2;
else if ((s > 0) && (b2 == 0))
lrs[i] += 20 * rand - s;
// backg. hypothesis
rate = fluctuated.getBackground().get(channel)
.getValue(bin);
rand = myrandom.next(rate);
tsb[i] += rand * fgTable[(channel * maxbins) + bin];
if ((s2 > 0) && (b > 0))
lrb[i] += LogLikelihood(s2, b2, b, rand) - s2 - b2
+ b;
else if ((s > 0) && (b == 0))
lrb[i] += 20 * rand - s;
}
}
}
lrs[i] = FastMath.exp(lrs[i] < 710 ? lrs[i] : 710);
lrb[i] = FastMath.exp(lrb[i] < 710 ? lrb[i] : 710);
}
/*
* for (i = 0; i fgSystNames= input.getErrorNames();
ArrayList fgSysErrorSignal= input.getErrorOnSignal();
ArrayList fgSysErrorBackground= input.getErrorOnBackground();
ArrayList fgSysErrorData= input.getErrorOnData();
DataSource output = new DataSource();
// if there are only stat, just fluctuate stats.
for (int channel = 0; channel < input.getSignal().size(); channel++) {
DataKeeper olddata = input.getCandidates().get(channel);
DataKeeper oldsignal = input.getSignal().get(channel);
DataKeeper newsignal = new DataKeeper(oldsignal.getSize());
newsignal.setTitle(oldsignal.getTitle());
DataKeeper newdata = olddata.copy();
double sysErrorSignal=fgSysErrorSignal.get(channel);
double sysErrorBackground=fgSysErrorBackground.get(channel);
double sysErrorData=fgSysErrorData.get(channel);
if (stat) {
for (int i = 0; i < oldsignal.getSize(); i++) {
newsignal.setValue(i,
oldsignal.getValue(i) + oldsignal.getError(i)*generator.nextGaussian());
newsignal.setError(i, oldsignal.getError(i));
}
} else newsignal = oldsignal.copy();
// take care of systematics
if (sysErrorSignal != 0) newsignal.scale(1+sysErrorSignal);
DataKeeper oldbackground = input.getBackground().get(channel);
DataKeeper newbackground = new DataKeeper(oldbackground.getSize());
newbackground.setTitle(oldbackground.getTitle());
if (stat) {
for (int i = 0; i < oldbackground.getSize(); i++) {
newbackground.setValue(i, oldbackground.getValue(i)
+ oldbackground.getError(i) * generator.nextGaussian());
newbackground.setError(i, oldbackground.getError(i));
}
} else newbackground= oldbackground.copy();
// take care of systematics on background
if (sysErrorBackground != 0) newbackground.scale(1+sysErrorBackground);
if (sysErrorData != 0) newdata.scale(1+sysErrorData);
output.addChannel(newsignal, newbackground, newdata);
}
return output;
}
/**
* Show online documentation.
*/
public void doc() {
String a = this.getClass().getName();
a = a.replace(".", "/") + ".html";
new HelpBrowser(HelpBrowser.JHPLOT_HTTP + a);
}
/**
* Compute LogLikelihood: d*FastMath.log((s+b)/b2);
*
* @param s
* signal number s
* @param b
* background number
* @param b2
* @param d
* @return log of likelihood function
*/
double LogLikelihood(double s, double b, double b2, double d) {
return d * FastMath.log((s + b) / b2);
}
}
Ads help maintain this website.