# Can the standard deviation be calculated for harmonic mean?

Can the standard deviation be calculated for the harmonic mean? I understand that the standard deviation can be calculated for arithmetic mean, but if you have harmonic mean, how do you calculate the standard deviation or CV?

The harmonic mean $$H$$ of random variables $$X_1,...,X_n$$ is defined as

$$H=\frac{1}{\frac{1}{n}\sum_{i=1}^n\frac{1}{X_i}}$$

Taking moments of fractions is a messy business, so instead I would prefer working with the $$1/H$$. Now

$$\frac{1}{H}=\frac{1}{n}\sum_{i=1}^n\frac{1}{X_i}$$.

Usin central limit theorem we immediately get that

$$\sqrt{n}\left(H^{-1}-EX_1^{-1}\right)\to N(0,VarX_1^{-1})$$

if of course $$VarX_1^{-1}<\infty$$ and $$X_i$$ are iid, since we simple work with arithmetic mean of variables $$Y_i=X_i^{-1}$$.

Now using delta method for function $$g(x)=x^{-1}$$ we get that

$$\sqrt{n}(H-(EX_1^{-1})^{-1})\to N\left(0, \frac{VarX_1^{-1}}{(EX_1^{-1})^4}\right)$$

This result is asymptotic, but for simple applications it might suffice.

Update As @whuber rightfully points out, simple applications is a misnomer. The central limit theorem holds only if $$VarX_1^{-1}$$ exists, which is quite a restrictive assumption.

Update 2 If you have a sample, then to calculate the standard deviation, simply plug sample moments into the formula. So for sample $$X_1,...,X_n$$, the estimate of harmonic mean is

\begin{align} \hat{H}=\frac{1}{\frac{1}{n}\sum_{i=1}^n\frac{1}{X_i}} \end{align}

the sample moments $$EX_1^{-1}$$ and $$Var(X_1^{-1})$$ respectively are:

\begin{align} \hat{\mu}_{R}&=\frac{1}{n}\sum_{i=1}^n\frac{1}{X_i}\\\\ \hat{\sigma}_{R}^2&=\frac{1}{n}\sum_{i=1}^n\left(\frac{1}{X_i}-\mu_R\right)^2 \end{align}

here $$R$$ stands for reciprocal.

Finally the approximate formula for standard deviation of $$\hat{H}$$ is

\begin{align*} sd(\hat{H})=\sqrt{\frac{\hat{\sigma}_R^2}{n\hat{\mu}_R^4}} \end{align*}

I ran some Monte-Carlo simulations for random variables uniformly distributed in interval $$[2,3]$$. Here is the code:

hm <- function(x)1/mean(1/x)
sdhm <- function(x)sqrt((mean(1/x))^(-4)*var(1/x)/length(x))

n<-1000

nn <- c(10,30,50,100,500,1000,5000,10000)

N<-1000

mc<-foreach(n=nn,.combine=rbind) %do% {

rr <- matrix(runif(n*N,min=2,max=3),nrow=N)

c(n,mean(apply(rr,1,sdhm)),sd(apply(rr,1,sdhm)),sd(apply(rr,1,hm)))

}
colnames(mc) <- c("n","DeltaSD","sdDeltaSD","trueSD")

> mc
n     DeltaSD    sdDeltaSD      trueSD
result.1    10 0.089879211 1.528423e-02 0.091677622
result.2    30 0.052870477 4.629262e-03 0.051738941
result.3    50 0.040915607 2.705137e-03 0.040257673
result.4   100 0.029017031 1.407511e-03 0.028284458
result.5   500 0.012959582 2.750145e-04 0.013200580
result.6  1000 0.009139193 1.357630e-04 0.009115592
result.7  5000 0.004094048 2.685633e-05 0.004070593
result.8 10000 0.002894254 1.339128e-05 0.002964259


I simulated N samples of n sized sample. For each n sized sample I calculated estimate of standard estimation (function sdhm). Then I compare the mean and standard deviation of these estimates with the sample standard deviation of harmonic mean estimated for each sample, which supposably should be the true standard deviation of harmonic mean.

As you can see the results are quite good even for moderate sample sizes. Of course uniform distribution is a very well behaved one, so it is not surprising that results are good. I'll leave for someone else to investigate the behaviour for other distributions, the code is very easy to adapt.

Note: In previous version of this answer there was an error in the result of delta method, incorrect variance.

• @mpiktas This is a nice start and provides some guidance when the CV is low. But even in practical, simple situations it is not clear that the CLT applies. I would expect reciprocals of many variables not to have finite second or even first moments when there's any appreciable probability their values could be close to zero. I would also expect the delta method not to apply due to the potentially large derivatives of the reciprocal near zero. Thus it could help to more precisely characterize the "simple applications" where your method might work. BTW, what is "D"? – whuber Feb 22 '11 at 14:58
• @whuber, D is for variance, $DX=E(X-EX)^2$. By simple applications I meant the ones for which variance and mean of reciprocal exists. As you say for random variables with appreciable probability that their values could be close for zero, reciprocal may not even have mean. But then the answer to original question is no. I assumed that the OP asked whether it is possible to calculate standard deviation when it exists. It clearly does not for a lot of random variables. – mpiktas Feb 22 '11 at 15:22
• @whuber, BTW out of curiosity $DX$ is pretty standard notation for me, but one might say that I come from Russian probability school. It is not so common in "capitalistic West"? :) – mpiktas Feb 22 '11 at 15:23
• @mpiktas I have never seen this notation for variance. My first reaction was that $D$ is a differential operator! The standard notations are mnemonic, such as $Var[X]$. – whuber Feb 22 '11 at 15:31
• The paper "Inverted Distributions" by E. L. Lehmann and Juliet Popper Shaffer is an interesting read regarding distributions of inverted random variables. – emakalic Feb 22 '11 at 23:27

