Apa nama metode estimasi kerapatan tempat semua pasangan kemungkinan digunakan untuk membuat distribusi campuran Normal?

12

Saya hanya memikirkan cara yang rapi (belum tentu bagus) untuk membuat perkiraan kepadatan satu dimensi dan pertanyaan saya adalah:

Apakah metode estimasi kerapatan ini memiliki nama? Jika tidak, apakah ini merupakan kasus khusus dari beberapa metode lain dalam literatur?

Berikut adalah metode: Kami memiliki vektor yang kami asumsikan diambil dari beberapa distribusi yang tidak diketahui yang ingin kami perkirakan. Cara melakukan ini adalah dengan mengambil semua pasangan nilai yang mungkin dalam X dan untuk setiap pasangan [ x i , x j ] i j cocok dengan distribusi Normal menggunakan kemungkinan maksimum. Estimasi kepadatan yang dihasilkan kemudian distribusi campuran yang terdiri dari semua Normals yang dihasilkan, di mana setiap Normal diberi bobot yang sama.X=[x1,x2,...,xn]X[xi,xj]ij

Gambar di bawah ini mengilustrasikan menggunakan metode ini pada vektor . Di sini lingkaran adalah titik data, Norma berwarna adalah distribusi kemungkinan maksimum yang diperkirakan menggunakan setiap pasangan yang mungkin dan garis hitam tebal menunjukkan perkiraan kerapatan yang dihasilkan (yaitu, distribusi campuran).[1.3,0.15,0.73,1.4]

masukkan deskripsi gambar di sini

Omong-omong, sangat mudah untuk menerapkan metode dalam R yang menarik sampel dari distribusi campuran yang dihasilkan:

# Generating some "data"
x <- rnorm(30)

# Drawing from the density estimate using the method described above.
density_estimate_sample <- replicate(9999, {
  pair <- sample(x, size = 2)
  rnorm(1, mean(pair), sd(pair))
})

# Plotting the density estimate compared with 
# the "data" and the "true" density.
hist(x ,xlim=c(-5, 5), main='The "data"')
hist(density_estimate_sample, xlim=c(-5, 5), main='Estimated density')
hist(rnorm(9999), xlim=c(-5, 5), main='The "true" density')

masukkan deskripsi gambar di sini

Rasmus Bååth
sumber
5
Cobalah metode Anda menggunakanx <- c(rnorm(30), rnorm(30, 10))
Dason
2
@ Alasan Yap, dalam hal ini metode ini tidak berfungsi sama sekali! :) Juga tidak menyatu dengan n besar.
Rasmus Bååth
4
Ini kedengarannya seperti versi estimasi kepadatan kernel yang rusak di mana bandwidth diperkirakan dengan validasi silang!
Xi'an
X=[x1,x2,,xn]n

Jawaban:

6

Ini adalah ide yang menarik, karena penduga standar deviasi tampaknya kurang sensitif terhadap outlier daripada pendekatan root-mean-square biasa. Namun, saya ragu estimator ini telah dipublikasikan. Ada tiga alasan mengapa: itu tidak efisien secara komputasi, itu bias, dan bahkan ketika bias diperbaiki, secara statistik tidak efisien (tetapi hanya sedikit). Ini dapat dilihat dengan sedikit analisis pendahuluan, jadi mari kita lakukan itu terlebih dahulu dan kemudian menarik kesimpulan.

Analisis

μσ(xi,xj)

μ^(xi,xj)=xi+xj2

dan

σ^(xi,xj)=|xixj|2.

Karena itu metode yang dijelaskan dalam pertanyaan adalah

μ^(x1,x2,,xn)=2n(n1)i>jxi+xj2=1ni=1nxi,

yang merupakan penduga rata-rata, dan

σ^(x1,x2,,xn)=2n(n1)i>j|xixj|2=1n(n1)i,j|xixj|.

E=E(|xixj|)ij

E(σ^(x1,x2,,xn))=1n(n1)i,jE(|xixj|)=E.

xixj2σ22σχ(1)2/π

E=2πσ.

2/π1.128 adalah bias dalam estimator ini.

σ^ , tapi - seperti yang akan kita lihat - ada tidak mungkin banyak minat ini, jadi saya hanya akan memperkirakan itu dengan simulasi cepat.

Kesimpulan

  1. σ^n=20,000

    Angka

  2. i,j|xixj|O(n2)O(n)n10,000R. (Pada platform lain, persyaratan RAM akan jauh lebih kecil, mungkin dengan sedikit biaya dalam waktu perhitungan.)

  3. Secara statistik tidak efisien. Untuk memberikan yang terbaik, mari kita pertimbangkan versi yang tidak bias dan bandingkan dengan versi yang tidak bias dari kuadrat terkecil atau penduga kemungkinan maksimum

    σ^OLS=(1n1i=1n(xiμ^)2)(n1)Γ((n1)/2)2Γ(n/2).

    Rn=3n=300σ^OLSσ

Kemudian

σ^


Kode

sigma <- function(x) sum(abs(outer(x, x, '-'))) / (2*choose(length(x), 2))
#
# sigma is biased.
#
y <- rnorm(1e3) # Don't exceed 2E4 or so!
mu.hat <- mean(y)
sigma.hat <- sigma(y)

hist(y, freq=FALSE,
     main="Biased (dotted red) and Unbiased (solid blue) Versions of the Estimator",
     xlab=paste("Sample size of", length(y)))
curve(dnorm(x, mu.hat, sigma.hat), col="Red", lwd=2, lty=3, add=TRUE)
curve(dnorm(x, mu.hat, sqrt(pi/4)*sigma.hat), col="Blue", lwd=2, add=TRUE)
#
# The variance of sigma is too large.
#
N <- 1e4
n <- 10
y <- matrix(rnorm(n*N), nrow=n)
sigma.hat <- apply(y, 2, sigma) * sqrt(pi/4)
sigma.ols <- apply(y, 2, sd) / (sqrt(2/(n-1)) * exp(lgamma(n/2)-lgamma((n-1)/2)))

message("Mean of unbiased estimator is ", format(mean(sigma.hat), digits=4))
message("Mean of unbiased OLS estimator is ", format(mean(sigma.ols), digits=4))
message("Variance of unbiased estimator is ", format(var(sigma.hat), digits=4))
message("Variance of unbiased OLS estimator is ", format(var(sigma.ols), digits=4))
message("Efficiency is ", format(var(sigma.ols) / var(sigma.hat), digits=4))
whuber
sumber
Literatur yang relevan kembali ke masa lalu misalnya Downton, F. 1966 Estimasi linear dengan koefisien polinomial. Biometrika 53: 129-141 doi: 10.1093 / biomet / 53.1-2.129
Nick Cox
Wow, saya mendapat lebih dari yang saya tawar! :)
Rasmus Bååth