RUNS – M_datapac

public subroutine RUNS(X, N)

NAME

runs(3f) - [M_datapac:ANALYSIS] perform a runs test

SYNOPSIS

   SUBROUTINE RUNS(X,N)

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

DESCRIPTION

RUNS(3f) performs a runs analysis of the data in the input vector x.

This runs analysis is a useful distribution-free test of the randomness
of a data set.

The analysis consists of first determining the observed number of
runs from the data, and then computing the expected number of runs,
the standard deviation of the number of runs, and the resulting
standardized statistic for the number of runs for runs of various
lengths.

This is done for runs up, runs down, and runs up and down.

INPUT ARGUMENTS

X  The precision vector of (unsorted or sorted) observations.

N  The integer number of observations in the vector x.

   restrictions-- The maximum allowable value of n for this subroutine
   is 15000.

OUTPUT

4 pages of automatic printout consisting of the observed number,
expected number, standard deviation and resulting standardized
statistic for runs of various lengths, and the cumulative frequency.

EXAMPLES

Sample program:

program demo_runs
use M_datapac, only : runs
implicit none
! call runs(x,y)
end program demo_runs

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

Levene and Wolfowitz, Annals of Mathematical Statistics, 1944, pages 58-69; especially pages 60, 63, and 64.
Bradley, Distribution-free Statistical Tests, 1968, Chapter 12, pages 271-282.

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in)			::	X(:)
integer,	intent(in)			::	N

Common Blocks

CAUPLT (subroutine)
CHSCDF (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)
SAMPP (subroutine)
SPCORR (subroutine)
TAIL (subroutine)
TPLT (subroutine)
TRIM (subroutine)
UNIPLT (subroutine)
WEIB (subroutine)
WEIPLT (subroutine)
WIND (subroutine)
"> common /BLOCK2_real64/

Type	Attributes		Name		Initial
real		::	WS(15000)

Source Code

SUBROUTINE RUNS(X,N)
REAL(kind=wp),intent(in) :: X(:)
INTEGER,intent(in)       :: N
REAL(kind=wp) :: ai, an, anrdl, anrdlg, anrtl, anrtlg, anrul, anrulg
REAL(kind=wp) :: c1, c2, c3, c4, den, enrtl, enrtlg, enrul, enrulg, hold, snrtl, snrtlg
REAL(kind=wp) :: snrul, snrulg, stat, WS, Y, znrdl, znrdlg, znrtl, znrtlg, znrul, znrulg
INTEGER       :: i, imax, ip1, iupper, j, jp1, lendn, lenup, maxlnd, maxlnt, maxlnu, nm1, nneg, npos, nrdl, nrdlg, nrtl, nrtlg
INTEGER       :: nrul, nrulg, nzer