My answer to a related question points out that the harmonic mean of a set of positive data $x_i$ is a weighted least squares (WLS) estimate (with weights $1/x_i$). You can therefore compute its standard error using WLS methods. This has some advantages, including simplicity, generality, and interpretability as well as being automatically produced by any statistical software that allows weights in its regression calculation.

The principal disadvantage is that the calculation does not produce good confidence intervals for highly skewed underlying distributions. That's likely to be a problem with any general-purpose method: the harmonic mean is sensitive to the presence of even a single tiny value in the dataset.

To illustrate, here are empirical distributions of $20$ independently generated samples of size $n=12$ from a Gamma(5) distribution (which is modestly skewed). The blue lines show the true harmonic mean (equal to $4$) while the red dashed lines show the weighted least squares estimates. The vertical gray bands around the blue lines are approximate two-sided 95% confidence intervals for the harmonic mean. In this case, in all $20$ samples the CI covers the true harmonic mean. Repetitions of this simulation (with random seeds) suggest the coverage is close to the intended 95% rate, even for these small datasets. Here is the R code for the simulation and figures.

k <- 5             # Gamma parameter
n <- 12            # Sample size
hm <- k-1          # True harmonic mean
set.seed(17)

t.crit <- -qt(0.05/2, n-1)
par(mfrow=c(4, 5))
for(i in 1:20) {
#
# Generate a random sample.
#
x <- rgamma(n, k)
#
# Estimate the harmonic mean.
#
fit <- lm(x ~ 1, weights=1/x)
beta <- coef(summary(fit))[1, ]
message("Harmonic mean estimate is ", signif(beta["Estimate"], 3),
" +/- ", signif(beta["Std. Error"], 3))
#
# Plot the results.
#
covers <- abs(beta["Estimate"] - hm) <= t.crit*beta["Std. Error"]
plot(ecdf(x), main="Empirical CDF", sub=ifelse(covers, "", "***"))
rect(beta["Estimate"] - t.crit*beta["Std. Error"], 0,
beta["Estimate"] + t.crit*beta["Std. Error"], 1.25,
border=NA, col=gray(0.5, alpha=0.10))
abline(v = hm, col="Blue", lwd=2)
abline(v = beta["Estimate"], col="Red", lty=3, lwd=2)
}


Here is an example for Exponential r.v's.

The harmonic mean for $n$ data points is defined as

$$S = \frac{1}{\frac{1}{n} \sum_{i=1}^n X_i}$$

Suppose you have $n$ iid samples of an Exponential random variable, $X_i \sim {\rm Exp}(\lambda)$. The sum of $n$ Exponential variables follows a Gamma distribution

$$\sum_{i=1}^n X_i \sim {\rm Gamma}(n, \theta)$$

where $\theta = \frac{1}{\lambda}$. We also know that

$$\frac{1}{n} {\rm Gamma}(n, \theta) \sim {\rm Gamma}(n, \frac{\theta}{n})$$

The distribution of $S$ is therefore

$$S \sim {\rm InvGamma}(n, \frac{n}{\theta})$$

The variance (and standard deviation) of this r.v. are well known, see, for example here.

• your definition for harmonic mean does not agree with wikipedia – mpiktas Feb 22 '11 at 7:20
• Using exponentials is a good approach to understanding the problem. – whuber Feb 23 '11 at 16:28
• All hope is not entirely lost. If Xi ~ Exp(\lambda) then Xi ~ Gamma(1, \lambda) so 1/Xi ~ InvGamma(1, 1/\lambda). Then use "V. Witkovsky (2001) Computing the distribution of a linear combination of inverted gamma variables, Kybernetika 37(1), 79-90" and see how far you get! – tristan Mar 20 '15 at 7:36

There is some concern that mpiktas's CLT requires a bounded variance on $1/X$. It is true that $1/X$ has crazy tails when $X$ has positive density around zero. However, in many applications using the harmonic mean, $X\ge1$. Here, $1/X$ is bounded by $1$, giving you all the moments that you want!

What I would suggest is to use the following formula as a substitute for the standard deviation: $$\sigma=\sqrt{\frac{N}{\sum_{i=1}^{N}{\left(\frac{1}{\hat{x}}-\frac{1}{x_i}\right)^2}}},$$

where $\hat{x} = \frac{N}{\sum \frac{1}{x_i}}$. The nice thing about this formula is that it is minimized when $\hat{x} = \frac{N}{\sum \frac{1}{x_i}}$, and it has the same units as the standard deviation would (which are the same units as $x$ has).

This is in analogy to the standard deviation, which is the value that $\sqrt{\frac{1}{N}\sum(\hat{x}-x_i)^2}$ takes when it is minimized over $\hat{x}$. It is minimized when $\hat{x}$ is the mean: $\hat{x}=\mu=\frac{1}{N}\sum x_i$.