kuiperone.pro

basic median filter simulation

 pro kuiperone, data, func_name, d, prob, PLOT = plot, _EXTRA = extra
;+
; NAME:
;       KUIPERONE
; PURPOSE:
;       Compute the one-sided Kuiper statistic (invariant Kolmogorov-Smirnov)
; EXPLANATION:
;       Returns the Kuiper statistic and associated probability
;       for an array of data values and a user-supplied cumulative distribution
;       function (CDF) of a single variable.   Algorithm adapted from KSONE
;       in "Numerical Recipes" by Press et al. 2nd edition (1992)
;
;       Kuiper's test is especially useful for data defined on a circle or 
;       to search for periodicity (see Paltani 2004, A&A, 420, 789). 
; CALLING SEQUENCE:
;       kuiperone, data, func_name, D, prob, [ /PLOT ]
;
; INPUT PARAMETERS:
;       data -  vector of data values, must contain at least 4 elements for the
;               Kuiper statistic to be meaningful
;       func_name - scalar string giving the name of the cumulative distribution
;               function.    The function must be defined to accept the data
;               vector as its only input (see example).
;
; OUTPUT PARAMETERS:
;       D - floating scalar giving the Kuiper statistic.   It
;               specifies the sum of positive and negative deviations between the
;               cumulative distribution of the data and the supplied function
;       prob - floating scalar between 0 and 1 giving the significance level of
;               the Kuiper statistic.   Small values of PROB show that the
;               cumulative distribution function of DATA is significantly
;               different from FUNC_NAME.
;
; OPTIONAL INPUT KEYWORD:
;       PLOT - If this keyword is set and non-zero, then KUIPERONE will display a
;               plot of the CDF of the data with the supplied function
;               superposed.   The data values where the Kuiper statistic is
;               computed (i.e. at the maximum difference between the data CDF
;               and the function) are indicated by vertical dashed lines.
;               KUIPERONE accepts the _EXTRA keyword, so that most plot keywords
;               (e.g. TITLE, XTITLE, XSTYLE) can also be passed to KUIPERONE.
;
; EXAMPLE:
;       Determine if a vector created by the RANDOMN function is really
;       consistent with a Gaussian distribution.
;       The CDF of a Gaussian is the error function except that a factor
;       of 2 is included in the error function.   So we must create a special
;       function:
;
;       function gauss_cdf, x
;       return, errorf( x/sqrt(2) )
;       end
;
;       IDL> data = randomn(seed, 50)          ;create data array to be tested
;       IDL> kuiperone, data, 'gauss_pdf', D, prob, /PLOT     ;Use Kuiper test
;
;       A small value of PROB indicates that the cumulative distribution of
;       DATA is significantly different from a Gaussian
;
; NOTES:
;       Note that the 2nd (1992) edition of Numerical Recipes includes
;       a more accurate computation of the K-S significance for small
;       values of N.
;
; PROCEDURE CALLS
;       procedure PROB_KUIPER - computes significance of Kuiper distribution
;
; REVISION HISTORY:
;       Written     W. Landsman                   August, 1992
;       Accept _EXTRA keywords   W. Landsman      September, 1995
;       Fixed possible bug in plot display showing position maximum difference
;       in histogram   M. Fardal/ W. Landsman      March, 1997
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Adapted from KSONE      J. Ballet     July 2003
;-
 On_error, 2

 if ( N_params() LT 3 ) then begin
    print,'Syntax - kuiperone, data, func_name, D, [prob ,/PLOT]'
    return
 endif

 N = N_elements( data )
 if N LT 3 then message, $
   'ERROR - Input data values (first param) must contain at least 3 values'

 sortdata = data[ sort( data ) ]

 f0 = findgen(N)/ N
 fn = ( findgen( N ) +1. ) / N
 ff = call_function( func_name, sortdata )

; Maximum distance above the reference
 D1 = max( fn-ff, subn )

; Maximum distance below the reference
 D2 = max( ff-f0, sub0 )

 D = D1 + D2

 if keyword_set(plot) then begin

; Prepare the step function
     xx = REBIN(sortdata,2*N,/SAMPLE)
     yy = REBIN(f0,2*N,/SAMPLE)
     yy = [yy[1:*],1.]

     plot, xx,yy,_EXTRA = extra
     plots, [sortdata[sub0], sortdata[sub0]], [0,ff[sub0]], linestyle=2
     plots, [sortdata[subn], sortdata[subn]], [ff[subn],1], linestyle=2

; Plot the expected cumulative distribution
     n2 = n > 100
     x2 = FINDGEN(n2+1)*(!X.CRANGE[1]-!X.CRANGE[0])/n2 + !X.CRANGE[0]
     y2 = call_function( func_name, x2 )
     oplot,x2,y2,lines=1,thick=2
 endif

 prob_kuiper, D, N, prob           ;Compute significance of Kuiper statistic

 return
 end