SUBROUTINE diptst(x, n, dip, xl, xu, ifault, gcm, lcm, mn, mj)
! Code converted using TO_F90 by Alan Miller
! Date: 2001-03-13 Time: 22:28:38
! ALGORITHM AS 217 APPL. STATIST. (1985) VOL.34, NO.3
! Does the dip calculation for an ordered vector X using the
! greatest convex minorant and the least concave majorant, skipping
! through the data using the change points of these distributions.
! It returns the dip statistic 'DIP' and the modal interval
! (XL, XU).
IMPLICIT NONE
REAL, INTENT(IN) :: x(:)
INTEGER, INTENT(IN) :: n
REAL, INTENT(OUT) :: dip
REAL, INTENT(OUT) :: xl
REAL, INTENT(OUT) :: xu
INTEGER, INTENT(OUT) :: ifault
INTEGER, INTENT(OUT) :: gcm(:)
INTEGER, INTENT(OUT) :: lcm(:)
INTEGER, INTENT(OUT) :: mn(:)
INTEGER, INTENT(OUT) :: mj(:)
INTEGER :: high, ic, icv, icva, icx, icxa, ig, igcm, igcm1, igcmx, &
ih, iv, ix, j, jb, je, jk, jr, k, kb, ke, kr, &
lcm1, lcmiv, lcmiv1, low, mjk, mjmjk, mnj, mnmnj, na
REAL :: a, b, const, d, dipnew, dl, du, dx, fn, t, temp
REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0
ifault = 1
IF (n <= 0) RETURN
ifault = 0
! Check if N = 1
IF (n /= 1) THEN
! Check that X is sorted
ifault = 2
DO k = 2, n
IF (x(k) < x(k-1)) RETURN
END DO
ifault = 0
! Check for all values of X identical,
! and for 1 < N < 4.
IF (x(n) > x(1) .AND. n >= 4) GO TO 20
END IF
xl = x(1)
xu = x(n)
dip = zero
RETURN
! LOW contains the index of the current estimate of the lower end
! of the modal interval, HIGH contains the index for the upper end.
20 fn = n
low = 1
high = n
dip = one / fn
xl = x(low)
xu = x(high)
! Establish the indices over which combination is necessary for the
! convex minorant fit.
mn(1) = 1
DO j = 2, n
mn(j) = j - 1
30 mnj = mn(j)
mnmnj = mn(mnj)
a = mnj - mnmnj
b = j - mnj
IF (mnj /= 1 .AND. (x(j)-x(mnj))*a >= (x(mnj)-x(mnmnj))*b) THEN
mn(j) = mnmnj
GO TO 30
END IF
END DO
! Establish the indices over which combination is necessary for the
! concave majorant fit.
mj(n) = n
na = n - 1
DO jk = 1, na
k = n - jk
mj(k) = k + 1
50 mjk = mj(k)
mjmjk = mj(mjk)
a = mjk - mjmjk
b = k - mjk
IF (mjk /= n .AND. (x(k)-x(mjk))*a >= (x(mjk)-x(mjmjk))*b) THEN
mj(k) = mjmjk
GO TO 50
END IF
END DO
! Start the cycling.
! Collect the change points for the GCM from HIGH to LOW.
70 ic = 1
gcm(1) = high
80 igcm1 = gcm(ic)
ic = ic + 1
gcm(ic) = mn(igcm1)
IF (gcm(ic) > low) GO TO 80
icx = ic
! Collect the change points for the LCM from LOW to HIGH.
ic = 1
lcm(1) = low
90 lcm1 = lcm(ic)
ic = ic + 1
lcm(ic) = mj(lcm1)
IF (lcm(ic) < high) GO TO 90
icv = ic
! ICX, IX, IG are counters for the convex minorant,
! ICV, IV, IH are counters for the concave majorant.
ig = icx
ih = icv
! Find the largest distance greater than 'DIP' between the GCM and
! the LCM from LOW to HIGH.
ix = icx - 1
iv = 2
d = zero
IF (icx == 2 .AND. icv == 2) THEN
d = one / fn
GO TO 120
END IF
100 igcmx = gcm(ix)
lcmiv = lcm(iv)
IF (igcmx <= lcmiv) THEN
! If the next point of either the GCM or LCM is from the LCM,
! calculate the distance here.
lcmiv1 = lcm(iv-1)
a = lcmiv - lcmiv1
b = igcmx - lcmiv1 - 1
dx = (x(igcmx)-x(lcmiv1)*a) / (fn*(x(lcmiv)-x(lcmiv1))) - b / fn
ix = ix - 1
IF (dx < d) GO TO 110
d = dx
ig = ix + 1
ih = iv
ELSE
! If the next point of either the GCM or LCM is from the GCM,
! calculate the distance here.
lcmiv = lcm(iv)
igcm = gcm(ix)
igcm1 = gcm(ix+1)
a = lcmiv - igcm1 + 1
b = igcm - igcm1
dx = a / fn - ((x(lcmiv)-x(igcm1))*b) / (fn*(x(igcm)-x(igcm1)))
iv = iv + 1
IF (dx >= d) THEN
d = dx
ig = ix + 1
ih = iv - 1
END IF
END IF
110 IF (ix < 1) ix = 1
IF (iv > icv) iv = icv
IF (gcm(ix) /= lcm(iv)) GO TO 100
120 IF (d >= dip) THEN
! Calculate the DIPs for the current LOW and HIGH.
! The DIP for the convex minorant.
dl = zero
IF (ig /= icx) THEN
icxa = icx - 1
DO j = ig, icxa
temp = one / fn
jb = gcm(j+1)
je = gcm(j)
IF (je-jb > 1) THEN
IF (x(je) /= x(jb)) THEN
a = je - jb
const = a / (fn*(x(je)-x(jb)))
DO jr = jb, je
b = jr - jb + 1
t = b / fn - (x(jr)-x(jb)) * const
IF (t > temp) temp = t
END DO
END IF
END IF
IF (dl < temp) dl = temp
END DO
END IF
! The DIP for the concave majorant.
du = zero
IF (ih /= icv) THEN
icva = icv - 1
DO k = ih, icva
temp = one / fn
kb = lcm(k)
ke = lcm(k+1)
IF (ke-kb > 1) THEN
IF (x(ke) /= x(kb)) THEN
a = ke - kb
const = a / (fn*(x(ke)-x(kb)))
DO kr = kb, ke
b = kr - kb - 1
t = (x(kr)-x(kb)) * const - b / fn
IF (t > temp) temp = t
END DO
END IF
END IF
IF (du < temp) du = temp
END DO
END IF
! Determine the current maximum.
dipnew = dl
IF (du > dl) dipnew = du
IF (dip < dipnew) dip = dipnew
low = gcm(ig)
high = lcm(ih)
! Recycle
GO TO 70
END IF
dip = half * dip
xl = x(low)
xu = x(high)
RETURN
END SUBROUTINE diptst