Re: Confidence level calculation

From: Rene Brun (Rene.Brun@cern.ch)
Date: Thu Nov 12 1998 - 16:29:10 MET


Hi Nick,
I have added a few new functions in the class TMath, including
  TMath::Probability
  TMath::Landau
where Probability is based on your function below.
Thanks for submitting the code.

Rene Brun


Nick van Eijndhoven wrote:
> 
> Dear fellow ROOTers,
> In performing some physics analysis I needed to compute the
> confidence level based on a chi-squared and number of statistical
> constraints (i.e. the functionality of the old cernlib routine prob).
> Since I couldn't find such a facility in the ROOT framework (in case
> there is already one, please let me know), I decided to dig into
> my old private fortran utilities.
> Below you find attached the function "conlev" which is a C++
> translation of the old fortran code (the original code is also
> included as comment) which does the job (I tested the reproduction
> of the CL values as noted in the particle data book.
> Note : the result for 1 constraint is incorrect, but normally one
> doesn't encounter this situation).
> I know it won't win a prize in a contest for nice programming,
> but at least it works.
> Feel free to use, modifiy etc.... the code and in case the ROOT
> team wants to include it in the standard ROOT distribution, just
> go ahead with it.
> Curently the easiest way I found to have it always available within
> ROOT is to just include the C++ code in the file rootalias.C.
> --
> 
>                                               Cheers,
> 
>                                _/_/      _/    _/   _/_/_/_/    _/   _/
>                               _/  _/    _/    _/   _/          _/  _/
>                              _/    _/  _/    _/   _/          _/_/
>                             _/      _/_/    _/   _/          _/  _/
>                            _/        _/    _/   _/_/_/_/    _/    _/
> 
> *----------------------------------------------------------------------*
>  Dr. Nick van Eijndhoven                Department of Subatomic Physics
>  email : nick@phys.uu.nl                Utrecht University / NIKHEF
>  tel. +31-30-2532331 (direct)           P.O. Box 80.000
>  tel. +31-30-2531492 (secr.)            NL-3508 TA Utrecht
>  fax. +31-30-2518689                    The Netherlands
>  WWW : http://www.phys.uu.nl/~nick      Office : Ornstein lab. 172
>  ----------------------------------------------------------------------
>  tel. +41-22-7679751 (direct)           CERN PPE Division / ALICE exp.
>  tel. +41-22-7675857 (secr.)            CH-1211 Geneva 23
>  fax. +41-22-7679480                    Switzerland
>  CERN beep : 13+7294                    Office : B 160 1-012
> *----------------------------------------------------------------------*
> 
>     ---------------------------------------------------------------
> ///////////////////////////////////////////////////////////////////////////////
> // Computation of the confidence level from Chi-squared (chi2)
> // and number of constraints (ncons).
> //
> // Note :
> // for even ncons ==> same result as the Cernlib function PROB
> // for odd  ncons ==> confidence level < result of the Cernlib function PROB
> //
> //--- Original Fortan routine copied from Lau Gatignon 1980
> //--- Modified by NvE 29-sep-1980 K.U. Nijmegen
> //--- Modified by NvE 24-apr-1985 NIKHEF-H Amsterdam
> //--- Converted to C++ by Nve 06-nov-1998 UU-SAP Utrecht
> ///////////////////////////////////////////////////////////////////////////////
> float conlev(float chi2,int ncons)
> {
>  const float a1=0.705230784e-1, a2=0.422820123e-1,
>              a3=0.92705272e-2,  a4=0.1520143e-3,
>              a5=0.2765672e-3,   a6=0.430638e-4;
> 
>  // Change to dummy variables
>  int n=ncons;
>  float y0=chi2;
> 
>  // Set CL to zero in case ncons<=0
>  if (n <= 0) return 0;
> 
>  if (y0 <= 0.)
>  {
>   if (y0 < 0.)
>   {
>     return 0;
>   }
>   else
>   {
>    return 1;
>   }
>  }
> 
>  if (n > 100)
>  {
>   float x=sqrt(2.*y0)-sqrt(float(n+n-1));
>   float t=fabs(x)/1.1412135;
>   float denom=1.0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6)))));
>   float vfun=1.0/denom
>   if (fabs(vfun) < 1.3e-5) vfun=0.; // Prevent underflow
>   vfun=pow(vfun,16);
>   float v=0.5*vfun;
>   if (x < 0.) v=1.-v;
>   if (v < 0.) v=0.;
>   return v;
>  }
>  else
>  {
>   float y1=0.5*y0;
>   int m=n/2;
>   if (2*m == n)
>   {
>    float sum=1.;
>    if (m > 1)
>    {
>     float term=1.;
>     for (int i=2; i<=m; i++)
>     {
>      if (term > 1.e20) break; // Prevent overflow
>      float fi=float(i-1);
>      term=term*y1/fi;
>      sum+=term;
>     }
>    }
>     float rnick=y1;
>     if (rnick >= 1.75e2) rnick=1.75e2; // Prevent underflow
>     float v=sum*exp(-rnick);
>     if (v < 0.) v=0.;
>     return v;
>   }
>   else
>   {
>    float x=sqrt(y0);
>    float t=fabs(x)/1.1412135;
>    float denom=1.0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6)))));
>    float vfun=1.0/denom
>    if (fabs(vfun) < 1.3e-5) vfun=0.; // Prevent underflow
>    vfun=pow(vfun,16);
>    float v=vfun;
>    if (n < 1) return 0;
>    if (n == 1)
>    {
>     if (v < 0.) v=0.;
>     return v;
>    }
>    float value=v;
>    float sum=1.;
>    if (n >= 4)
>    {
>     float term=1.;
>     int k=m-1;
>     for (int j=1; j<=k; j++)
>     {
>      if (term > 1.e20) break; // Prevent overflow
>      float fj=float(j);
>      term=term*y0/(fj+fj+1.);
>      sum+=term;
>     }
>    }
>    if (y1 > 1.75e2) y1=1.75e2; // Prevent underflow
>    float vi=sum*0.797846*sqrt(y0)*exp(-y1);
>    v=vi+value;
>    if (v < 0.) v=0.;
>    return v;
>   }
>  }
> }
> ///////////////////////////////////////////////////////////////////////////////
> //        SUBROUTINE CONLEV (CL,NCONS,CHI2)
> //C *** ROUTINE COPIED FROM LAU GATIGNON ***
> //C *** MODIFIED BY NVE 29-SEPT-1980 K.U. NIJMEGEN ***
> //C --- LAST MOD. NVE 24-APR-1985 NIKHEF-H AMSTERDAM ---
> //C
> //C --- THE COMPUTATION OF THE CONFIDENCE LEVEL (CL) FROM CHI-SQAURED (CHI2) ---
> //C --- AND NUMBER OF CONSTRAINTS (NCONS) ---
> //C NOTE ;
> //C FOR EVEN NCONS ==> SAME RESULT AS THE CERN FUNCTION PROB
> //C FOR ODD  NCONS ==> CL < RESULT OF THE CERN FUNCTION PROB
> //C
> //*      DATA A1 /0.705230784E-1/, A2 /0.422820123E-1/
> //*      DATA A3 /0.92705272 E-2/, A4 /0.1520143  E-3/
> //*      DATA A5 /0.2765672  E-3/, A6 /0.430638   E-4/
> //*C
> //*C --- CHANGE TO DUMMY VARIABLES ---
> //*      N=NCONS
> //*      Y0=CHI2
> //C
> //*C --- SET CL TO ZERO IN CASE OF NCONS .LE. 0 ---
> //*      IF (N .GT. 0) GO TO 1
> //*      V=0.0
> //*      GO TO 130
> //*C
> //* 1    IF (Y0 .GT. 0.0) GO TO 10
> //*      V=0.0
> //*      IF (Y0 .EQ. 0.0) V=1.0
> //*      GO TO 130
> //*C
> //* 10   IF (N .LT. 101) GO TO 30
> //*      X=SQRT(Y0+Y0)-SQRT(FLOAT(N+N-1))
> //*      ASSIGN 20 TO JUMP
> //*      GO TO 140
> //*C
> //* 20   V=0.5*VFUN
> //*      IF (X .LT. 0.0) V=1.0-V
> //*      GO TO 110
> //*C
> //* 30   Y1=0.5*Y0
> //*      M = N/2
> //*      IF (M+M .NE. N) GO TO 60
> //*      SUM=1.0
> //*      IF (M .LE. 1) GO TO 50
> //*      TERM=1.0
> //*      DO 40 I=2,M
> //*C
> //*C *** PREVENT OVERFLOW ***
> //*      IF (TERM .GT. 1.0E20) GO TO 50
> //*      FI=I-1
> //*      TERM=TERM*Y1/FI
> //* 40   SUM=SUM+TERM
> //*C *** PREVENT UNDERFLOW ***
> //* 50   RNICK=Y1
> //*      IF (RNICK .GE. 1.75E2) RNICK=1.75E2
> //*      V=SUM*EXP(-RNICK)
> //*      GO TO 110
> //C
> //* 60   X=SQRT(Y0)
> //*      ASSIGN 70 TO JUMP
> //*      GO TO 140
> //*C
> //* 70   V=VFUN
> //*      IF (N-1) 120,110,80
> //* 80   VALUE=V
> //*      SUM=1.0
> //*      IF (N .LT. 4) GO TO 100
> //*      TERM=1.0
> //*      K=M-1
> //*      DO 90 I=1,K
> //*      IF (TERM .GT. 1.0E20) GO TO 100
> //*      FI=I
> //*      TERM=TERM*Y0/(FI+FI+1.0)
> //* 90   SUM=SUM+TERM
> //*C *** PREVENT UNDERFLOW ***
> //* 100  IF (Y1 .GT. 1.75E2) Y1=1.75E2
> //*      VI=SUM*0.797846*SQRT(Y0)*EXP(-Y1)
> //*      V=VI+VALUE
> //* 110  IF (V) 120,130,130
> //* 120  V=0.0
> //* 130  CONTINUE
> //*      CL=V
> //*      RETURN
> //*C
> //* 140  T=ABS(X)/1.1412135
> //*      DENOM=1.0+T*(A1+T*(A2+T*(A3+T*(A4+T*(A5+T*A6)))))
> //*      VFUN=1.0/DENOM
> //*C
> //*C *** PREVENT UNDERFLOW ***
> //*      IF (ABS(VFUN) .LT. 1.3E-5) VFUN=0.0
> //*      VFUN=VFUN**16
> //*      GO TO JUMP,(20,70)
> //*C
> //*      END
> //////////////////////////////////////////////////////////////////////////



This archive was generated by hypermail 2b29 : Tue Jan 04 2000 - 00:34:39 MET