;+ ; NAME: ; MPNORMLIM ; ; 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 confidence limits for normally distributed variable ; ; MAJOR TOPICS: ; Curve and Surface Fitting, Statistics ; ; CALLING SEQUENCE: ; Z = MPNORMLIM(PROB, [/CLEVEL, /SLEVEL ]) ; ; DESCRIPTION: ; ; The function MPNORMLIM() computes confidence limits of a normally ; distributed variable (with zero mean and unit variance), for a ; desired probability level. The returned values, Z, are the ; limiting values: a the magnitude of a normally distributed value ; is greater than Z by chance with a probability PROB: ; ; P_NORM(ABS(X) > Z) = PROB ; ; In specifying the probability level the user has two choices: ; ; * give the confidence level (default); ; ; * give the significance level (i.e., 1 - confidence level) and ; pass the /SLEVEL keyword; OR ; ; Note that /SLEVEL and /CLEVEL are mutually exclusive. ; ; INPUTS: ; ; PROB - scalar or vector number, giving the desired probability ; level as described above. ; ; RETURNS: ; ; Returns a scalar or vector of normal confidence limits. ; ; KEYWORD PARAMETERS: ; ; SLEVEL - if set, then PROB describes the significance level. ; ; CLEVEL - if set, then PROB describes the confidence level ; (default). ; ; EXAMPLE: ; ; print, mpnormlim(0.99d, /clevel) ; ; Print the 99% confidence limit for a normally distributed ; variable. In this case it is about 2.58 sigma. ; ; 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 ; ; $Id: mpnormlim.pro,v 1.6 2006/11/25 01:44:13 craigm Exp $ ;- ; Copyright (C) 1997-2001, 2006, 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_polevl, cephes_ndtri, mpnormlim ;; 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_polevl, x, coef COMPILE_OPT strictarr ans = coef[0] nc = n_elements(coef) for i = 1L, nc-1 do ans = ans * x + coef[i] return, ans end function cephes_ndtri, y0 ; ; Inverse of Normal distribution function ; ; ; ; SYNOPSIS: ; ; double x, y, ndtri(); ; ; x = ndtri( y ); ; ; ; ; DESCRIPTION: ; ; Returns the argument, x, for which the area under the ; Gaussian probability density function (integrated from ; minus infinity to x) is equal to y. ; ; ; For small arguments 0 < y < exp(-2), the program computes ; z = sqrt( -2.0 * log(y) ); then the approximation is ; x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). ; There are two rational functions P/Q, one for 0 < y < exp(-32) ; and the other for y up to exp(-2). For larger arguments, ; w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). ; ; ; ACCURACY: ; ; Relative error: ; arithmetic domain # trials peak rms ; DEC 0.125, 1 5500 9.5e-17 2.1e-17 ; DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17 ; IEEE 0.125, 1 20000 7.2e-16 1.3e-16 ; IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 ; ; ; ERROR MESSAGES: ; ; message condition value returned ; ndtri domain x <= 0 -MAXNUM ; ndtri domain x >= 1 MAXNUM COMPILE_OPT strictarr common cephes_ndtri_data, s2pi, p0, q0, p1, q1, p2, q2 if n_elements(s2pi) EQ 0 then begin s2pi = sqrt(2.D*!dpi) p0 = [ -5.99633501014107895267D1, 9.80010754185999661536D1, $ -5.66762857469070293439D1, 1.39312609387279679503D1, $ -1.23916583867381258016D0 ] q0 = [ 1.D, $ 1.95448858338141759834D0, 4.67627912898881538453D0, $ 8.63602421390890590575D1, -2.25462687854119370527D2, $ 2.00260212380060660359D2, -8.20372256168333339912D1, $ 1.59056225126211695515D1, -1.18331621121330003142D0 ] p1 = [ 4.05544892305962419923D0, 3.15251094599893866154D1, $ 5.71628192246421288162D1, 4.40805073893200834700D1, $ 1.46849561928858024014D1, 2.18663306850790267539D0, $ -1.40256079171354495875D-1,-3.50424626827848203418D-2,$ -8.57456785154685413611D-4 ] q1 = [ 1.D, $ 1.57799883256466749731D1, 4.53907635128879210584D1, $ 4.13172038254672030440D1, 1.50425385692907503408D1, $ 2.50464946208309415979D0, -1.42182922854787788574D-1,$ -3.80806407691578277194D-2,-9.33259480895457427372D-4 ] p2 = [ 3.23774891776946035970D0, 6.91522889068984211695D0, $ 3.93881025292474443415D0, 1.33303460815807542389D0, $ 2.01485389549179081538D-1, 1.23716634817820021358D-2,$ 3.01581553508235416007D-4, 2.65806974686737550832D-6,$ 6.23974539184983293730D-9 ] q2 = [ 1.D, $ 6.02427039364742014255D0, 3.67983563856160859403D0, $ 1.37702099489081330271D0, 2.16236993594496635890D-1,$ 1.34204006088543189037D-2, 3.28014464682127739104D-4,$ 2.89247864745380683936D-6, 6.79019408009981274425D-9] endif common cephes_machar, machvals MAXNUM = machvals.maxnum if y0 LE 0 then begin message, 'ERROR: domain', /info return, -MAXNUM endif if y0 GE 1 then begin message, 'ERROR: domain', /info return, MAXNUM endif code = 1 y = y0 exp2 = exp(-2.D) if y GT (1.D - exp2) then begin y = 1.D - y code = 0 endif if y GT exp2 then begin y = y - 0.5 y2 = y * y x = y + y * y2 * cephes_polevl(y2, p0) / cephes_polevl(y2, q0) x = x * s2pi return, x endif x = sqrt( -2.D * alog(y)) x0 = x - alog(x)/x z = 1.D/x if x LT 8. then $ x1 = z * cephes_polevl(z, p1) / cephes_polevl(z, q1) $ else $ x1 = z * cephes_polevl(z, p2) / cephes_polevl(z, q2) x = x0 - x1 if code NE 0 then x = -x return, x end ; MPNORMLIM - given a probability level, return the corresponding ; "sigma" level. ; ; p - Either the significance level (if SLEVEL is set) or the ; confidence level (if CLEVEL is set). This should be the ; two-tailed level, ie: ; ; * SLEVEL: p = Prob(|z| > z0) ; * CLEVEL: p = Prob(|z| < z0) ; function mpnormlim, p, clevel=clevel, slevel=slevel COMPILE_OPT strictarr if n_params() EQ 0 then begin message, 'USAGE: Z = MPNORMLIM(PROB, [/CLEVEL, /SLEVEL ])', /info return, !values.d_nan endif cephes_setmachar ;; Set machine constants ;; Default is to assume the confidence level if n_elements(clevel) EQ 0 then clevel = 1 y = 0 * p ;; cephes_ndtri accepts the integrated probability from negative ;; infinity to z, so we have to compute. if keyword_set(slevel) then begin p1 = 0.5D * p ;; Take only one of the two tails for i = 0L, n_elements(y)-1 do begin y[i] = - cephes_ndtri(p1[i]) endfor endif else if keyword_set(clevel) then begin p1 = 0.5D + 0.5D * p ;; On binary computers this computation is ;; exact (to the machine precision), so don't worry about it. ;; This computation shaves off the top half of the confidence ;; region, and then adds the "negative infinity to zero part. for i = 0L, n_elements(y)-1 do begin y[i] = cephes_ndtri(p1[i]) endfor endif else begin message, 'ERROR: must specify one of CLEVEL or SLEVEL' endelse return, y end