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