Standar deviasi pengamatan binned

24

Saya memiliki dataset pengamatan sampel, disimpan sebagai jumlah dalam jangkauan sampah. misalnya:

min/max  count
40/44    1
45/49    2
50/54    3
55/59    4
70/74    1

Sekarang, menemukan perkiraan rata-rata dari ini cukup mudah. Cukup gunakan rata-rata (atau median) dari masing-masing rentang bin sebagai pengamatan dan menghitung sebagai bobot dan temukan rata-rata tertimbang:

x¯=1i=1Nwii=1Nwixi

Untuk kasus pengujian saya, ini memberi saya 53,82.

Pertanyaan saya sekarang adalah, apa metode yang benar untuk menemukan standar deviasi (atau varian)?

Melalui pencarian saya, saya telah menemukan beberapa jawaban, tetapi saya tidak yakin yang mana, jika ada, sebenarnya sesuai untuk dataset saya. Saya dapat menemukan rumus berikut pada pertanyaan lain di sini dan dokumen NIST acak .

s2=i=1Nwi(xix¯)2(M1)Mi=1Nwi

Yang memberikan standar deviasi 8,35 untuk test case saya. Namun, artikel Wikipedia tentang cara tertimbang memberikan kedua rumus:

s2=i=1Nwi(i=1Nwi)2i=1Nwi2i=1Nwi(xix¯)2

dan

s2=1(i=1Nwi)1i=1Nwi(xix¯)2

Yang memberikan standar deviasi 8,66 dan 7,83, masing-masing, untuk kasus uji saya.

Memperbarui

Terima kasih kepada @whuber yang menyarankan untuk melihat Koreksi Sheppard, dan komentar Anda yang bermanfaat terkait dengannya. Sayangnya, saya mengalami kesulitan memahami sumber daya yang dapat saya temukan (dan saya tidak dapat menemukan contoh yang bagus). Singkatnya, saya mengerti bahwa berikut ini adalah estimasi varians yang bias:

s2=1i=1Nwii=1Nwi(xix¯)2

