¿Se puede calcular la desviación estándar para la media armónica?

12

¿Se puede calcular la desviación estándar para la media armónica? Entiendo que la desviación estándar se puede calcular para la media aritmética, pero si tiene una media armónica, ¿cómo calcula la desviación estándar o CV?

kjetil b halvorsen
fuente

Respuestas:

13

La media armónica de las variables aleatorias se define comoHX1,...,Xn

H=11ni=1n1Xi

Tomando momentos de las fracciones es un negocio sucio, así que en vez yo preferiría trabajar con el . Ahora1/H

1H=1ni=1n1Xi
.

Usando el teorema del límite central, obtenemos que

n(H1EX11)N(0,VarX11)

si, por supuesto, y son iid, ya que simplemente trabajamos con la media aritmética de las variables .VarX11<XiYi=Xi1

Ahora usando el método delta para la función obtenemos queg(x)=x1

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

Este resultado es asintótico, pero para aplicaciones simples puede ser suficiente.

Actualización Como @whuber señala con razón, las aplicaciones simples son un nombre inapropiado. El teorema del límite central se mantiene solo si existe, lo cual es una suposición bastante restrictiva.VarX11

Actualización 2 Si tiene una muestra, para calcular la desviación estándar, simplemente conecte los momentos de muestra a la fórmula. Entonces, para la muestra , la estimación de la media armónica esX1,...,Xn

H^=11ni=1n1Xi

Los momentos de muestra y respectivamente son:EX11Var(X11)

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

aquí significa recíproco.R

Finalmente, la fórmula aproximada para la desviación estándar de esH^

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

Ejecuté algunas simulaciones de Monte-Carlo para variables aleatorias distribuidas uniformemente en intervalos . Aquí está el código:[2,3]

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

Simulé Nmuestras de muestras ndimensionadas. Para cada nmuestra de tamaño, calculé la estimación de la estimación estándar (función sdhm). Luego comparo la media y la desviación estándar de estas estimaciones con la desviación estándar de la muestra de la media armónica estimada para cada muestra, que supuestamente debería ser la verdadera desviación estándar de la media armónica.

Como puede ver, los resultados son bastante buenos incluso para tamaños de muestra moderados. Por supuesto, la distribución uniforme es muy buena, por lo que no sorprende que los resultados sean buenos. Dejaré que otra persona investigue el comportamiento de otras distribuciones, el código es muy fácil de adaptar.

Nota: En la versión anterior de esta respuesta hubo un error en el resultado del método delta, variación incorrecta.

mpiktas
fuente
2
@mpiktas Este es un buen comienzo y proporciona una guía cuando el CV es bajo. Pero incluso en situaciones prácticas y simples, no está claro si se aplica el CLT. Esperaría que los recíprocos de muchas variables no tengan segundos o incluso primeros momentos finitos cuando haya una probabilidad apreciable de que sus valores puedan estar cerca de cero. También esperaría que el método delta no se aplique debido a los derivados potencialmente grandes del recíproco cercano a cero. Por lo tanto, podría ayudar a caracterizar con mayor precisión las "aplicaciones simples" donde su método podría funcionar. Por cierto, ¿qué es "D"?
whuber
@whuber, D es por varianza, . Por aplicaciones simples me refería a aquellas para las que existe varianza y media de reciprocidad. Como dice para las variables aleatorias con una probabilidad apreciable de que sus valores podrían estar cerca de cero, el recíproco puede no tener una media. Pero entonces la respuesta a la pregunta original es no. Supuse que el OP preguntó si es posible calcular la desviación estándar cuando existe. Claramente no lo hace para muchas variables aleatorias. DX=E(XEX)2
mpiktas
@whuber, BTW por curiosidad, es una notación bastante estándar para mí, pero uno podría decir que vengo de la escuela de probabilidad rusa. ¿No es tan común en el "Occidente capitalista"? :)DX
mpiktas
@mpiktas Nunca he visto esta notación de variación. ¡Mi primera reacción fue que es un operador diferencial! Las notaciones estándar son mnemotécnicas, como . DVar[X]
whuber
1
El artículo "Distribuciones invertidas" de EL Lehmann y Juliet Popper Shaffer es una lectura interesante sobre las distribuciones de variables aleatorias invertidas.
emakalic
2

