Kann die Standardabweichung für das harmonische Mittel berechnet werden?

12

Kann die Standardabweichung für das harmonische Mittel berechnet werden? Ich verstehe, dass die Standardabweichung für das arithmetische Mittel berechnet werden kann, aber wenn Sie ein harmonisches Mittel haben, wie berechnen Sie die Standardabweichung oder den CV?

kjetil b halvorsen
quelle

Antworten:

13

Das harmonische Mittel H der Zufallsvariablen X1,...,Xn ist definiert als

H=11ni=1n1Xi

Momente von Brüchen zu nehmen ist eine unordentliche Angelegenheit, daher würde ich lieber mit dem 1/H . Jetzt

1H=1ni=1n1Xi
.

Mit dem zentralen Grenzwertsatz bekommen wir das sofort

n(H1EX11)N(0,VarX.1- -1)

wenn natürlich VarX11< und Xi iid sind, da wir einfach mit dem arithmetischen Mittel der Variablen Yi=Xi1 .

Mit der Delta-Methode für die Funktion g(x)=x1 wir das

n(H(EX11)1)N(0,VarX11(EX11)4)

Dieses Ergebnis ist asymptotisch, aber für einfache Anwendungen kann es ausreichen.

Update Wie @whuber zu Recht betont, sind einfache Anwendungen eine Fehlbezeichnung. Der zentrale Grenzwertsatz gilt nur, wenn VarX11 existiert, was eine recht restriktive Annahme ist.

Update 2 Wenn Sie eine Stichprobe haben, fügen Sie zur Berechnung der Standardabweichung einfach Stichprobenmomente in die Formel ein. Also für Probe X1,...,Xn ist die Schätzung des harmonischen Mittelwerts

H^=11ni=1n1Xi

Die Abtastmomente EX11 und Var(X11) sind:

μ^R=1ni=1n1Xiσ^R2=1ni=1n(1XiμR)2

hier steht R für wechselseitig.

Schließlich wird die Näherungsformel für die Standardabweichung von H ist ,H^

sd(H^)=σ^R2nμ^R4

Ich habe einige Monte-Carlo-Simulationen für Zufallsvariablen durchgeführt, die gleichmäßig im Intervall verteilt sind [2,3] . Hier ist der 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

Ich habe NProben von ngroßen Proben simuliert . Für jede nStichprobe berechnete ich die Schätzung der Standardschätzung (Funktionsdhm ). Dann vergleiche ich den Mittelwert und die Standardabweichung dieser Schätzungen mit der für jede Probe geschätzten Standardabweichung des harmonischen Mittelwerts, die vermutlich die wahre Standardabweichung des harmonischen Mittelwerts sein sollte.

Wie Sie sehen können, sind die Ergebnisse auch bei moderaten Stichprobengrößen recht gut. Natürlich ist eine gleichmäßige Verteilung sehr gut, daher ist es nicht verwunderlich, dass die Ergebnisse gut sind. Ich überlasse es jemand anderem, das Verhalten für andere Distributionen zu untersuchen. Der Code ist sehr einfach anzupassen.

Hinweis: In der vorherigen Version dieser Antwort ist ein Fehler im Ergebnis der Delta-Methode aufgetreten, falsche Varianz.

mpiktas
quelle
2
@mpiktas Dies ist ein guter Anfang und bietet eine Anleitung, wenn der Lebenslauf niedrig ist. Aber auch in praktischen, einfachen Situationen ist nicht klar, dass das CLT gilt. Ich würde erwarten, dass Kehrwerte vieler Variablen keine endlichen zweiten oder sogar ersten Momente haben, wenn es eine nennenswerte Wahrscheinlichkeit gibt, dass ihre Werte nahe Null sein könnten. Ich würde auch erwarten, dass die Delta-Methode aufgrund der potenziell großen Ableitungen des Kehrwerts nahe Null nicht angewendet wird. Daher könnte es hilfreich sein, die "einfachen Anwendungen", in denen Ihre Methode möglicherweise funktioniert, genauer zu charakterisieren. Übrigens, was ist "D"?
whuber
DX=E(XEX)2
@whuber, BTW out of curiosity D.X. 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
@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
1
The paper "Inverted Distributions" by E. L. Lehmann and Juliet Popper Shaffer is an interesting read regarding distributions of inverted random variables.
emakalic
2

My answer to a related question points out that the harmonic mean of a set of positive data xi is a weighted least squares (WLS) estimate (with weights 1/xi). 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=12aus einer Gamma (5) -Verteilung (die leicht verzerrt ist). Die blauen Linien zeigen den wahren harmonischen Mittelwert (gleich4), während die roten gestrichelten Linien die Schätzungen der gewichteten kleinsten Quadrate zeigen. Die vertikalen grauen Bänder um die blauen Linien sind ungefähre zweiseitige 95% -Konfidenzintervalle für das harmonische Mittel. In diesem Fall insgesamt20Proben Der CI deckt den wahren harmonischen Mittelwert ab. Wiederholungen dieser Simulation (mit zufälligen Samen) legen nahe, dass die Abdeckung selbst für diese kleinen Datensätze nahe an der beabsichtigten 95% -Rate liegt.

Figures

Hier ist der RCode für die Simulation und die Abbildungen.

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)
}
whuber
quelle
1

Hier ist ein Beispiel für Exponential r.v.

Das harmonische Mittel für n Datenpunkte ist definiert als

S.=11nich=1nX.ich

Angenommen, Sie haben n iid Stichproben einer exponentiellen Zufallsvariablen, X.ichE.xp((λ). The sum of n Exponential variables follows a Gamma distribution

i=1nXiGamma(n,θ)

where θ=1λ. We also know that

1nGamma(n,θ)Gamma(n,θn)

The distribution of S is therefore

SInvGamma(n,nθ)

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

emakalic
quelle
3
your definition for harmonic mean does not agree with wikipedia
mpiktas
Using exponentials is a good approach to understanding the problem.
whuber
1
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
0

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, X1. Here, 1/X is bounded by 1, giving you all the moments that you want!

karl
quelle
0

What I would suggest is to use the following formula as a substitute for the standard deviation:

σ=Ni=1N(1x^1xi)2,

where x^=N1xi. The nice thing about this formula is that it is minimized when x^=N1xi, 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 1N(x^xi)2 takes when it is minimized over x^. It is minimized when x^ is the mean: x^=μ=1Nxi.

Gil Wolff
quelle