Saya juga mengerti bahwa sebagian besar koreksi standar untuk bias adalah untuk sampel acak langsung dari distribusi normal. Karena itu, saya melihat dua masalah potensial bagi saya:

  1. Ini adalah sampel acak binned (yang, saya cukup yakin, adalah tempat Sheppard's Corrections masuk.)
  2. Tidak diketahui apakah data tersebut untuk distribusi normal (jadi saya berasumsi tidak, yang, saya yakin, membatalkan koreksi Sheppard.)

Jadi, pertanyaan saya yang diperbarui adalah; Apa metode yang tepat untuk menangani bias yang diberlakukan oleh rumus deviasi / varian standar tertimbang "sederhana" pada distribusi yang tidak normal? Paling khusus berkaitan dengan data yang dibuang.

Catatan: Saya menggunakan ketentuan berikut:

  • adalah varian tertimbang s2
  • adalah jumlah pengamatan. (yaitu jumlah tempat sampah)N
  • adalah jumlah bobot bukan nol. (yaitu jumlah nampan dengan jumlah)M
  • adalah bobot (yaitu jumlah)wi
  • adalah pengamatan. (Yaitu tempat sampah)xi
  • adalah rata-rata tertimbang.x¯
chezy525
sumber
5
Google "Koreksi Sheppard" untuk solusi standar untuk masalah ini.
whuber
@whuber, saya khawatir google-foo saya gagal ... Saya tidak menemukan banyak tentang cara menggunakan koreksi Sheppard. Sejauh yang saya tahu, ini adalah koreksi untuk sifat binned dari data, dan dalam kasus pengujian saya akan digunakan seperti , di manacadalah ukuran nampan (dalam test case saya, 4). Apakah ini benar? Bagaimanapun, apa yang saya temukan masih nampaknya tidak membantu saya dalam menghitungs2. s2c212cs2
chezy525
The hit kedua di saya pencarian Google memberikan rumus eksplisit (persamaan 9).
whuber
@whuber, sudah beberapa bulan, dan saya sudah mencoba membaca dokumen yang Anda tautkan beberapa kali. Saya pikir saya masih kehilangan sesuatu, tetapi yang terbaik yang saya temukan adalah bahwa persamaan terakhir yang saya cantumkan benar sebagai estimator yang tidak bias. Apakah ini benar?
chezy525
Koreksi Sheppard tidak menganggap normal.
Glen_b -Reinstate Monica

Jawaban:

18

Balasan ini menyajikan dua solusi: Koreksi Sheppard dan estimasi kemungkinan maksimum. Keduanya sangat setuju pada estimasi standar deviasi: untuk yang pertama dan 7,69 untuk yang kedua (bila disesuaikan agar sebanding dengan penduga "tidak bias" yang biasa).7.707.69


Koreksi Sheppard

"Koreksi Sheppard" adalah rumus yang mengatur momen yang dihitung dari data yang dibuang (seperti ini) di mana

  • data diasumsikan diatur oleh distribusi yang didukung pada interval terbatas [a,b]

  • interval tersebut dibagi secara berurutan menjadi nampan yang sama dengan lebar umum yang relatif kecil (tidak ada nampan berisi sebagian besar dari semua data)h

  • distribusi memiliki fungsi kepadatan kontinu.

Mereka diturunkan dari rumus jumlah Euler-Maclaurin, yang mendekati integral dalam hal kombinasi linear dari nilai-nilai integrand pada titik-titik yang berjarak secara teratur, dan oleh karena itu umumnya berlaku (dan bukan hanya untuk distribusi Normal).

Meskipun secara tegas distribusi Normal tidak didukung pada interval terbatas, hingga perkiraan yang sangat dekat. Pada dasarnya semua probabilitasnya terkandung dalam tujuh standar deviasi rata-rata. Oleh karena itu koreksi Sheppard berlaku untuk data yang diasumsikan berasal dari distribusi Normal.

Dua koreksi Sheppard pertama adalah

  1. Gunakan rata-rata dari data yang dibuang untuk rata-rata data (yaitu, tidak diperlukan koreksi untuk rata-rata).

  2. Kurangi h2/12 dari varians dari data binned untuk mendapatkan (perkiraan) varians dari data.

h2/12hh/2h/2h2/12 .

Mari kita lakukan perhitungan. Saya gunakan Runtuk mengilustrasikannya, mulai dengan menentukan jumlah dan nampan:

counts <- c(1,2,3,4,1)
bin.lower <- c(40, 45, 50, 55, 70)
bin.upper <- c(45, 50, 55, 60, 75)

Rumus yang tepat untuk digunakan untuk penghitungan berasal dari mereplikasi lebar bin dengan jumlah yang diberikan oleh penghitungan; yaitu, data yang dikosongkan setara dengan

42.5, 47.5, 47.5, 52.5, 52.5, 57.5, 57.5, 57.5, 57.5, 72.5

xkkx2

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)

mu1195/2254.32sigma2675/1161.367.83h=5h2/12=25/122.08675/1152/127.70 untuk standar deviasi.


Estimasi Kemungkinan Maksimum

Fθθ(x0,x1] contains k values out of a set of independent, identically distributed values from Fθ, then the (additive) contribution to the log likelihood of this bin is

logi=1k(Fθ(x1)Fθ(x0))=klog(Fθ(x1)Fθ(x0))

(see MLE/Likelihood of lognormally distributed interval).

Summing over all bins gives the log likelihood Λ(θ) for the dataset. As usual, we find an estimate θ^ which minimizes Λ(θ). This requires numerical optimization and that is expedited by supplying good starting values for θ. 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 (μ^,σ^)=(54.32,7.33).

Remember, though, that for Normal distributions the maximum likelihood estimate of σ (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/(n1). Let us then (for comparison) correct the MLE of σ, finding n/(n1)σ^=11/10×7.33=7.69. This compares favorably with the result of Sheppard's correction, which was 7.70.

Verifying the Assumptions

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)

Figure

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 χ2 test: the estimated parameters indicate the expected amount of data in each bin; the χ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.

whuber
sumber