Mi respuesta a una pregunta relacionada señala que la media armónica de un conjunto de datos positivos es una estimación de mínimos cuadrados ponderados (WLS) (con pesos ). Por lo tanto, puede calcular su error estándar utilizando métodos WLS. Esto tiene algunas ventajas, que incluyen simplicidad, generalidad e interpretabilidad, además de ser producido automáticamente por cualquier software estadístico que permita ponderaciones en su cálculo de regresión.xi1/xi

La desventaja principal es que el cálculo no produce buenos intervalos de confianza para distribuciones subyacentes muy sesgadas. Es probable que eso sea un problema con cualquier método de propósito general: la media armónica es sensible a la presencia de incluso un pequeño valor en el conjunto de datos.

Para ilustrar, aquí hay distribuciones empíricas de muestras generadas independientemente de tamaño partir de una distribución Gamma (5) (que está ligeramente sesgada). Las líneas azules muestran la media armónica verdadera (igual a ) mientras que las líneas discontinuas rojas muestran las estimaciones de mínimos cuadrados ponderados. Las bandas grises verticales alrededor de las líneas azules son intervalos de confianza aproximados de 95% para la media armónica. En este caso, en las muestras, el IC cubre la media armónica verdadera. Las repeticiones de esta simulación (con semillas aleatorias) sugieren que la cobertura está cerca de la tasa prevista del 95%, incluso para estos pequeños conjuntos de datos.20n=12420

Cifras

Aquí está el Rcódigo para la simulación y las figuras.

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
fuente
1

Aquí hay un ejemplo de exponencial r.v's.

La media armónica para puntos de datos se define comon

S=11ni=1nXi

Suponga que tiene iid muestras de una variable aleatoria exponencial, . La suma de variables exponenciales sigue una distribución gammanXiExp(λ)n

i=1nXiGamma(n,θ)

donde . También sabemos queθ=1λ

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

La distribución de es por lo tantoS

SInvGamma(n,nθ)

La varianza (y la desviación estándar) de este rv son bien conocidas, véase, por ejemplo, aquí .

emakalic
fuente
3
su definición de media armónica no concuerda con wikipedia
mpiktas
Usar exponenciales es un buen enfoque para comprender el problema.
whuber
1
Toda esperanza no está completamente perdida. Si Xi ~ Exp (\ lambda) entonces Xi ~ Gamma (1, \ lambda) entonces 1 / Xi ~ InvGamma (1, 1 / \ lambda). Luego use "V. Witkovsky (2001) Calculando la distribución de una combinación lineal de variables gamma invertidas, Kybernetika 37 (1), 79-90" y vea hasta dónde llega.
Tristan
0

Existe cierta preocupación de que CLT de mpiktas requiere una variación acotada en . Es cierto que tiene colas locas cuando tiene una densidad positiva alrededor de cero. Sin embargo, en muchas aplicaciones que usan la media armónica, . Aquí, está limitado por , ¡dándote todos los momentos que quieras!1/X1/XXX11/X1

karl
fuente
0

Lo que sugeriría es utilizar la siguiente fórmula como sustituto de la desviación estándar:

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

donde . Lo bueno de esta fórmula es que se minimiza cuando , y tiene las mismas unidades que la desviación estándar (que son las mismas unidades que tiene).x^=N1xix^=N1xix

Esto está en analogía con la desviación estándar, que es el valor que toma cuando se minimiza sobre . Se minimiza cuando es la media: .1N(x^xi)2x^x^x^=μ=1Nxi

Gil Wolff
fuente