@@ -661,9 +661,9 @@ cdef inline void add_var(double val, double *nobs, double *mean_x,
661
661
if val == val:
662
662
nobs[0 ] = nobs[0 ] + 1
663
663
664
- delta = ( val - mean_x[0 ])
664
+ delta = val - mean_x[0 ]
665
665
mean_x[0 ] = mean_x[0 ] + delta / nobs[0 ]
666
- ssqdm_x[0 ] = ssqdm_x[0 ] + delta * (val - mean_x [0 ])
666
+ ssqdm_x[0 ] = ssqdm_x[0 ] + ((nobs[ 0 ] - 1 ) * delta ** 2 ) / nobs [0 ]
667
667
668
668
669
669
cdef inline void remove_var(double val, double * nobs, double * mean_x,
@@ -675,9 +675,9 @@ cdef inline void remove_var(double val, double *nobs, double *mean_x,
675
675
if val == val:
676
676
nobs[0 ] = nobs[0 ] - 1
677
677
if nobs[0 ]:
678
- delta = ( val - mean_x[0 ])
678
+ delta = val - mean_x[0 ]
679
679
mean_x[0 ] = mean_x[0 ] - delta / nobs[0 ]
680
- ssqdm_x[0 ] = ssqdm_x[0 ] - delta * (val - mean_x [0 ])
680
+ ssqdm_x[0 ] = ssqdm_x[0 ] - ((nobs[ 0 ] + 1 ) * delta ** 2 ) / nobs [0 ]
681
681
else :
682
682
mean_x[0 ] = 0
683
683
ssqdm_x[0 ] = 0
@@ -689,7 +689,7 @@ def roll_var(ndarray[double_t] input, int64_t win, int64_t minp,
689
689
Numerically stable implementation using Welford's method.
690
690
"""
691
691
cdef:
692
- double val, prev, mean_x = 0 , ssqdm_x = 0 , nobs = 0 , delta
692
+ double val, prev, mean_x = 0 , ssqdm_x = 0 , nobs = 0 , delta, mean_x_old
693
693
int64_t s, e
694
694
bint is_variable
695
695
Py_ssize_t i, j, N
@@ -760,10 +760,12 @@ def roll_var(ndarray[double_t] input, int64_t win, int64_t minp,
760
760
761
761
# Adding one observation and removing another one
762
762
delta = val - prev
763
- prev -= mean_x
763
+ mean_x_old = mean_x
764
+
764
765
mean_x += delta / nobs
765
- val -= mean_x
766
- ssqdm_x += (val + prev) * delta
766
+ ssqdm_x += ((nobs - 1 ) * val
767
+ + (nobs + 1 ) * prev
768
+ - 2 * nobs * mean_x_old) * delta / nobs
767
769
768
770
else :
769
771
add_var(val, & nobs, & mean_x, & ssqdm_x)
0 commit comments