File: mpchitest.pro

package info (click to toggle)
mpfit 1.84%2B2016.05.19-1
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 784 kB
  • ctags: 45
  • sloc: python: 126; makefile: 16; sh: 13
file content (314 lines) | stat: -rw-r--r-- 8,232 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
;+
; 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