DIMENSION Y(15000)
DIMENSION nrul(16) , nrdl(16) , nrtl(16) , nrulg(16) , nrdlg(16)
DIMENSION nrtlg(16)
DIMENSION enrul(16) , enrtl(16) , enrulg(16) , enrtlg(16)
DIMENSION snrul(16) , snrtl(16) , snrulg(16) , snrtlg(16)
DIMENSION znrul(16) , znrdl(16) , znrtl(16) , znrulg(16) , znrdlg(16)
DIMENSION znrtlg(16)
DIMENSION c1(15) , c2(15) , c3(15) , c4(15)
DIMENSION anrul(16) , anrdl(16) , anrtl(16)
DIMENSION anrulg(16) , anrdlg(16) , anrtlg(16)
COMMON /BLOCK2_real64/ WS(15000)
EQUIVALENCE (Y(1),WS(1))

      DATA c1(1) , c1(2) , c1(3) , c1(4) , c1(5) , c1(6) , c1(7) ,      &
     &     c1(8) , c1(9) , c1(10) , c1(11) , c1(12) , c1(13) , c1(14) , &
     &     c1(15)/0.4236111111E+00_wp , .1126675485E+00_wp , .4191688713E-01_wp , &
     &     .1076912487E-01_wp , .2003959238E-02_wp , .3023235799E-03_wp ,        &
     &     .3911555473E-04_wp , .4459038843E-05_wp , .4551105210E-06_wp ,        &
     &     .4207466837E-07_wp , .3555930927E-08_wp , .2768273257E-09_wp ,        &
     &     .1997821524E-10_wp , .1343876568E-11_wp , .8465610177E-13_wp/
      DATA c2(1) , c2(2) , c2(3) , c2(4) , c2(5) , c2(6) , c2(7) ,      &
     &     c2(8) , c2(9) , c2(10) , c2(11) , c2(12) , c2(13) , c2(14) , &
     &     c2(15)/ - .4819444444E+00_wp , -.1628284832E+00_wp ,               &
     &     -.9690696649E-01_wp , -.3778106786E-01_wp , -.9289228716E-02_wp ,     &
     &     -.1724429252E-02_wp , -.2638557888E-03_wp , -.3466965096E-04_wp ,     &
     &     -.4004129153E-05_wp , -.4130382587E-06_wp , -.3851876069E-07_wp ,     &
     &     -.3279103786E-08_wp , -.2568491117E-09_wp , -.1863433868E-10_wp ,     &
     &     -.1259220466E-11_wp/
      DATA c3(1) , c3(2) , c3(3) , c3(4) , c3(5) , c3(6) , c3(7) ,      &
     &     c3(8) , c3(9) , c3(10) , c3(11) , c3(12) , c3(13) , c3(14) , &
     &     c3(15)/.1777777778E+00_wp , .7916666667E-01_wp , .4738977072E-01_wp , &
     &     .1274801587E-01_wp , .2338606059E-02_wp , .3461358734E-03_wp ,        &
     &     .4407121770E-04_wp , .4960020603E-05_wp , .5010387575E-06_wp ,        &
     &     .4592883352E-07_wp , .3854170274E-08_wp , .2982393839E-09_wp ,        &
     &     .2141205844E-10_wp , .1433843200E-11_wp , .8996663214E-13_wp/
      DATA c4(1) , c4(2) , c4(3) , c4(4) , c4(5) , c4(6) , c4(7) ,      &
     &     c4(8) , c4(9) , c4(10) , c4(11) , c4(12) , c4(13) , c4(14) , &
     &     c4(15)/ - .3222222222E+00_wp , -.5972222222E-01_wp ,               &
     &     -.1130268959E+00_wp , -.4696428571E-01_wp , -.1123273065E-01_wp ,     &
     &     -.2025170849E-02_wp , -.3029410411E-03_wp , -.3912824548E-04_wp ,     &
     &     -.4459234519E-05_wp , -.4551128785E-06_wp , -.4207469124E-07_wp ,     &
     &     -.3555931110E-08_wp , -.2768273269E-09_wp , -.1997821525E-10_wp ,     &
     &     -.1343876568E-11_wp/
!
      iupper = 15000
!
!     CHECK THE INPUT ARGUMENTS FOR ERRORS
!
      IF ( N<1 .OR. N>iupper ) THEN
         WRITE (G_IO,99001) iupper
99001    FORMAT (' ***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO RUNS(3f) IS OUTSIDE THE ALLOWABLE (1,',I0,') INTERVAL *****')
         WRITE (G_IO,99002) N
99002    FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****')
         RETURN
      ELSEIF ( N==1 ) THEN
         WRITE (G_IO,99003)
99003    FORMAT (' ***** FATAL ERROR--         THE SECOND INPUT ARGUMENT TO THE RUNS SUBROUTINE HAS THE VALUE 1 *****')
         RETURN
      ELSE
         hold = X(1)
         DO i = 2 , N
            IF ( X(i)/=hold ) GOTO 50
         ENDDO
         WRITE (G_IO,99004) hold
99004    FORMAT (' ***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUMENT (A VECTOR) TO RUNS(3f) HAS ALL ELEMENTS = ',E15.8,' *****')
!
!-----START POINT-----------------------------------------------------
!
 50      an = N
!
!     FORM THE SEQUENTIAL DIFFERENCE TABLE
!
         nm1 = N - 1
         DO i = 1 , nm1
            ip1 = i + 1
            Y(i) = X(ip1) - X(i)
         ENDDO
!
!     ZERO-OUT THE 6 'NUMBER OF RUNS' VECTORS
!
         DO i = 1 , 16
            nrul(i) = 0
            nrdl(i) = 0
            nrtl(i) = 0
            nrulg(i) = 0
            nrdlg(i) = 0
            nrtlg(i) = 0
         ENDDO
