CHSCDF Subroutine

public subroutine CHSCDF(X, Nu, Cdf)

NAME

chscdf(3f) - [M_datapac:CUMULATIVE_DISTRIBUTION] compute the chi-square cumulative
distribution function

SYNOPSIS

   SUBROUTINE CHSCDF(X,Nu,Cdf)

    REAL(kind=wp),intent(in) :: X
    REAL(kind=wp),intent(out) :: Cdf
    INTEGER,intent(in) :: Nu

DESCRIPTION

CHSCDF(3f) computes the cumulative distribution function value for
the chi-squared distribution with integer degrees of freedom parameter
= NU.

This distribution is defined for all non-negative X.

The probability density function is given in the references below.

INPUT ARGUMENTS

X      The value at which the cumulative distribution function is to
       be evaluated. X should be non-negative.

NU     The integer number of degrees of freedom. NU should be positive.

OUTPUT ARGUMENTS

CDF    The cumulative distribution function value.

EXAMPLES

Sample program:

program demo_chscdf
use M_datapac, only : chscdf
implicit none
! call chscdf(x,y)
end program demo_chscdf

Results:

AUTHOR

The original DATAPAC library was written by James Filliben of the Statistical
Engineering Division, National Institute of Standards and Technology.

MAINTAINER

John Urban, 2022.05.31

LICENSE

CC0-1.0

REFERENCES

  • National Bureau of Standards Applied Mathematics Series 55, 1964, page 941, Formulae 26.4.4 and 26.4.5.
  • Johnson and Kotz, Continuous Univariate Distributions–1, 1970, page 176, Formula 28, and page 180, Formula 33.1.
  • Owen, Handbook of Statistical Tables, 1962, pages 50-55.
  • Pearson and Hartley, Biometrika Tables for Statisticians, Volume 1, 1954, pages 122-131.

NAME

chsplt(3f) - [M_datapac:LINE_PLOT] generate a Chi-square probability
plot

SYNOPSIS

   SUBROUTINE CHSPLT(X,N,Nu)

    REAL(kind=wp),intent(in) :: X(:)
    INTEGER,intent(in) :: N
    INTEGER,intent(in) :: Nu

DESCRIPTION

Chsplt(3f) generates a Chi-squared probability plot (with integer
degrees of freedom parameter value = NU).

The prototype Chi-squared distribution used herein is defIned for all
non-negative X, and its probability density function is given in the
references below.

As used herein, a probability plot for a distribution is a plot
of the ordered observations versus the order statistic medians for
that distribution.

The Chi-squared probability plot is useful in graphically testing
the composite (that is, location and scale parameters need not be
specified) hypothesis that the underlying distribution from which
the data have been randomly drawn is the Chi-squared distribution
with degrees of freedom parameter value = NU.

If the hypothesis is true, the probability plot should be near-linear.

a measure of such linearity is given by the calculated probability
plot correlation coefficient.

INPUT ARGUMENTS

X      The vector of (unsorted or sorted) observations.

N      The integer number of observations in the vector X.
       NU should be positive. The maximum allowable value of N for
       this subroutine is 7500.

NU     The integer number of degrees of freedom. NU should be positive.

OUTPUT

A one-page Chi-squared probability plot.

EXAMPLES

Sample program:

program demo_chsplt
use M_datapac, only : chsplt
implicit none
! call chsplt(x,y)
end program demo_chsplt

Results:

AUTHOR

The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.

MAINTAINER

John Urban, 2022.05.31

LICENSE

CC0-1.0

REFERENCES

  • Wilk, Gnanadesikan, and Huyett, ‘Probability Plots for the Gamma Distribution’, Technometrics, 1962, pages 1-15.
  • Filliben, ‘Techniques for Tail Length Analysis’, Proceedings of the Eighteenth Conference on the Design of Experiments in Army Research Development and Testing (Aberdeen, Maryland, October, 1972), pages 425-450.
  • Hahn and Shapiro, Statistical Methods in Engineering, 1967, pages 260-308.
  • Johnson and Kotz, Continuous Univariate Distributions–1, 1970, pages 166-206.
  • Hastings and Peacock, Statistical Distributions–A Handbook for Students and Practitioners, 1975, pages 46-51.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: X
integer, intent(in) :: Nu
real(kind=wp), intent(out) :: Cdf

Common Blocks

