;+
; NAME:
;   MPCHITEST
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Compute the probability of a given chi-squared value
;
; MAJOR TOPICS:
;   Curve and Surface Fitting, Statistics
;
; CALLING SEQUENCE:
;   PROB = MPCHITEST(CHI, DOF, [/SIGMA, /CLEVEL, /SLEVEL ])
;
; DESCRIPTION:
;
;  The function MPCHITEST() computes the probability for a value drawn
;  from the chi-square distribution to equal or exceed the given value
;  CHI.  This can be used for confidence testing of a measured value
;  obeying the chi-squared distribution.
;
;    P_CHI(X > CHI; DOF) = PROB
;
;  In specifying the returned probability level the user has three
;  choices:
;
;    * return the confidence level when the /CLEVEL keyword is passed;
;      OR
;
;    * return the significance level (i.e., 1 - confidence level) when
;      the /SLEVEL keyword is passed (default); OR
;
;    * return the "sigma" of the probability (i.e., compute the
;      probability based on the normal distribution) when the /SIGMA
;      keyword is passed.
;
;  Note that /SLEVEL, /CLEVEL and /SIGMA are mutually exclusive.
;
; INPUTS:
;
;   CHI - chi-squared value to be tested.
;
;   DOF - scalar or vector number, giving the number of degrees of
;         freedom in the chi-square distribution.
;
; RETURNS:
;
;  Returns a scalar or vector of probabilities, as described above,
;  and according to the /SLEVEL, /CLEVEL and /SIGMA keywords.
;
; KEYWORD PARAMETERS:
;
;   SLEVEL - if set, then PROB describes the significance level
;            (default).
;
;   CLEVEL - if set, then PROB describes the confidence level.
;
;   SIGMA - if set, then PROB is the number of "sigma" away from the
;           mean in the normal distribution.
;
; EXAMPLES:
;
;   print, mpchitest(1300d,1252d)
;
;   Print the probability for a chi-squared value with 1252 degrees of
;   freedom to exceed a value of 1300, as a confidence level.
;
; REFERENCES:
;
;   Algorithms taken from CEPHES special function library, by Stephen
;   Moshier. (http://www.netlib.org/cephes/)
;
; MODIFICATION HISTORY:
;   Completed, 1999, CM
;   Documented, 16 Nov 2001, CM
;   Reduced obtrusiveness of common block and math error handling, 18
;     Nov 2001, CM
;   Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
;   Move STRICTARR compile option inside each function/procedure, 9 Oct 2006
;   Add usage message, 24 Nov 2006, CM
;   Really add usage message, with /CONTINUE, 23 Sep 2009, CM
;
;  $Id: mpchitest.pro,v 1.10 2009/10/05 16:22:44 craigm Exp $
;-
; Copyright (C) 1997-2001, 2006, 2009, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

forward_function cephes_igamc, cephes_igam

;; Set machine constants, once for this session.  Double precision
;; only.
pro cephes_setmachar
  COMPILE_OPT strictarr
  common cephes_machar, cephes_machar_vals
  if n_elements(cephes_machar_vals) GT 0 then return

  if (!version.release) LT 5 then dummy = check_math(1, 1)

  mch = machar(/double)
  machep = mch.eps
  maxnum = mch.xmax
  minnum = mch.xmin
  maxlog = alog(mch.xmax)
  minlog = alog(mch.xmin)
  maxgam = 171.624376956302725D

  cephes_machar_vals = {machep: machep, maxnum: maxnum, minnum: minnum, $
                        maxlog: maxlog, minlog: minlog, maxgam: maxgam}

  if (!version.release) LT 5 then dummy = check_math(0, 0)
  return
end

function cephes_igam, a, x
;   
;   	Incomplete gamma integral
;   
;   
;   
;    SYNOPSIS:
;   
;    double a, x, y, igam();
;   
;    y = igam( a, x );
;   
;    DESCRIPTION:
;   
;    The function is defined by
;   
;                              x
;                               -
;                      1       | |  -t  a-1
;     igam(a,x)  =   -----     |   e   t   dt.
;                     -      | |
;                    | (a)    -
;                              0
;   
;   
;    In this implementation both arguments must be positive.
;    The integral is evaluated by either a power series or
;    continued fraction expansion, depending on the relative
;    values of a and x.
;   
;    ACCURACY:
;   
;                         Relative error:
;    arithmetic   domain     # trials      peak         rms
;       IEEE      0,30       200000       3.6e-14     2.9e-15
;       IEEE      0,100      300000       9.9e-14     1.5e-14
  COMPILE_OPT strictarr
  common cephes_machar, machvals
  MAXLOG = machvals.maxlog
  MACHEP = machvals.machep

  if x LE 0 OR a LE 0 then return, 0.D
  if x GT 1. AND x GT a then return, 1.D - cephes_igamc(a, x)
  
  ax = a * alog(x) - x - lngamma(a)
  if ax LT -MAXLOG then begin
;      message, 'WARNING: underflow', /info
      return, 0.D
  endif
  ax = exp(ax)
  r = a
  c = 1.D
  ans = 1.D
  
  repeat begin
      r = r + 1
      c = c * x/r
      ans = ans + c
  endrep until (c/ans LE MACHEP)

  return, ans*ax/a
end

function cephes_igamc, a, x
;   
;   	Complemented incomplete gamma integral
;   
;   
;   
;    SYNOPSIS:
;   
;    double a, x, y, igamc();
;   
;    y = igamc( a, x );
;   
;    DESCRIPTION:
;   
;    The function is defined by
;   
;   
;     igamc(a,x)   =   1 - igam(a,x)
;   
;                               inf.
;                                 -
;                        1       | |  -t  a-1
;                  =   -----     |   e   t   dt.
;                       -      | |
;                      | (a)    -
;                                x
;   
;   
;    In this implementation both arguments must be positive.
;    The integral is evaluated by either a power series or
;    continued fraction expansion, depending on the relative
;    values of a and x.
;   
;    ACCURACY:
;   
;    Tested at random a, x.
;                   a         x                      Relative error:
;    arithmetic   domain   domain     # trials      peak         rms
;       IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
;       IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15

  COMPILE_OPT strictarr

  if n_params() EQ 0 then begin
      message, 'USAGE: PROB = MPCHITEST(CHI, DOF, [/SIGMA, /CLEVEL, /SLEVEL ])', /info
      return, !values.d_nan
  endif


  common cephes_machar, machvals
  MAXLOG = machvals.maxlog
  MACHEP = machvals.machep

  big = 4.503599627370496D15
  biginv = 2.22044604925031308085D-16

  if x LE 0 OR a LE 0 then return, 1.D
  if x LT 1. OR x LT a then return, 1.D - cephes_igam(a, x)
  ax = a * alog(x) - x - lngamma(a)

  if ax LT -MAXLOG then begin
;      message, 'ERROR: underflow', /info
      return, 0.D
  endif
  
  ax = exp(ax)
  y = 1.D - a
  z = x + y + 1.D
  c = 0.D
  pkm2 = 1.D
  qkm2 = x
  pkm1 = x + 1.D
  qkm1 = z * x
  ans = pkm1 / qkm1

  repeat begin
      c = c + 1.D
      y = y + 1.D
      z = z + 2.D
      yc = y * c
      pk = pkm1 * z - pkm2 * yc
      qk = qkm1 * z - qkm2 * yc
      if qk NE 0 then begin
          r = pk/qk
          t = abs( (ans-r)/r )
          ans = r
      endif else begin
          t = 1.D
      endelse
      pkm2 = pkm1
      pkm1 = pk
      qkm2 = qkm1
      qkm1 = qk
      if abs(pk) GT big then begin
          pkm2 = pkm2 * biginv
          pkm1 = pkm1 * biginv
          qkm2 = qkm2 * biginv
          qkm1 = qkm1 * biginv
      endif
  endrep until t LE MACHEP

  return, ans * ax
end

; MPCHITEST
;  compute the probability for a chi-squared value to exceed x give
;  the number of degrees of freedom dof.
function mpchitest, x, dof, slevel=slevel, clevel=clevel, sigma=sigma
  COMPILE_OPT strictarr

  if n_params() LT 2 then begin
      message, 'USAGE: PROB = MPCHITEST(CHI, DOF, [/SIGMA, /CLEVEL, /SLEVEL ])', /cont
      return, !values.d_nan
  endif

  cephes_setmachar   ;; Set machine constants
  
  p = double(x) * 0
  for i = 0, n_elements(x)-1 do begin
      p[i] = cephes_igamc(0.5D * dof, 0.5D * double(x[i]))
  endfor
  if keyword_set(clevel) then return, 1D - double(p)
  if keyword_set(sigma) then return, mpnormlim(p, /slevel)

  return, p
end