!
!     DETERMINE THE NUMBER OF RUNS UP OF LENGTH EXACTLY I
!     AND THE NUMBER OF RUNS DOWN OF LENGTH EXACTLY I
!     DETERMINE THE LENGTH OF THE LONGEST RUN UP
!     AND THE LENGTH OF THE LONGEST RUN DOWN
!
         lenup = 0
         lendn = 0
         maxlnu = 0
         maxlnd = 0
         DO i = 1 , nm1
            IF ( Y(i)==0.0_wp .AND. lenup>=1 ) lenup = lenup + 1
            IF ( Y(i)==0.0_wp .AND. lendn>=1 ) lendn = lendn + 1
            IF ( Y(i)==0.0_wp .AND. lenup==0 .AND. lendn==0 )              lenup = lenup + 1
            IF ( Y(i)>0.0_wp .AND. lendn>=1 .AND. lendn<=15 ) nrdl(lendn)  = nrdl(lendn) + 1
            IF ( Y(i)>0.0_wp .AND. lendn>=1 .AND. lendn>=16 ) nrdl(16)     = nrdl(16) + 1
            IF ( Y(i)>0.0_wp ) lendn = 0
            IF ( Y(i)>0.0_wp ) lenup = lenup + 1
            IF ( Y(i)<0.0_wp .AND. lenup>=1 .AND. lenup<=15 ) nrul(lenup)  = nrul(lenup) + 1
            IF ( Y(i)<0.0 .AND. lenup>=1 .AND. lenup>=16 ) nrul(16)     = nrul(16) + 1
            IF ( Y(i)<0.0 ) lenup = 0
            IF ( Y(i)<0.0 ) lendn = lendn + 1
            IF ( i==nm1 .AND. lendn>=1 .AND. lendn<=15 ) nrdl(lendn)    = nrdl(lendn) + 1
            IF ( i==nm1 .AND. lendn>=1 .AND. lendn>=16 ) nrdl(16)       = nrdl(16) + 1
            IF ( i==nm1 .AND. lenup>=1 .AND. lenup<=15 ) nrul(lenup)    = nrul(lenup) + 1
            IF ( i==nm1 .AND. lenup>=1 .AND. lenup>=16 ) nrul(16)       = nrul(16) + 1
            IF ( lenup>maxlnu ) maxlnu = lenup
            IF ( lendn>maxlnd ) maxlnd = lendn
         ENDDO
!
!     DETERMINE THE NUMBER OF RUNS TOTAL OF LENGTH EXACTLY I
!     AND THE LENGTH OF THE LONGEST RUN UP OR DOWN
!
         DO i = 1 , 16
            nrtl(i) = nrul(i) + nrdl(i)
         ENDDO
         maxlnt = maxlnu
         IF ( maxlnd>maxlnu ) maxlnt = maxlnd
!
!     DETERMINE THE NUMBER OF RUNS UP OF LENGTH I OR MORE
!     AND THE NUMBER OF RUNS DOWN OF LENGTH I OR MORE
!     AND THE NUMBER OF RUNS TOTAL OF LENGTH I OR MORE
!
         nrulg(16) = nrul(16)
         nrdlg(16) = nrdl(16)
         nrtlg(16) = nrtl(16)
         DO i = 1 , 15
            j = 16 - i
            jp1 = j + 1
            nrulg(j) = nrulg(jp1) + nrul(j)
            nrdlg(j) = nrdlg(jp1) + nrdl(j)
            nrtlg(j) = nrtlg(jp1) + nrtl(j)
         ENDDO
!
!     DETERMINE THE NUMBER OF POSITIVE, ZERO, AND NEGATIVE ENTRIES
!     IN THE DIFFERENCE TABLE.  IF RANDOM, THE NUMBER OF POSITIVE SHOULD BE
!     APPROXIMATELY EQUAL TO THE NUMBER OF NEGATIVE
!
         nneg = 0
         nzer = 0
         npos = 0
         DO i = 1 , nm1
            IF ( Y(i)<0.0_wp ) nneg = nneg + 1
            IF ( Y(i)==0.0_wp ) nzer = nzer + 1
            IF ( Y(i)>0.0_wp ) npos = npos + 1
         ENDDO
