Bias penaksir momen dari distribusi lognormal

25

Saya melakukan beberapa percobaan numerik yang terdiri dari pengambilan sampel distribusi lognormal XLN(μ,σ) , dan mencoba memperkirakan momen dengan dua metode:E[Xn]

  1. Melihat rata-rata sampelXn
  2. Memperkirakan dan dengan menggunakan mean sampel untuk , dan kemudian menggunakan fakta bahwa untuk distribusi lognormal, kita memiliki .σ 2μσ2E [ X n ] = exp ( n μ + ( n σ ) 2 / 2 )log(X),log2(X)E[Xn]=exp(nμ+(nσ)2/2)

Pertanyaannya adalah :

Saya menemukan secara eksperimental bahwa metode kedua berkinerja jauh lebih baik daripada yang pertama, ketika saya menjaga jumlah sampel tetap, dan meningkatkan oleh beberapa faktor T. Apakah ada penjelasan sederhana untuk fakta ini?μ,σ2

Saya melampirkan gambar di mana sumbu x adalah T, sedangkan sumbu y adalah nilai membandingkan nilai sebenarnya dari (garis oranye), ke nilai yang diestimasi. metode 1 - titik biru, metode 2 - titik hijau. sumbu y dalam skala logE [ X 2 ] = exp ( 2 μ + 2 σ 2 )E[X2]E[X2]=exp(2μ+2σ2)

Nilai true dan estimasi untuk $ \ mathbb {E} [X ^ 2] $.  Titik biru adalah rata-rata sampel untuk $ \ mathbb {E} [X ^ 2] $ (metode 1), sedangkan titik hijau adalah nilai yang diestimasi menggunakan metode 2. Garis oranye dihitung dari $ \ mu $, $ \ yang dikenal sigma $ dengan persamaan yang sama seperti pada metode 2. sumbu y adalah dalam skala log

EDIT:

Di bawah ini adalah kode Mathematica minimal untuk menghasilkan hasil untuk satu T, dengan output:

   ClearAll[n,numIterations,sigma,mu,totalTime,data,rmomentFromMuSigma,rmomentSample,rmomentSample]
(* Define variables *)
n=2; numIterations = 10^4; sigma = 0.5; mu=0.1; totalTime = 200;
(* Create log normal data*)
data=RandomVariate[LogNormalDistribution[mu*totalTime,sigma*Sqrt[totalTime]],numIterations];

(* the moment by theory:*)
rmomentTheory = Exp[(n*mu+(n*sigma)^2/2)*totalTime];

(*Calculate directly: *)
rmomentSample = Mean[data^n];

(*Calculate through estimated mu and sigma *)
muNumerical = Mean[Log[data]]; (*numerical \[Mu] (gaussian mean) *)
sigmaSqrNumerical = Mean[Log[data]^2]-(muNumerical)^2; (* numerical gaussian variance *)
rmomentFromMuSigma = Exp[ muNumerical*n + (n ^2sigmaSqrNumerical)/2];

(*output*)
Log@{rmomentTheory, rmomentSample,rmomentFromMuSigma}

Keluaran:

(*Log of {analytic, sample mean of r^2, using mu and sigma} *)
{140., 91.8953, 137.519}

di atas, hasil kedua adalah rerata sampel r2 , yang di bawah dua hasil lainnya

pengguna29918
sumber
2
Estimator yang tidak bias tidak menyiratkan bahwa titik-titik biru harus mendekati nilai yang diharapkan (kurva oranye). Estimator dapat tidak memihak jika memiliki probabilitas tinggi terlalu rendah dan probabilitas kecil (mungkin sangat kecil) terlalu tinggi. Itulah yang terjadi ketika T meningkat dan varians menjadi sangat besar (lihat jawaban saya).
Matthew Gunn
Untuk cara mendapatkan estimator yang tidak bias, silakan lihat stats.stackexchange.com/questions/105717 . UMVUE dari mean dan varians diberikan dalam jawaban dan komentar di dalamnya.
whuber

Jawaban:

22

