|
|
|
00001 #include "DQMServices/Core/src/QStatisticalTests.h" 00002 #include <TMath.h> 00003 00004 using namespace TMath; 00005 00006 //--------------------------------------------------------------------------- 00007 void BinLogLikelihoodRatio(long Nentries, long Nfailures, double epsilon_max, 00008 double* S_fail_obs, double* S_pass_obs ) 00009 { 00010 /*--------------------------------------------------------------------------+ 00011 | Description: Log-likelihood Ratio for Binomial PDF | 00012 +--------------------------------------------------------------------------+ 00013 | Input to this function | 00014 +--------------------------------------------------------------------------+ 00015 |int Nentries : The number of attempts | 00016 |int Nfailures : The number of failures | 00017 |double epsilon_max, : maximum allowed failure rate fraction | 00018 |double* S_fail_obs : uninitialised Significance of failure | 00019 |double* S_pass_obs : uninitialised Significance of Success | 00020 +--------------------------------------------------------------------------+ 00021 | Return values for this function | 00022 +--------------------------------------------------------------------------+ 00023 |double* S_fail_obs : the observed Significance of failure | 00024 |double* S_pass_obs : the observed Significance of Success | 00025 +--------------------------------------------------------------------------+ 00026 | Author: Richard Cavanaugh, University of Florida | 00027 | email: Richard.Cavanaugh@cern.ch | 00028 | Creation Date: 11.July.2005 | 00029 | Last Modified: 17.Jan.2006 | 00030 | Comments: | 00031 +--------------------------------------------------------------------------*/ 00032 long N = Nentries, n = Nfailures; 00033 if( n == 0 ) n = 1; //protect against no failures: approx. 0 by 1 00034 if( n == N ) n -= 1; //protect against all failures: approx. n by (n - 1) 00035 double epsilon_meas = (double) n / (double) N; 00036 00037 double LogQ = 00038 ((double) n ) * ( Log(epsilon_meas) - Log(epsilon_max) ) + 00039 ((double) (N - n)) * ( Log(1.0 - epsilon_meas) - Log(1.0 - epsilon_max) ); 00040 00041 //x-check: var of binomial = epsilon_max * (1 - epsilon_max) / N 00042 if( Nentries <= 1 ) //guard against insufficient entries 00043 { 00044 *S_fail_obs = 0.0; 00045 *S_pass_obs = 0.0; 00046 } 00047 else if( Nfailures == 0 && ( epsilon_max <= 1.0 / (double) Nentries ) ) 00048 { 00049 *S_fail_obs = 0.0; 00050 *S_pass_obs = 0.0; 00051 } 00052 else if( Nfailures == 0 ) 00053 { 00054 *S_fail_obs = 0.0; 00055 *S_pass_obs = sqrt( 2.0 * LogQ ); 00056 } 00057 else if( Nfailures == Nentries ) 00058 { 00059 *S_fail_obs = sqrt( 2.0 * LogQ ); 00060 *S_pass_obs = 0.0; 00061 } 00062 else if( epsilon_meas >= epsilon_max ) 00063 { 00064 *S_fail_obs = sqrt( 2.0 * LogQ ); 00065 *S_pass_obs = 0.0; 00066 } 00067 else 00068 { 00069 *S_fail_obs = 0.0; 00070 *S_pass_obs = sqrt( 2.0 * LogQ ); 00071 } 00072 } 00073 //--------------------------------------------------------------------------- 00074 00075 //--------------------------------------------------------------------------- 00076 void PoissionLogLikelihoodRatio(double data, double hypothesis, 00077 double epsilon_max, double epsilon_min, 00078 double* S_fail_obs, double* S_pass_obs ) 00079 { 00080 /*--------------------------------------------------------------------------+ 00081 | Description: Log-likelihood Ratio for Poission PDF | 00082 +--------------------------------------------------------------------------+ 00083 | Input to this function | 00084 +--------------------------------------------------------------------------+ 00085 |double data, : The observed number of entries | 00086 |double sigma, : The uncertainty on, data, the observed entries | 00087 |double hypothesis, : The assumed hypothese, tested against data | 00088 |double epsilon_max, : Maximum tolerance above fraction of fitted line | 00089 |double epsilon_min, : Minimum tolerance below fraction of fitted line | 00090 |double* S_fail_obs : uninitialised Significance of failure | 00091 |double* S_pass_obs : uninitialised Significance of Success | 00092 +--------------------------------------------------------------------------+ 00093 | Return values for this function | 00094 +--------------------------------------------------------------------------+ 00095 |double* S_fail_obs : the observed Significance of failure | 00096 |double* S_pass_obs : the observed Significance of Success | 00097 +--------------------------------------------------------------------------+ 00098 | Author: Richard Cavanaugh, University of Florida | 00099 | email: Richard.Cavanaugh@cern.ch | 00100 | Creation Date: 14.Jan.2006 | 00101 | Last Modified: 16.Jan.2006 | 00102 | Comments: | 00103 +--------------------------------------------------------------------------*/ 00104 double tolerance_min = hypothesis*(1.0-epsilon_min); 00105 double tolerance_max = hypothesis*(1.0+epsilon_max); 00106 *S_pass_obs = 0.0; 00107 *S_fail_obs = 0.0; 00108 if( data > tolerance_max ) 00109 { 00110 double Nsig = data - tolerance_max; 00111 double Nbak = tolerance_max; 00112 double LogQ = (double) (Nsig + Nbak) * 00113 Log( 1.0 + (double) Nsig / (double) Nbak ) - (double) Nsig; 00114 *S_fail_obs = sqrt( 2.0 * LogQ ); 00115 } 00116 else if( tolerance_min < data && data < tolerance_max ) 00117 { 00118 if( data - hypothesis > 0.0 ) 00119 { 00120 double Nsig = tolerance_max - data; 00121 double Nbak = tolerance_max; 00122 double LogQ = (double) (Nsig + Nbak) * 00123 Log( 1.0 + (double) Nsig / (double) Nbak ) - (double) Nsig; 00124 *S_pass_obs = sqrt( 2.0 * LogQ ); 00125 } 00126 else 00127 { 00128 double Nsig = data - tolerance_min; 00129 double Nbak = tolerance_min; 00130 double LogQ = (double) (Nsig + Nbak) * 00131 Log( 1.0 + (double) Nsig / (double) Nbak ) - (double) Nsig; 00132 *S_pass_obs = sqrt( 2.0 * LogQ ); 00133 } 00134 } 00135 else // data < tolerance_min 00136 { 00137 double Nsig = tolerance_min - data; 00138 double Nbak = tolerance_min; 00139 double LogQ = (double) (Nsig + Nbak) * 00140 Log( 1.0 + (double) Nsig / (double) Nbak ) - (double) Nsig; 00141 *S_fail_obs = sqrt( 2.0 * LogQ ); 00142 } 00143 } 00144 //---------------------------------------------------------------------------
1.5.4