!
!     COMPUTE THE EXPECTED NUMBER OF RUNS UP OF LENGTH EXACTLY I =
!     THE EXPECTED NUMBER OF RUNS DOWN OF LENGTH EXACTLY I =
!     ONE HALF THE EXPECTED NUMBER OF RUNS TOTAL OF LENGTH EXACTLY I
!
         den = 6.0_wp
         DO i = 1 , 15
            ai = i
            enrul(i) = an*(ai*ai+3.0_wp*ai+1.0_wp) - (ai*ai*ai+3.0_wp*ai*ai-ai-4.0_wp)
            den = den*(ai+3.0_wp)
            enrul(i) = enrul(i)/den
            enrtl(i) = 2.0_wp*enrul(i)
         ENDDO
!
!     COMPUTE THE EXPECTED NUMBER OF RUNS UP OF LENGTH I OR MORE =
!     THE EXPECTED NUMBER OF RUNS DOWN OF LENGTH I OR MORE =
!     ONE HALF THE EXPECTED NUMBER OF RUNS TOTAL OF LENGTH I OR MORE
!
         den = 2.0_wp
         DO i = 1 , 15
            ai = i
            enrulg(i) = an*(ai+1.0_wp) - (ai*ai+ai-1.0_wp)
            den = den*(ai+2.0_wp)
            enrulg(i) = enrulg(i)/den
            enrtlg(i) = 2.0_wp*enrulg(i)
         ENDDO
!
!     COMPUTE THE STANDARD DEV. OF THE NUMBER OF RUNS UP OF LENGTH EXACTLY I =
!     THE STANDARD DEV. OF THE NUMBER OF RUNS DOWN OF LENGTH EXACTLY I =
!     SQRT(0.5)* THE STAND. DEV. OF THE NUMBER OF RUNS TOTAL OF LENGTH EXACTLY I
!
         DO i = 1 , 15
            snrtl(i) = SQRT(c1(i)*an+c2(i))
            snrul(i) = SQRT(0.5_wp)*snrtl(i)
         ENDDO
!
!     COMPUTE THE STAND. DEV. OF THE NUMBER OF RUNS UP OF LENGTH I OR MORE =
!     THE STAND. DEV. OF THE NUMBER OF RUNS DOWN OF LENGTH I OR MORE =
!     SQRT(0.5)* THE STAND. DEV. OF THE NUMBER OF RUNS TOTAL OF LENGTH I OR MORE
!
         DO i = 1 , 15
            snrtlg(i) = SQRT(c3(i)*an+c4(i))
            snrulg(i) = SQRT(0.5_wp)*snrtlg(i)
         ENDDO
!
!     FORM Z STATISTICS
!
         DO i = 1 , 15
            stat = nrul(i)
            znrul(i) = (stat-enrul(i))/snrul(i)
            stat = nrdl(i)
            znrdl(i) = (stat-enrul(i))/snrul(i)
            stat = nrtl(i)
            znrtl(i) = (stat-enrtl(i))/snrtl(i)
            stat = nrulg(i)
            znrulg(i) = (stat-enrulg(i))/snrulg(i)
            stat = nrdlg(i)
            znrdlg(i) = (stat-enrulg(i))/snrulg(i)
            stat = nrtlg(i)
            znrtlg(i) = (stat-enrtlg(i))/snrtlg(i)
         ENDDO
!
         DO i = 1 , 15
            anrul(i) = nrul(i)
            anrdl(i) = nrdl(i)
            anrtl(i) = nrtl(i)
            anrulg(i) = nrulg(i)
            anrdlg(i) = nrdlg(i)
            anrtlg(i) = nrtlg(i)
         ENDDO
!
!     WRITE EVERYTHING OUT
!
         imax = 15
         WRITE (G_IO,99024)
         WRITE (G_IO,99005)
