!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! wav.f - weighted averages and weighted circular averages ! by Andy Allinger, 2016, released to the public domain ! This program may be used by any person for any purpose. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! weighted mean !__Name____Type_____In/Out________Description___________________________ ! X[N] Real In data array ! W[N] Real In importance weights ! N Integer In length of arrays ! XWM Real Out weighted median !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE WMEAN (X, W, N, XWM) IMPLICIT NONE INTEGER N REAL X, W, XWM DIMENSION X(N), W(N) ! local variables INTEGER I REAL SW, SWX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! begin !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (N < 1) RETURN SW = 0. SWX = 0. DO I = 1, N SW = SW + W(I) SWX = SWX + W(I) * X(I) END DO IF (SW > 0.) XWM = SWX / SW RETURN END ! of WMEAN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! weighted median !__Name____Type_____In/Out________Description___________________________ ! X[N] Real In data array ! W[N] Real In importance weights ! IND[N] Integer Work index array ! N Integer In length of arrays ! XWM Real Out weighted median !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE WMED (X, W, IND, N, XWM) IMPLICIT NONE INTEGER IND, N REAL X, W, XWM DIMENSION X(N), W(N), IND(N) ! local variables INTEGER I REAL SW, HSW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! begin !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (N < 1) RETURN IF (N .EQ. 1) THEN XWM = X(1) RETURN END IF SW = 0. DO I = 1, N SW = SW + W(I) END DO IF (SW .LE. 0.) RETURN HSW = SW / 2. CALL QSORTR (X, IND, N) CALL RORDER (W, IND, N) CALL RORDER (X, IND, N) SW = 0. DO I = 1, N SW = SW + W(I) IF (SW > HSW) THEN XWM = X(I) RETURN ELSE IF (SW .EQ. HSW) THEN XWM = 0.5 * (X(I) + X(I+1)) RETURN END IF END DO END ! of WMED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! weighted cyclic mean for phase only ! !__Name____Type_____In/Out__Description_________________________________ ! X[N] Real In data array ! W[N] Real In importance weights ! IND[N] Integer Work permutation vector ! N Integer In length of arrays ! T Real In period ! AVE[N] Real Out average ! TIES Integer Out how many values in AVE[] !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE WCYMEN (X, W, IND, N, T, AVE, TIES) IMPLICIT NONE INTEGER IND, N, TIES REAL X, W, T, AVE DIMENSION X(N), W(N), IND(N), AVE(N) INTEGER I REAL & MACHINE, ! function for EPSILON and HUGE & R, ! temporary scalar X[i] & TOL ! a small number DOUBLE PRECISION & SW, SWX, SWX2, ! sum of weights, sum of squares & V, VMIN, ! n @ variance, minimal n @ variance & TT, T2, ! period, period squared & RR, ! double precision of R & WW ! double precision of W(I) ! begin IF (N .EQ. 1) THEN TIES = 1 AVE(1) = X(1) RETURN END IF TOL = SQRT(T * MACHINE(1)) TT = DBLE(T) T2 = TT **2 SW = 0. SWX = 0. SWX2 = 0. ! move all data to the range 0...T, sort DO I = 1, N R = MOD(X(I), T) IF (R < 0.) R = R + T WW = DBLE(W(I)) RR = DBLE(R) SW = SW + WW ! accumulate sums SWX = SWX + WW * RR SWX2 = SWX2 + WW * RR **2 X(I) = R END DO CALL QSORTR (X, IND, N) CALL RORDER (W, IND, N) CALL RORDER (X, IND, N) ! find phase of seam of lowest weighted deviation VMIN = MACHINE(3) TIES = 0 DO I = 1, N V = SWX2 - SWX **2 / SW ! SW @ weighted variance IF (VMIN - V > TOL) THEN ! new low deviation TIES = 0 VMIN = V END IF IF (ABS(V - VMIN) < TOL) THEN ! tie for new low deviation TIES = TIES + 1 AVE(TIES) = SNGL(SWX / SW) ! weighted mean END IF ! apply update formula to sum and sum of squares WW = DBLE(W(I)) RR = DBLE(X(I)) SWX = SWX + WW * TT SWX2 = SWX2 + WW * (2D0 * RR * TT + T2) END DO DO I = 1, TIES R = AVE(I) IF (R .GE. T) R = R - T ! subtract if out of 0...T AVE(I) = R END DO END ! of WCYMEN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! weighted cyclic median for phase only ! !NOTE THAT! unlike the unweighted cyclic median, there is not ! an update formula available. The mean deviation is computed ! N repetitions, in O(N^2) time ! ! !__Name____Type_____In/Out________Description___________________________ ! X[N] Real In data array ! W[N] Real In importance weights ! IND[N] Integer Work permuatation vector ! N Integer In length of X ! T Real In period ! AVE[N] Real Out average ! TIES Integer Out how many values in AVE[] !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE WCYMED (X, W, IND, N, T, AVE, TIES) IMPLICIT NONE INTEGER IND, N, TIES REAL X, W, T, AVE DIMENSION X(N), W(N), IND(N), AVE(N) INTEGER I, J, K, ! counters & J1 ! position where median was found REAL & MACHINE, ! function for EPSILON and HUGE & R, ! temporary copy of X & TOL ! a small number DOUBLE PRECISION & SW, HSW, ! sum of W[], half the sum & S1, SMIN ! deviation, minimal deviation ! begin IF (N .EQ. 1) THEN TIES = 1 AVE(1) = X(1) RETURN END IF TOL = T * MACHINE(1) ! move all data to the range 0...t, sort SW = 0. R = 0. ! avoid compiler warning DO I = 1, N R = MOD(X(I), T) IF (R < 0.) R = R + T X(I) = R SW = SW + DBLE(W(I)) END DO IF (SW .LE. 0.) RETURN HSW = SW / 2. CALL QSORTR (X, IND, N) CALL RORDER (W, IND, N) CALL RORDER (X, IND, N) ! initial value of the cyclical median SW = 0. J1 = 1 ! avoid compiler warning DO I = 1, N SW = SW + DBLE(W(I)) IF (SW > HSW) THEN R = X(I) J1 = I GO TO 10 ! done seeking initial median ELSE IF (SW .EQ. HSW) THEN R = 0.5 * (X(I) + X(I+1)) J1 = I GO TO 10 ! done seeking initial median END IF END DO 10 CONTINUE ! sum of weighted absolute deviations S1 = 0. DO I = 1, N S1 = S1 + DPROD(W(I), ABS(X(I) - R)) END DO ! find phase of seam of lowest deviation SMIN = MACHINE(3) TIES = 0 DO I = 1, N IF (SMIN - S1 > TOL) THEN ! new low deviation TIES = 0 SMIN = S1 END IF IF (ABS(S1 - SMIN) < TOL) THEN ! tie for new low deviation TIES = TIES + 1 AVE(TIES) = R END IF ! recompute median and deviation X(I) = X(I) + T ! advance seam SW = SW - DBLE(W(I)) ! don't include in total K = J1 ! possible for median to stay same DO J = 1, N IF (SW > HSW) THEN R = X(K) J1 = K GO TO 20 ELSE IF (SW .EQ. HSW) THEN IF (K .EQ. N) THEN R = 0.5 * (X(N) + X(1)) ELSE R = 0.5 * (X(K) + X(K+1)) END IF J1 = K GO TO 20 END IF K = K + 1 IF (K > N) K = 1 SW = SW + DBLE(W(K)) END DO ! next J 20 CONTINUE ! sum of absolute deviations S1 = 0. DO J = 1, N S1 = S1 + DPROD(W(J), ABS(X(J) - R)) END DO END DO ! next I DO I = 1, TIES R = AVE(I) IF (R .GE. T) R = R - T AVE(I) = R END DO END ! of WCYMED *!!!!!!!!!!!!!!!!!!!!!!! End of file wav.f !!!!!!!!!!!!!!!!!!!!!!!!!!!!!