!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!