99005    FORMAT (' ',48X,'RUNS UP')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99006)
99006    FORMAT (' ',27X, 'STATISTIC = NUMBER OF RUNS UP    OF LENGTH EXACTLY I')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99022)
         WRITE (G_IO,99025)
         DO i = 1 , imax
            WRITE (G_IO,99023) i , anrul(i) , enrul(i) , snrul(i) , znrul(i)
         ENDDO
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99007)
99007    FORMAT (' ',27X, 'STATISTIC = NUMBER OF RUNS UP    OF LENGTH I OR MORE')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99022)
         WRITE (G_IO,99025)
         DO i = 1 , imax
            WRITE (G_IO,99023) i , anrulg(i) , enrulg(i) , snrulg(i) ,   znrulg(i)
         ENDDO
         WRITE (G_IO,99024)
         WRITE (G_IO,99008)
99008    FORMAT (' ',48X,'RUNS DOWN')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99009)
99009    FORMAT (' ',27X, 'STATISTIC = NUMBER OF RUNS DOWN  OF LENGTH EXACTLY I')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99022)
         WRITE (G_IO,99025)
         DO i = 1 , imax
            WRITE (G_IO,99023) i , anrdl(i) , enrul(i) , snrul(i) , znrdl(i)
         ENDDO
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99010)
99010    FORMAT (' ',27X, 'STATISTIC = NUMBER OF RUNS DOWN  OF LENGTH I OR MORE')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99022)
         WRITE (G_IO,99025)
         DO i = 1 , imax
            WRITE (G_IO,99023) i , anrdlg(i) , enrulg(i) , snrulg(i) , znrdlg(i)
         ENDDO
         WRITE (G_IO,99024)
         WRITE (G_IO,99011)
99011    FORMAT (' ',40X,'RUNS TOTAL = RUNS UP + RUNS DOWN')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99012)
99012    FORMAT (' ',27X, 'STATISTIC = NUMBER OF RUNS TOTAL OF LENGTH EXACTLY I')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99022)
         WRITE (G_IO,99025)
         DO i = 1 , imax
            WRITE (G_IO,99023) i , anrtl(i) , enrtl(i) , snrtl(i) , znrtl(i)
         ENDDO
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99013)
99013    FORMAT (' ',27X,'STATISTIC = NUMBER OF RUNS TOTAL OF LENGTH I OR MORE')
         WRITE (6,99025)
         WRITE (6,99025)
         WRITE (G_IO,99022)
         WRITE (G_IO,99025)
         DO i = 1 , imax
            WRITE (G_IO,99023) i , anrtlg(i) , enrtlg(i) , snrtlg(i) ,   znrtlg(i)
         ENDDO
         WRITE (G_IO,99024)
         WRITE (G_IO,99014) maxlnu
99014    FORMAT (' ','LENGTH OF THE LONGEST RUN UP         = ',I0)
         WRITE (G_IO,99015) maxlnd
99015    FORMAT (' ','LENGTH OF THE LONGEST RUN DOWN       = ',I0)
         WRITE (G_IO,99016) maxlnt
99016    FORMAT (' ','LENGTH OF THE LONGEST RUN UP OR DOWN = ',I0)
         WRITE (G_IO,99025)
         WRITE (G_IO,99017) npos
99017    FORMAT (' ','NUMBER OF POSITIVE DIFFERENCES = ',I0)
         WRITE (G_IO,99018) nneg
99018    FORMAT (' ','NUMBER OF NEGATIVE DIFFERENCES = ',I0)
         WRITE (G_IO,99019) nzer
99019    FORMAT (' ','NUMBER OF ZERO     DIFFERENCES = ',I0)
99020    FORMAT (' ',2(I4,2X,F7.1,2X,F8.4,2X,F8.4,2X,F8.2,8X))
99021    FORMAT (' ',I0,2X,I0,2X,I0)
      ENDIF
99022 FORMAT (' ', 'I = LENGTH OF RUN         VALUE OF STAT        EXP(STAT)            SD(STAT)    (STAT-EXP(STAT))/SD(STAT)')
99023 FORMAT (' ',4X,I4,13X,6X,F7.1,13X,F8.4,12X,F8.4,11X,F8.2)
99024 FORMAT ('1')
99025 FORMAT (' ')
END SUBROUTINE RUNS

RUNS Subroutine

Contents

Common Blocks

Source Code

public subroutine RUNS(X, N)

NAME

SYNOPSIS

DESCRIPTION

INPUT ARGUMENTS

OUTPUT

EXAMPLES

AUTHOR

MAINTAINER

LICENSE

REFERENCES

Arguments

Common Blocks

Source Code