Sejak itu ada sesuatu yang membingungkan

  1. metode pertama memberikan penduga yang tidak bias dari , yaitu 1E[X2] memilikiE[X2]sebagai artinya. Karenanya titik-titik biru harus berada di sekitar nilai yang diharapkan (kurva oranye);
    1Ni=1NXi2
    E[X2]
  2. metode kedua menyediakan estimator bias dari , yaitu E [ exp ( n μ + n 2 σ 2 / 2 ) ] > exp ( n μ + ( n σ ) 2 / 2 ) ketika μ dan σ ² adalah estimator berisi dari μ dan σ ²E[X2]
    E[exp(nμ^+n2σ^2/2)]>exp(nμ+(nσ)2/2)
    μ^σ^²μσ² masing-masing, dan dengan demikian aneh bahwa titik-titik hijau selaras dengan kurva oranye.

tetapi mereka disebabkan oleh masalah dan bukan karena perhitungan numerik: Saya mengulangi percobaan dalam R dan mendapatkan gambar berikut dengan kode warna yang sama dan urutan yang sama dari dan σ T , yang mewakili masing-masing penduga yang dibagi dengan harapan sebenarnya:μTσT

Dua momen empiris kedua, berdasarkan pada 10⁶ simulasi log-normal

Berikut adalah kode R yang sesuai:

moy1=moy2=rep(0,200)
mus=0.14*(1:200)
sigs=sqrt(0.13*(1:200))
tru=exp(2*mus+2*sigs^2)
for (t in 1:200){
x=rnorm(1e5)
moy1[t]=mean(exp(2*sigs[t]*x+2*mus[t]))
moy2[t]=exp(2*mean(sigs[t]*x+mus[t])+2*var(sigs[t]*x+mus[t]))}

plot(moy1/tru,col="blue",ylab="relative mean",xlab="T",cex=.4,pch=19)
abline(h=1,col="orange")
lines((moy2/tru),col="green",cex=.4,pch=19)

Oleh karena itu memang ada keruntuhan momen empiris kedua sebagai dan σ meningkat yang saya akan atribut untuk peningkatan besar dalam varian dari momen empiris kedua tersebut sebagai μ dan σ meningkat.μσμσ

E[X2]X2X2e2μX2exp{2μ+2σϵ}ϵN(0,1)σσϵσ2XL.N(μ,σ)

P(X2>E[X2])=P(log{X2}>2μ+2σ2)=P(μ+σϵ>μ+σ2)=P(ϵ>σ)=1-Φ(σ)
Xi'an
sumber
1
Saya juga bingung. Saya menambahkan kode minimal dengan hasil (Mathematica)
user29918
Baik. Terima kasih! Menempatkan beberapa angka, saya melihat sekarang bahwa ukuran sampel saya yang sedikit benar-benar tidak cocok untuk tugas itu!
user29918
2
σ
2
P(X2>E[X2])=1-Φ(σ)σσ
2
σ
13

Saya pikir saya akan memunculkan beberapa buah ara yang menunjukkan bahwa plot user29918 dan Xi'an konsisten. Gambar 1 memplot apa yang user29918 lakukan, dan Gambar 2 (berdasarkan data yang sama), melakukan apa yang Xi'an lakukan untuk plotnya. Hasil yang sama, presentasi berbeda.

1nsayaxsaya2

Komentar Lebih Lanjut:

  1. Penaksir yang tidak bias tidak berarti penaksir diharapkan akan menutup! Titik-titik biru tidak harus mendekati harapan. Misalnya. sebuah pengamatan tunggal yang dipilih secara acak memberikan perkiraan rata-rata populasi, tetapi estimator itu tidak akan diharapkan menjadi dekat.
  2. Masalahnya muncul karena varians menjadi benar-benar astronomi. Sebagai varians pergi batshit, estimasi untuk metode pertama sedang didorong menjadi hanya beberapa pengamatan. Anda juga mulai memiliki probabilitas kecil, sangat kecil dari Gila, GILA-GILA, GILA-GILA besar ...
  3. P(X2>E[X2])=1-Φ(σ)σX2>E[X2]masukkan deskripsi gambar di sini

masukkan deskripsi gambar di sini

Matthew Gunn
sumber