MODULE Update_EWMA ! Update exponentially weighted moving averages (EWMA) & variances. ! New average = lambda.x + (1 - lambda).(old average) ! The variance of the estimated mean = [lambda/(2 - lambda)].variance(x) ! Code by Alan Miller ! amiller@bigpond.net.au ! Web sites: ! http://www.ozemail.com.au/~milleraj/ ! http://users.bigpond.net.au/amiller ! Latest revision - 24 October 2001 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), SAVE, PRIVATE :: lambda, one_lambda, deg_freedom LOGICAL, SAVE, PRIVATE :: initialized = .FALSE. CONTAINS SUBROUTINE initialize(l) REAL (dp), INTENT(IN) :: l lambda = l initialized = .TRUE. one_lambda = 1.0_dp - lambda deg_freedom = 2.0_dp * one_lambda / (lambda*(2.0_dp - lambda)) RETURN END SUBROUTINE initialize SUBROUTINE update(x, average, sumsq) ! New average = lambda.x + (1 - lambda).(old average) REAL (dp), INTENT(IN) :: x REAL (dp), INTENT(IN OUT) :: average, sumsq ! Local variable REAL (dp) :: dev ! For speed, the following lines are not used ! IF (.NOT. initialized) THEN ! WRITE(*, *) ' *** EWMA not initialized ***' ! STOP ! END IF dev = x - average average = average + lambda * dev sumsq = one_lambda * sumsq + dev * (x - average) RETURN END SUBROUTINE update FUNCTION variance(sumsq) RESULT(fn_val) ! Calculate the variance from the sum of squares. REAL (dp), INTENT(IN) :: sumsq REAL (dp) :: fn_val fn_val = sumsq / deg_freedom RETURN END FUNCTION variance END MODULE Update_EWMA PROGRAM test_EWMA ! Generate random normal values with mean = 7.0 and variance = 1.0 ! Calculate and output the weighted average and variance after every 10 values. USE Update_EWMA IMPLICIT NONE REAL (dp) :: lambda, x, average, sumsq INTEGER :: i WRITE(*, '(a)', ADVANCE='NO') ' Enter lambda between 0 and 1: ' READ(*, *) lambda CALL initialize(lambda) ! Set arbitrary starting values for average & sumsq average = 5.0_dp sumsq = 1.0_dp WRITE(*, '(a, f6.2)') ' Lambda = ', lambda WRITE(*, *) ' I Avge. Variance' DO i = 1, 1000 x = rnorm() + 7.0 CALL update(x, average, sumsq) IF (i == 10*(i/10)) THEN WRITE(*, '(i5, f7.1, f7.1)') i, average, variance(sumsq) END IF END DO STOP CONTAINS FUNCTION rnorm() RESULT( fn_val ) ! Generate a random normal deviate using the polar method. ! Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating ! normal variables', Siam Rev., vol.6, 260-264, 1964. IMPLICIT NONE REAL :: fn_val ! Local variables REAL :: u, sum REAL, SAVE :: v, sln LOGICAL, SAVE :: second = .FALSE. REAL, PARAMETER :: one = 1.0, vsmall = TINY( one ) IF (second) THEN ! If second, use the second random number generated on last call second = .false. fn_val = v*sln ELSE ! First call; generate a pair of random normals second = .true. DO CALL RANDOM_NUMBER( u ) CALL RANDOM_NUMBER( v ) u = SCALE( u, 1 ) - one v = SCALE( v, 1 ) - one sum = u*u + v*v + vsmall ! vsmall added to prevent LOG(zero) / zero IF(sum < one) EXIT END DO sln = SQRT(- SCALE( LOG(sum), 1 ) / sum) fn_val = u*sln END IF RETURN END FUNCTION rnorm END PROGRAM test_EWMA