exponential-smoothingpythonstandard deviation
I wrote a simple function in Python to calculate the exponentially weighted mean:
def test():
x = [1,2,3,4,5]
alpha = 0.98
s_old = x[0]
for i in range(1, len(x)):
s = alpha * x[i] + (1- alpha) * s_old
s_old = s
return s
However, how can I calculate the corresponding SD?
This reply presents two solutions: Sheppard's corrections and a maximum likelihood estimate. Both closely agree on an estimate of the standard deviation: $7.70$ for the first and $7.69$ for the second (when adjusted to be comparable to the usual "unbiased" estimator).
"Sheppard's corrections" are formulas that adjust moments computed from binned data (like these) where
the data are assumed to be governed by a distribution supported on a finite interval $[a,b]$
that interval is divided sequentially into equal bins of common width $h$ that is relatively small (no bin contains a large proportion of all the data)
the distribution has a continuous density function.
They are derived from the Euler-Maclaurin sum formula, which approximates integrals in terms of linear combinations of values of the integrand at regularly spaced points, and therefore generally applicable (and not just to Normal distributions).
Although strictly speaking a Normal distribution is not supported on a finite interval, to an extremely close approximation it is. Essentially all its probability is contained within seven standard deviations of the mean. Therefore Sheppard's corrections are applicable to data assumed to come from a Normal distribution.
The first two Sheppard's corrections are
Use the mean of the binned data for the mean of the data (that is, no correction is needed for the mean).
Subtract $h^2/12$ from the variance of the binned data to obtain the (approximate) variance of the data.
Where does $h^2/12$ come from? This equals the variance of a uniform variate distributed over an interval of length $h$. Intuitively, then, Sheppard's correction for the second moment suggests that binning the data--effectively replacing them by the midpoint of each bin--appears to add an approximately uniformly distributed value ranging between $-h/2$ and $h/2$, whence it inflates the variance by $h^2/12$.
Let's do the calculations. I use R to illustrate them, beginning by specifying the counts and the bins:
counts <- c(1,2,3,4,1)
bin.lower <- c(40, 45, 50, 55, 70)
bin.upper <- c(45, 50, 55, 60, 75)
The proper formula to use for the counts comes from replicating the bin widths by the amounts given by the counts; that is, the binned data are equivalent to
42.5, 47.5, 47.5, 52.5, 52.5, 57.5, 57.5, 57.5, 57.5, 72.5
Their number, mean, and variance can be directly computed without having to expand the data in this way, though: when a bin has midpoint $x$ and a count of $k$, then its contribution to the sum of squares is $kx^2$. This leads to the second of the Wikipedia formulas cited in the question.
bin.mid <- (bin.upper + bin.lower)/2
n <- sum(counts)
mu <- sum(bin.mid * counts) / n
sigma2 <- (sum(bin.mid^2 * counts) - n * mu^2) / (n-1)
The mean (mu) is $1195/22 \approx 54.32$ (needing no correction) and the variance (sigma2) is $675/11 \approx 61.36$. (Its square root is $7.83$ as stated in the question.) Because the common bin width is $h=5$, we subtract $h^2/12 = 25/12 \approx 2.08$ from the variance and take its square root, obtaining $\sqrt{675/11 - 5^2/12} \approx 7.70$ for the standard deviation.
An alternative method is to apply a maximum likelihood estimate. When the assumed underlying distribution has a distribution function $F_\theta$ (depending on parameters $\theta$ to be estimated) and the bin $(x_0, x_1]$ contains $k$ values out of a set of independent, identically distributed values from $F_\theta$, then the (additive) contribution to the log likelihood of this bin is
$$\log \prod_{i=1}^k \left(F_\theta(x_1) - F_\theta(x_0)\right) = k\log\left(F_\theta(x_1) - F_\theta(x_0)\right)$$
(see MLE/Likelihood of lognormally distributed interval).
Summing over all bins gives the log likelihood $\Lambda(\theta)$ for the dataset. As usual, we find an estimate $\hat\theta$ which minimizes $-\Lambda(\theta)$. This requires numerical optimization and that is expedited by supplying good starting values for $\theta$. The following R code does the work for a Normal distribution:
sigma <- sqrt(sigma2) # Crude starting estimate for the SD
likelihood.log <- function(theta, counts, bin.lower, bin.upper) {
mu <- theta[1]; sigma <- theta[2]
-sum(sapply(1:length(counts), function(i) {
counts[i] *
log(pnorm(bin.upper[i], mu, sigma) - pnorm(bin.lower[i], mu, sigma))
}))
}
coefficients <- optim(c(mu, sigma), function(theta)
likelihood.log(theta, counts, bin.lower, bin.upper))$par
The resulting coefficients are $(\hat\mu, \hat\sigma) = (54.32, 7.33)$.
Remember, though, that for Normal distributions the maximum likelihood estimate of $\sigma$ (when the data are given exactly and not binned) is the population SD of the data, not the more conventional "bias corrected" estimate in which the variance is multiplied by $n/(n-1)$. Let us then (for comparison) correct the MLE of $\sigma$, finding $\sqrt{n/(n-1)} \hat\sigma = \sqrt{11/10}\times 7.33 = 7.69$. This compares favorably with the result of Sheppard's correction, which was $7.70$.
To visualize these results we can plot the fitted Normal density over a histogram:
hist(unlist(mapply(function(x,y) rep(x,y), bin.mid, counts)),
breaks = breaks, xlab="Values", main="Data and Normal Fit")
curve(dnorm(x, coefficients[1], coefficients[2]),
from=min(bin.lower), to=max(bin.upper),
add=TRUE, col="Blue", lwd=2)
To some this might not look like a good fit. However, because the dataset is small (only $11$ values), surprisingly large deviations between the distribution of the observations and the true underlying distribution can occur.
Let's more formally check the assumption (made by the MLE) that the data are governed by a Normal distribution. An approximate goodness of fit test can be obtained from a $\chi^2$ test: the estimated parameters indicate the expected amount of data in each bin; the $\chi^2$ statistic compares the observed counts to the expected counts. Here is a test in R:
breaks <- sort(unique(c(bin.lower, bin.upper)))
fit <- mapply(function(l, u) exp(-likelihood.log(coefficients, 1, l, u)),
c(-Inf, breaks), c(breaks, Inf))
observed <- sapply(breaks[-length(breaks)], function(x) sum((counts)[bin.lower <= x])) -
sapply(breaks[-1], function(x) sum((counts)[bin.upper < x]))
chisq.test(c(0, observed, 0), p=fit, simulate.p.value=TRUE)
The output is
Chi-squared test for given probabilities with simulated p-value (based on 2000 replicates)
data: c(0, observed, 0)
X-squared = 7.9581, df = NA, p-value = 0.2449
The software has performed a permutation test (which is needed because the test statistic does not follow a chi-squared distribution exactly: see my analysis at How to Understand Degrees of Freedom). Its p-value of $0.245$, which is not small, shows very little evidence of departure from normality: we have reason to trust the maximum likelihood results.
You're confusing addition of random variables with concatenation of samples (easy to do, it took me a while to realize why your code didn't work!). So for independent random variables $X$ and $Y$ you can write $\rm{Var}(aX+bY) = a^2\rm{Var}(X) + b^2\rm{Var}(Y)$, noting the coefficients are squared. This would apply (approximately) to samples as follows
set.seed(1)
n <- 100
mu = 1
sigma = 1
a<-rnorm(n,mu,sigma)
b<-rnorm(n,mu,sigma)
weight_a <- 1/4
weight_b <- 3/4
sqrt(sd(a)^2*weight_a^2+sd(b)^2*weight_b^2)
[1] 0.7526849
sd(weight_a*a + weight_b*b)
[1] 0.7524718
For concatenation of samples you need to apply the formula for variance $\rm{Var}(X) = E(X^2) - (E(X))^2$, adjusted appropriately for the fact that you're using samples, i.e. multiply by $n/(n-1)$. This is illustrated with R code below.
suma2 = sum(a^2)
sumb2 = sum(b^2)
suma = sum(a)
sumb = sum(b)
sumc2 = suma2 + sumb2
sumc = suma + sumb
nc = n + n
(sumb2/n - (sumb/n)^2) * n/(n-1)
[1] 0.9175323
sd(b)^2
[1] 0.9175323
sumc2 = suma2 + sumb2
sumc = suma + sumb
(sumc2/nc - (sumc/nc)^2) * nc/(nc-1)
[1] 0.8632217
sd(c(a,b))^2
[1] 0.8632217
When you are concatenating weighted series, the weights can be incorporated naturally to give the formula.
$$ \rm{Var}(\mathbf{w}) = \left(\frac{w_a^2 S_{aa} + w_b^2 S_{bb}}{n_a + n_b} - \left(\frac{w_a S_a + w_b S_b}{n_a + n_b}\right)^2 \right) \times \frac{n_a + n_b}{n_a + n_b -1}, $$
where $\mathbf{w}$ is the series formed by concatenating series $\mathbf{a}$ of length $n_a$ weighted by $w_a$ and series $\mathbf{b}$ of length $n_b$ weighted by $w_b$. $S_{aa}$ is the sum of squares of series $\mathbf{a}$, $S_{a}$ is the sum of series $\mathbf{a}$, and likewise for $\mathbf{b}$ .
This is illustrated in the following R code.
sumw2 = weight_a^2*suma2 + weight_b^2*sumb2
sumw = weight_a*suma + weight_b*sumb
(sumw2/nc - (sumw/nc)^2) * nc/(nc-1)
[1] 0.3314697
sd(c(weight_a*a, weight_b*b))^2
[1] 0.3314697
Note that in the first case, addition, I had to use the same size of series. For convenience, I carried on with the same series. But the code used for the second case, concatenation, should work for different-sized series.
However, as @Peter Ellis states in his answer, the problem you are ultimately trying to solve may be something like a two-sample t test, where you assume that the two samples are from the same population, and so the underlying variance of each sample is the same. In this case, you want to estimate the population variance. The formula for the pooled variance is well-known and can be found on wikipedia (wikipedia also provides a generalization to multiple samples).
Best Answer
You can use the following recurrent formula:
$\sigma_i^2 = S_i = (1 - \alpha) (S_{i-1} + \alpha (x_i - \mu_{i-1})^2)$
Here $x_i$ is your observation in the $i$-th step, $\mu_{i-1}$ is the estimated EWM, and $S_{i-1}$ is the previous estimate of the variance. See Section 9 here for the proof and pseudo-code.