*-----------------------------------------------------------------------
* average.f - routines for ordinary, cyclic, and weighted averages
* by Andy Allinger, 2013-2017, released to the public domain
* This file may be used by any person for any purpose.
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
* MEAN - arithmetic mean
* X[N] Real In data array
* N Integer In length of X
* AVE Real Out average
*-----------------------------------------------------------------------
SUBROUTINE MEAN (X, N, AVE)
IMPLICIT NONE
INTEGER N
REAL X, AVE
DIMENSION X(N)
IF (N < 1) RETURN
CALL ADD (X, N, AVE)
AVE = AVE / FLOAT(N)
RETURN
END ! of MEAN
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* CYMEAN - cyclic mean
* by Andy Allinger, 2013, released to the public domain
* origin: David Radcliffe, 2005 (or earlier)
*
* He proposed the algorithm:
* 1. Find the average and standard deviation of the given numbers.
* 2. Increase the smallest number by 360.
* 3. Repeat steps 1 and 2 until all numbers are greater than 360.
* 4. Choose the average that yields the smallest standard deviation.
* 5. If the average is greater than 360 then subtract 360.
*
* Notice that unlike the linear mean, the cyclic mean does not
* always have a distinct value. Instead, the array ave[] contains in
* its beginning entries the various averages, and the integer ties denotes
* how many averages were found.
* WARNING: X is returned shifted and sorted
*
* X [n] Real In data array
* n Integer In length of X
* t Real In period
* AVE [n] Real Out average
* ties Integer Out how many values in ave[]
* dev Real Out standard deviation
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE CYMEAN (X, n, t, AVE, ties, dev)
IMPLICIT NONE
INTEGER n, ties
REAL X, t, AVE, dev
DIMENSION X(n), AVE(n)
INTEGER i
REAL
& rmachine, ! function for EPSILON()
& r, ! temporary scalar X[i]
& tol ! a small number
DOUBLE PRECISION
& nn, ! double precision of n
& s, s2, ! sum, sum of squares
& v, vmin, ! n @ variance, minimal n @ variance
& tt, t2, ! period, period squared
& rr ! double precision of r
* begin
IF (n .EQ. 1) THEN
ties = 1
AVE(1) = X(1)
dev = 0.
RETURN
END IF
nn = DBLE(n)
tt = DBLE(t)
! tol = SQRT(t * rmachine(1))
tol = t * rmachine(1)
t2 = tt **2
s = 0.
s2 = 0.
* move all data to the range 0...t, sort
DO i = 1, n
r = X(i)
r = MOD(r, t)
IF (r < 0.) r = r + t
rr = DBLE(r)
s = s + rr
s2 = s2 + rr **2 ! accumulate sums
X(i) = r
END DO
CALL SSORT (X, n)
* find phase of seam of lowest deviation
vmin = rmachine(3)
ties = 0
DO i = 1, n
v = s2 - s **2 / nn ! n @ variance
IF (vmin - v > tol) THEN ! new low deviation
ties = 0
vmin = v
END IF
IF (ABS(v - vmin) .LE. tol) THEN ! tie for new low deviation
ties = ties + 1
AVE(ties) = SNGL(s / nn)
END IF
* apply update formula to sum and sum of squares
s = s + tt
s2 = s2 + 2D0 * DBLE(X(i)) * tt + t2
END DO
* subtract if out of 0...t
DO i = 1, ties
IF (AVE(i) .GE. t) AVE(i) = AVE(i) - t
END DO
dev = SNGL(SQRT(vmin / nn)) ! (n@variance / n)^(1/2)
RETURN
END ! of CYMEAN
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* CYMED - cyclic median, a variation of Radcliffe's algorithm
* by Andy Allinger, 2013, released to the public domain
* The cyclic median does not always have a distinct value.
* WARNING: X is returned shifted and sorted
*
* X [n] Real In data array
* n Integer In length of X
* t Real In period
* AVE [n] Real Out average
* ties Integer Out how many values in AVE[]
* dev Real Out deviation
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE CYMED (X, n, t, AVE, ties, dev)
IMPLICIT NONE
INTEGER n, ties
REAL X, t, AVE, dev
DIMENSION X(n), AVE(n)
INTEGER i, ! count
& j0, j1 ! positions where to find current estimate of median
REAL
& rmachine, ! function when EPSILON() not available
& r, ! temporary copy of X
& tol ! a small number
DOUBLE PRECISION
& s1, smin, ! sum of deviations, minimum sum deviation
& tt ! period
LOGICAL odd ! even or odd n
* begin
IF (n .EQ. 1) THEN
ties = 1
AVE(1) = X(1)
dev = 0.
RETURN
END IF
j1 = n / 2 + 1
odd = MOD(n, 2) .EQ. 1 ! is n an odd number?
tt = DBLE(t)
tol = t * rmachine(1)
* move all data to the range 0...t, sort
DO i = 1, n
r = MOD(X(i), t)
IF (r < 0.) r = r + t
X(i) = r
END DO
CALL SSORT (X, n)
* initial value of cyclic median
IF (odd) THEN
r = X(j1)
ELSE
j0 = j1 - 1
r = (X(j0) + X(j1)) / 2. ! split the difference
END IF
* sum of absolute deviations
s1 = 0.
DO i = 1, n
s1 = s1 + DBLE(ABS(X(i) - r))
END DO
* find phase of seam of lowest deviation
smin = rmachine(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) .LE. tol) THEN ! tie for new low deviation
ties = ties + 1
IF (odd) THEN
AVE(ties) = X(j1)
ELSE ! two points required for even case
j0 = j1 - 1
IF (j0 .EQ. 0) j0 = n ! j0 is the prior point
AVE(ties) = (X(j0) + X(j1)) / 2.
END IF
END IF
* apply update formula
r = X(i)
IF (odd) THEN
j0 = j1 + 1
IF (j0 > n) j0 = 1 ! new and current medians
s1 = s1 + 2. * r + tt - DBLE(X(j0) + X(j1)) ! j0 is the new median
ELSE
s1 = s1 + 2. * r + tt - 2. * X(j1)
END IF
X(i) = r + t ! advance seam
j1 = j1 + 1
IF (j1 > n) j1 = 1
END DO
* bring into 0...t
DO i = 1, ties
IF (AVE(i) .GE. t) AVE(i) = AVE(i) - t
END DO
dev = SNGL(smin) / FLOAT(n) ! deviation
RETURN
END ! of CYMED
*-----------------------------------------------------------------------
* imode - the mode of an array of integers
*
*___Name_______Type______In/Out____Description__________________________
* JX(N) Integer In An array
* JAVE(N) Integer Out Most frequent value in JX
* JPERM(N) Integer Neither Permutation vector
* N Integer In Length of array
* JTIES Integer Out How many ties for mode
*-----------------------------------------------------------------------
SUBROUTINE IMODE (JX, JAVE, JPERM, N, JTIES)
IMPLICIT NONE
INTEGER JX, JAVE, JPERM, N, JTIES
DIMENSION JX(N), JAVE(N), JPERM(N)
INTEGER I, J, JI, JJ, L, M
* begin
IF (N < 1) RETURN
IF (N .EQ. 1) THEN ! special case
JAVE(1) = JX(1)
JTIES = 1
RETURN
END IF
* find duplicates, mark positions and populations
CALL QSORTI (JX, JPERM, N)
M = 0
J = 1
JJ = JPERM(1)
JI = JPERM(2)
DO I = 2, N
JI = JPERM(I)
IF (JX(JI) .NE. JX(JJ)) THEN
L = I - J
IF (L > M) THEN
M = L
JTIES = 0
END IF
IF (L .EQ. M) THEN
JTIES = JTIES + 1
JAVE(JTIES) = JX(JJ)
END IF
J = I
JJ = JI
END IF
END DO
L = N - J + 1
IF (L > M) THEN
M = L
JTIES = 0
END IF
IF (L .EQ. M) THEN
JTIES = JTIES + 1
JAVE(JTIES) = JX(JI)
END IF
RETURN
END ! of IMODE
************************************************************************
* discrete mode of Real approximate integers
*
* I/O list
*_Name__________Type_________I/O________Description______________________
* X[N] Real In data array
* N Integer In length of array
* R Integer In maximum value in X
* AVE[R] Real Out modes
* TIES Integer Out how many entries in AVE[]
************************************************************************
SUBROUTINE RMODE (X, N, R, AVE, TIES)
IMPLICIT NONE
INTEGER N, R, TIES
REAL X, AVE
DIMENSION X(N), AVE(R)
* local variables
INTEGER I, J
* begin
IF (N .EQ. 1) THEN
TIES = 1
AVE(1) = X(1)
RETURN
END IF
DO J = 1, R
AVE(J) = 0.
END DO
DO I = 1, N ! count
J = NINT(X(I))
AVE(J) = AVE(J) + 1.
END DO
I = 0
TIES = 0
DO J = 1, R ! find max
IF (NINT(AVE(J)) > I) THEN
TIES = 0
I = NINT(AVE(J))
END IF
IF (NINT(AVE(J)) .EQ. I) THEN
TIES = TIES + 1
AVE(TIES) = FLOAT(J)
END IF
END DO
RETURN
END ! of RMODE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Weighted averages and weighted circular averages
! by Andy Allinger, 2016, released to the public domain
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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 mean
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE WMEAN (X, W, N, XWM)
IMPLICIT NONE
INTEGER N
REAL X, W, XWM
DIMENSION X(N), W(N)
! local variables
INTEGER I
DOUBLE PRECISION SW, SWX, C, Y, T
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! begin
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF (N < 1) RETURN
SW = 0.
C = 0.
DO I = 1, N ! compensated addition for SW
Y = W(I) - C
T = SW + Y
C = (T - SW) - Y
SW = T
END DO
SWX = 0.
C = 0.
DO I = 1, N ! compensated addition for SWX
Y = W(I) * X(I) - C
T = SWX + Y
C = (T - SWX) - Y
SWX = T
END DO
IF (SW > 0.) XWM = SNGL(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
& RMACHINE, ! 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 = T * RMACHINE(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 = RMACHINE(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) .LE. 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
RETURN
END ! of WCYMEN
*!!!!!!!!!!!!!!!!!!!!!! End of file average.f !!!!!!!!!!!!!!!!!!!!!!!!!!