CAUPLT (subroutine)
CODE (subroutine)
DECOMP (subroutine)
DEMOD (subroutine)
DEXPLT (subroutine)
EV1PLT (subroutine)
EV2PLT (subroutine)
EXPPLT (subroutine)
EXTREM (subroutine)
EXTREM (subroutine)
FREQ (subroutine)
GAMPLT (subroutine)
GEOPLT (subroutine)
HFNPLT (subroutine)
INVXWX (subroutine)
LAMPLT (subroutine)
LGNPLT (subroutine)
LOGPLT (subroutine)
MEDIAN (subroutine)
norout (subroutine)
NORPLT (subroutine)
PARPLT (subroutine)
PLOTU (subroutine)
POIPLT (subroutine)
RUNS (subroutine)
SAMPP (subroutine)
SPCORR (subroutine)
TAIL (subroutine)
TPLT (subroutine)
TRIM (subroutine)
UNIPLT (subroutine)
WEIB (subroutine)
WEIPLT (subroutine)
WIND (subroutine)
"> common /BLOCK2_real32/

Type Attributes Name Initial
real :: WS(15000)

Source Code

SUBROUTINE CHSCDF(X,Nu,Cdf)
REAL(kind=wp),intent(in) :: X
REAL(kind=wp),intent(out) :: Cdf
INTEGER,intent(in) :: Nu

REAL(kind=wp) :: amean , anu , cdfn , danu , sd , spchi , u , z
INTEGER :: i , ibran , ievodd , imax , imin , nucut
DOUBLE PRECISION dx , chi , sum , term , ai , dcdfn
DOUBLE PRECISION dnu
DOUBLE PRECISION DSQRT , DEXP
DOUBLE PRECISION DLOG
DOUBLE PRECISION dfact , dpower
DOUBLE PRECISION dw
DOUBLE PRECISION d1 , d2 , d3
DOUBLE PRECISION term0 , term1 , term2 , term3 , term4
DOUBLE PRECISION b11
DOUBLE PRECISION b21
DOUBLE PRECISION b31 , b32
DOUBLE PRECISION b41 , b42 , b43
DATA nucut/1000/
DATA dpower/0.33333333333333D0/
DATA b11/0.33333333333333D0/
DATA b21/ - 0.02777777777778D0/
DATA b31/ - 0.00061728395061D0/
DATA b32/ - 13.0D0/
DATA b41/0.00018004115226D0/
DATA b42/6.0D0/
DATA b43/17.0D0/
!
!     CHECK THE INPUT ARGUMENTS FOR ERRORS
!
      IF ( Nu<=0 ) THEN
         WRITE (G_IO,99001)
