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