Dapatkah standar deviasi dihitung untuk rata-rata harmonik?

12

Bisakah standar deviasi dihitung untuk rata-rata harmonik? Saya mengerti bahwa standar deviasi dapat dihitung untuk rata-rata aritmatika, tetapi jika Anda memiliki rata-rata harmonis, bagaimana Anda menghitung standar deviasi atau CV?

kjetil b halvorsen
sumber

Jawaban:

13

Rata-rata harmonik dari variabel acak didefinisikan sebagaiHX1,...,Xn

H=11ni=1n1Xi

Mengambil momen pecahan adalah bisnis berantakan, jadi bukan saya lebih suka bekerja dengan . Sekarang1/H

1H=1ni=1n1Xi
.

Dalam teorema limit pusat, kita segera mendapatkan itu

n(H1EX11)N(0,VarX11)

jika tentu saja dan adalah iid, karena kami bekerja sederhana dengan rata-rata aritmatika variabel .VarX11<XiYi=Xi1

Sekarang menggunakan metode delta untuk fungsi kita dapatkan itug(x)=x1

n(H-(EX1-1)-1)N(0,VSebuahrX1-1(EX1-1)4)

Hasil ini asimptotik, tetapi untuk aplikasi sederhana mungkin cukup.

Pembaruan Seperti yang ditunjukkan @whuber dengan tepat, aplikasi sederhana adalah nama yang salah. Teorema batas pusat hanya berlaku jika ada, yang merupakan asumsi yang sangat membatasi.VSebuahrX1-1

Perbarui 2 Jika Anda memiliki sampel, maka untuk menghitung standar deviasi, cukup masukkan momen sampel ke dalam rumus. Jadi untuk sampel , estimasi rata-rata harmonik adalahX1,...,Xn

H^=11nsaya=1n1Xsaya

contoh momen dan masing-masing adalah:EX1-1VSebuahr(X1-1)

μ^R=1nsaya=1n1Xsayaσ^R2=1nsaya=1n(1Xsaya-μR)2

di sini berarti timbal balik.R

Akhirnya rumus perkiraan untuk standar deviasi adalahH^

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

Saya menjalankan beberapa simulasi Monte-Carlo untuk variabel acak yang terdistribusi secara merata dalam interval . Ini kodenya:[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

Saya mensimulasikan Nsampel sampel nberukuran. Untuk setiap nsampel berukuran saya menghitung estimasi estimasi standar (fungsi sdhm). Kemudian saya membandingkan rata-rata dan standar deviasi dari estimasi ini dengan standar deviasi sampel dari rata-rata harmonik yang diperkirakan untuk setiap sampel, yang seharusnya merupakan standar deviasi rata-rata rata-rata harmonik.

Seperti yang Anda lihat hasilnya cukup baik bahkan untuk ukuran sampel sedang. Tentu saja distribusi seragam adalah yang berperilaku sangat baik, sehingga tidak mengherankan bahwa hasilnya baik. Saya akan pergi ke orang lain untuk menyelidiki perilaku untuk distribusi lain, kode ini sangat mudah untuk beradaptasi.

Catatan: Di versi sebelumnya dari jawaban ini ada kesalahan dalam hasil metode delta, varians salah.

mpiktas
sumber
2
@mpiktas Ini adalah awal yang baik dan memberikan beberapa panduan ketika CV rendah. Tetapi bahkan dalam situasi praktis dan sederhana tidak jelas bahwa CLT berlaku. Saya berharap timbal balik dari banyak variabel tidak memiliki momen hingga saat kedua atau bahkan saat ketika ada kemungkinan yang cukup nilai mereka bisa mendekati nol. Saya juga berharap metode delta tidak berlaku karena turunan resiprokal yang berpotensi besar mendekati nol. Dengan demikian dapat membantu untuk lebih tepatnya menandai "aplikasi sederhana" tempat metode Anda bekerja. BTW, apa itu "D"?
whuber
DX=E(X-EX)2
DX
DVSebuahr[X]
1
Makalah "Distribusi Terbalik" oleh EL Lehmann dan Juliet Popper Shaffer adalah bacaan yang menarik tentang distribusi variabel acak terbalik.
emakalic
2

xsaya1/xsaya

Kerugian utama adalah bahwa perhitungan tidak menghasilkan interval kepercayaan yang baik untuk distribusi dasar yang sangat miring. Itu kemungkinan akan menjadi masalah dengan metode tujuan umum: rata-rata harmonik peka terhadap keberadaan bahkan nilai kecil tunggal dalam dataset.

20n=12420

Tokoh

Berikut adalah Rkode untuk simulasi dan gambar.

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

Berikut adalah contoh untuk Exponential r.v's.

n

S=11nsaya=1nXsaya

nXsayaExhal(λ)n

saya=1nXsayaGSebuahmmSebuah(n,θ)

θ=1λ

1nGSebuahmmSebuah(n,θ)GSebuahmmSebuah(n,θn)

S

SsayanvGSebuahmmSebuah(n,nθ)

Varian (dan standar deviasi) rv ini sudah dikenal, lihat, misalnya di sini .

emakalic
sumber
3
definisi Anda untuk rata-rata harmonis tidak setuju dengan wikipedia
mpiktas
Menggunakan eksponensial adalah pendekatan yang baik untuk memahami masalah.
whuber
1
Semua harapan tidak sepenuhnya hilang. Jika Xi ~ Exp (\ lambda) maka Xi ~ Gamma (1, \ lambda) jadi 1 / Xi ~ InvGamma (1, 1 / \ lambda). Kemudian gunakan "V. Witkovsky (2001) Menghitung distribusi kombinasi linear dari variabel gamma terbalik, Kybernetika 37 (1), 79-90" dan lihat seberapa jauh Anda dapatkan!
tristan
0

1/X1/XXX11/X1

karl
sumber
0

σ=Nsaya=1N(1x^-1xsaya)2,

x^=N1xsayax^=N1xsayax

1N(x^-xsaya)2x^x^x^=μ=1Nxsaya

Gil Wolff
sumber