99001    FORMAT (' ',                                                   &
     &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE CHSCDF SUBROU&
     &TINE IS NON-POSITIVE *****')
         WRITE (G_IO,99002) Nu
99002    FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****')
         Cdf = 0.0_wp
         RETURN
      ELSE
         IF ( X<0.0_wp ) THEN
            WRITE (G_IO,99003)
99003       FORMAT (' ',                                                &
     &'***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUMENT TO THE CHSC&
     &DF SUBROUTINE IS NEGATIVE *****')
            WRITE (G_IO,99004) X
99004       FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,    &
     &              ' *****')
            Cdf = 0.0_wp
            RETURN
         ELSE
!
!-----START POINT-----------------------------------------------------
!
            dx = X
            anu = Nu
            dnu = Nu
!
!     IF X IS NON-POSITIVE, SET CDF = 0.0 AND RETURN.
!     IF NU IS SMALLER THAN 10 AND X IS MORE THAN 200
!     STANDARD DEVIATIONS BELOW THE MEAN,
!     SET CDF = 0.0 AND RETURN.
!     IF NU IS 10 OR LARGER AND X IS MORE THAN 100
!     STANDARD DEVIATIONS BELOW THE MEAN,
!     SET CDF = 0.0 AND RETURN.
!     IF NU IS SMALLER THAN 10 AND X IS MORE THAN 200
!     STANDARD DEVIATIONS ABOVE THE MEAN,
!     SET CDF = 1.0 AND RETURN.
!     IF NU IS 10 OR LARGER AND X IS MORE THAN 100
!     STANDARD DEVIATIONS ABOVE THE MEAN,
!     SET CDF = 1.0 AND RETURN.
!
            IF ( X>0.0_wp ) THEN
               amean = anu
               sd = SQRT(2.0_wp*anu)
               z = (X-amean)/sd
               IF ( Nu>=10 .OR. z>=-200.0_wp ) THEN
                  IF ( Nu<10 .OR. z>=-100.0_wp ) THEN
                     IF ( Nu<10 .AND. z>200.0_wp ) GOTO 50
                     IF ( Nu>=10 .AND. z>100.0_wp ) GOTO 50
!
!     DISTINGUISH BETWEEN 3 SEPARATE REGIONS
!     OF THE (X,NU) SPACE.
!     BRANCH TO THE PROPER COMPUTATIONAL METHOD
!     DEPENDING ON THE REGION.
!     NUCUT HAS THE VALUE 1000.
!
                     IF ( Nu<nucut ) THEN
!
!     TREAT THE SMALL AND MODERATE DEGREES OF FREEDOM CASE
!     (THAT IS, WHEN NU IS SMALLER THAN 1000).
!     METHOD UTILIZED--EXACT FINITE SUM
!     (SEE AMS 55, page 941, FORMULAE 26.4.4 AND 26.4.5).
!
                        chi = DSQRT(dx)
                        ievodd = Nu - 2*(Nu/2)
                        IF ( ievodd==0 ) THEN
!
                           sum = 1.0D0
                           term = 1.0D0
                           imin = 2
                           imax = Nu - 2
                        ELSE
!
                           sum = 0.0D0
                           term = 1.0_wp/chi
                           imin = 1
                           imax = Nu - 1
                        ENDIF
!
                        IF ( imin<=imax ) THEN
                           DO i = imin , imax , 2
                              ai = i
                              term = term*(dx/ai)
                              sum = sum + term
                           ENDDO
                        ENDIF
!
                        sum = sum*DEXP(-dx/2.0D0)
                        IF ( ievodd/=0 ) THEN
                           sum = (DSQRT(2.0D0/G_pi_dp))*sum
                           spchi = chi
                           CALL NORCDF(spchi,cdfn)
                           dcdfn = cdfn
                           sum = sum + 2.0D0*(1.0D0-dcdfn)
                        ENDIF
                        GOTO 100
                     ELSEIF ( Nu>=nucut .AND. X<=anu ) THEN
!
!     TREAT THE CASE WHEN NU IS LARGE
!     (THAT IS, WHEN NU IS EQUAL TO OR GREATER THAN 1000)
!     AND X IS LESS THAN OR EQUAL TO NU.
!     METHOD UTILIZED--WILSON-HILFERTY APPROXIMATION
!     (SEE JOHNSON AND KOTZ, VOLUME 1, page 176, FORMULA 28).
!
                        dfact = 4.5D0*dnu
                        u = (((dx/dnu)**dpower)-1.0D0+(1.0D0/dfact))    &
     &                      *DSQRT(dfact)
                        CALL NORCDF(u,cdfn)
                        Cdf = cdfn
                        RETURN
                     ELSEIF ( Nu>=nucut .AND. X>anu ) THEN
!
!     TREAT THE CASE WHEN NU IS LARGE
!     (THAT IS, WHEN NU IS EQUAL TO OR GREATER THAN 1000)
!     AND X IS LARGER THAN NU.
!     METHOD UTILIZED--HILL'S ASYMPTOTIC EXPANSION
!     (SEE JOHNSON AND KOTZ, VOLUME 1, page 180, FORMULA 33.1).
!
                        dw = DSQRT(dx-dnu-dnu*DLOG(dx/dnu))
                        danu = DSQRT(2.0D0/dnu)
                        d1 = dw
                        d2 = dw**2
                        d3 = dw**3
                        term0 = dw
                        term1 = b11*danu
                        term2 = b21*d1*(danu**2)
                        term3 = b31*(d2+b32)*(danu**3)
                        term4 = b41*(b42*d3+b43*d1)*(danu**4)
                        u = term0 + term1 + term2 + term3 + term4
                        CALL NORCDF(u,cdfn)
                        Cdf = cdfn
                        GOTO 99999
                     ELSE
                        ibran = 1
                        WRITE (G_IO,99005) ibran
99005                   FORMAT (' ',                                    &
     &                      '*****INTERNAL ERROR IN CHSCDF SUBROUTINE--'&
     &                      ,                                           &
     &                  'IMPOSSIBLE BRANCH CONDITION AT BRANCH POINT = '&
     &                  ,I0)
                        RETURN
                     ENDIF
                  ENDIF
               ENDIF
            ENDIF
            Cdf = 0.0_wp
            RETURN
         ENDIF
 50      Cdf = 1.0_wp
         RETURN
      ENDIF
 100  Cdf = 1.0D0 - sum
      RETURN
!
99999 END SUBROUTINE CHSCDF