;+ ; NAME: ; MPNORMTEST ; ; 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 normally distributed Z value ; ; MAJOR TOPICS: ; Curve and Surface Fitting, Statistics ; ; CALLING SEQUENCE: ; PROB = MPNORMTEST(Z, [/CLEVEL, /SLEVEL ]) ; ; DESCRIPTION: ; ; The function MPNORMTEST() computes the probability for the ; magnitude of a value drawn from the normal distribution to equal or ; exceed the given value Z. This can be used for confidence testing ; of a measured value obeying the normal distribution. ; ; P_NORM(ABS(X) > Z) = PROB ; ; In specifying the returned probability level the user has two ; 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). ; ; Note that /SLEVEL and /CLEVEL are mutually exclusive. ; ; INPUTS: ; ; Z - the value to best tested. Z should be drawn from a normal ; distribution with zero mean and unit variance. If a given ; quantity Y has mean MU and standard deviation STD, then Z can ; be computed as Z = (Y-MU)/STD. ; ; RETURNS: ; ; Returns a scalar or vector of probabilities, as described above, ; and according to the /SLEVEL and /CLEVEL keywords. ; ; KEYWORD PARAMETERS: ; ; SLEVEL - if set, then PROB describes the significance level ; (default). ; ; CLEVEL - if set, then PROB describes the confidence level. ; ; EXAMPLES: ; ; print, mpnormtest(5d, /slevel) ; ; Print the probability for the magnitude of a randomly distributed ; variable with zero mean and unit variance to exceed 5, as a ; significance 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 ; Corrected error in handling of CLEVEL keyword, 05 Sep 2003 ; 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 ; Usage message with /CONTINUE, 23 Sep 2009, CM ; ; $Id: mpnormtest.pro,v 1.9 2009/09/23 20:12:46 craigm Exp $ ;- ; Copyright (C) 1997-2001, 2003, 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_polevl, cephes_erfc, cephes_erf, mpnormtest ;; 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 pro cephes_set_erf_common COMPILE_OPT strictarr common cephes_erf_data, pp, qq, rr, ss, tt, uu, uthresh pp = [ 2.46196981473530512524D-10, 5.64189564831068821977D-1, $ 7.46321056442269912687D0, 4.86371970985681366614D1, $ 1.96520832956077098242D2, 5.26445194995477358631D2, $ 9.34528527171957607540D2, 1.02755188689515710272D3, $ 5.57535335369399327526D2 ] qq = [ 1.00000000000000000000D0, 1.32281951154744992508D1, $ 8.67072140885989742329D1, 3.54937778887819891062D2, $ 9.75708501743205489753D2, 1.82390916687909736289D3, $ 2.24633760818710981792D3, 1.65666309194161350182D3, $ 5.57535340817727675546D2 ] rr = [ 5.64189583547755073984D-1, 1.27536670759978104416D0, $ 5.01905042251180477414D0, 6.16021097993053585195D0, $ 7.40974269950448939160D0, 2.97886665372100240670D0 ] ss = [ 1.00000000000000000000D0, 2.26052863220117276590D0, $ 9.39603524938001434673D0, 1.20489539808096656605D1, $ 1.70814450747565897222D1, 9.60896809063285878198D0, $ 3.36907645100081516050D0 ] tt = [ 9.60497373987051638749D0, 9.00260197203842689217D1, $ 2.23200534594684319226D3, 7.00332514112805075473D3, $ 5.55923013010394962768D4 ] uu = [ 1.00000000000000000000D0, 3.35617141647503099647D1, $ 5.21357949780152679795D2, 4.59432382970980127987D3, $ 2.26290000613890934246D4, 4.92673942608635921086D4 ] uthresh = 37.519379347D return end ; erfc.c ; ; Complementary error function ; ; ; ; SYNOPSIS: ; ; double x, y, erfc(); ; ; y = erfc( x ); ; ; ; ; DESCRIPTION: ; ; ; 1 - erf(x) = ; ; inf. ; - ; 2 | | 2 ; erfc(x) = -------- | exp( - t ) dt ; sqrt(pi) | | ; - ; x ; ; ; For small x, erfc(x) = 1 - erf(x); otherwise rational ; approximations are computed. ; ; ; ; ACCURACY: ; ; Relative error: ; arithmetic domain # trials peak rms ; DEC 0, 9.2319 12000 5.1e-16 1.2e-16 ; IEEE 0,26.6417 30000 5.7e-14 1.5e-14 ; ; ; ERROR MESSAGES: ; ; message condition value returned ; erfc underflow x > 9.231948545 (DEC) 0.0 ; ; ; / function cephes_erfc, a COMPILE_OPT strictarr common cephes_erf_data if n_elements(p) EQ 0 then cephes_set_erf_common common cephes_machar, machvals MAXLOG = machvals.maxlog if a LT 0 then x = -a else x = a if x LT 1. then return, 1.D - cephes_erf(a) z = -a * a if z LT -MAXLOG then begin under: ; message, 'ERROR: underflow', /info if a LT 0 then return, 2.D else return, 0.D endif z = exp(z) if x LT 8. then begin p = cephes_polevl(x, pp) q = cephes_polevl(x, qq) endif else begin p = cephes_polevl(x, rr) q = cephes_polevl(x, ss) endelse y = (z*p)/q if a LT 0 then y = 2.D - y if y EQ 0 then goto, under return, y end ; erf.c ; ; Error function ; ; ; ; SYNOPSIS: ; ; double x, y, erf(); ; ; y = erf( x ); ; ; ; ; DESCRIPTION: ; ; The integral is ; ; x ; - ; 2 | | 2 ; erf(x) = -------- | exp( - t ) dt. ; sqrt(pi) | | ; - ; 0 ; ; The magnitude of x is limited to 9.231948545 for DEC ; arithmetic; 1 or -1 is returned outside this range. ; ; For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise ; erf(x) = 1 - erfc(x). ; ; ; ; ACCURACY: ; ; Relative error: ; arithmetic domain # trials peak rms ; DEC 0,1 14000 4.7e-17 1.5e-17 ; IEEE 0,1 30000 3.7e-16 1.0e-16 ; ; function cephes_erf, x COMPILE_OPT strictarr common cephes_erf_data if abs(x) GT 1. then return, 1.D - cephes_erfc(x) if n_elements(p) EQ 0 then cephes_set_erf_common z = x * x y = x * cephes_polevl(z, tt) / cephes_polevl(z, uu) return, y end function mpnormtest, a, clevel=clevel, slevel=slevel COMPILE_OPT strictarr if n_params() EQ 0 then begin message, 'USAGE: PROB = MPNORMTEST(Z, [/CLEVEL, /SLEVEL ])', /cont return, !values.d_nan endif cephes_setmachar ;; Set machine constants y = a*0 sqrth = sqrt(2.D)/2.D x = a * sqrth ;; Default is to return the significance level if n_elements(slevel) EQ 0 AND n_elements(clevel) EQ 0 then slevel = 1 if keyword_set(slevel) then begin for i = 0L, n_elements(y)-1 do begin if abs(x[i]) LT sqrth then y[i] = 1.D - cephes_erf(abs(x[i])) $ else y[i] = cephes_erfc(abs(x[i])) endfor endif else if keyword_set(clevel) then begin for i = 0L, n_elements(y)-1 do begin if abs(x[i]) LT sqrth then y[i] = cephes_erf(abs(x[i])) $ else y[i] = 1.D - cephes_erfc(x[i]) endfor endif else begin message, 'ERROR: must specify one of CLEVEL, SLEVEL' endelse return, y end