MLE / Kemungkinan interval terdistribusi secara lognormal

8

Saya memiliki serangkaian tanggapan yang dinyatakan sebagai interval seperti contoh di bawah ini.

> head(left)
[1]  860  516  430 1118  860  602
> head(right)
[1]  946  602  516 1204  946  688

di mana kiri adalah batas bawah dan kanan adalah batas atas dari respons. Saya ingin memperkirakan parameter sesuai dengan distribusi lognormal.

Untuk sementara ketika saya mencoba menghitung kemungkinan secara langsung, saya bergumul dengan fakta bahwa karena kedua batas tersebut didistribusikan di antara seperangkat paramat yang berbeda, saya mendapatkan beberapa nilai negatif seperti di bawah ini:

> Pr_high=plnorm(wta_high,meanlog_high,sdlog_high)
> Pr_low=plnorm(wta_low, meanlog_low,sdlog_low)
> Pr=Pr_high-Pr_low
> 
> head(Pr)
[1] -0.0079951419  0.0001207749  0.0008002343 -0.0009705125 -0.0079951419 -0.0022395514

Saya tidak bisa benar-benar mencari cara untuk menyelesaikannya dan memutuskan untuk menggunakan titik tengah interval sebagai gantinya yang merupakan kompromi yang baik sampai saya menemukan fungsi mledist yang mengekstrak kemungkinan log dari respons interval, ini adalah ringkasan yang saya dapatkan:

> mledist(int, distr="lnorm")
$estimate
meanlog     sdlog 
6.9092257 0.3120138 

$convergence
[1] 0

$loglik
[1] -152.1236

$hessian
         meanlog       sdlog
meanlog 570.760358    7.183723
sdlog     7.183723 1112.098031

$optim.function
[1] "optim"

$fix.arg
NULL

Warning messages:
1: In plnorm(q = c(946L, 602L, 516L, 1204L, 946L, 688L, 1376L, 1376L,  :
NaNs produced
2: In plnorm(q = c(860L, 516L, 430L, 1118L, 860L, 602L, 1290L, 1290L,  :
NaNs produced

Nilai parameter tampaknya masuk akal dan kemungkinan loglikel lebih besar daripada metode lain yang saya gunakan (distribusi titik tengah atau distribusi salah satu dari batas).

Ada pesan peringatan yang tidak saya mengerti sehingga dapatkah seseorang memberi tahu saya jika saya melakukan hal yang benar dan apa arti pesan ini?

Hargai bantuannya!

Elio Druml
sumber
Pertanyaan Anda sama dengan "Bagaimana cara saya menggunakan fungsi R tertentu dan apa arti pesan Peringatan ini?". Itu pertanyaan untuk StackOverflow daripada CrossValidated. Selanjutnya, ketika Anda merujuk ke fungsi dari suatu paket, Anda harus menyebutkan dari mana paket itu berasal . Dalam hal ini saya kira maksud Anda fungsi dari paket fitdistrplus.
Glen_b -Reinstate Monica
Selamat datang di situs ini, @ElioDruml. Saya tidak tahu apakah pertanyaan utama Anda adalah tentang cara memperkirakan parameter ini, atau apa arti pesan peringatan itu. Yang pertama akan menjadi pertanyaan yang bagus untuk CV, tetapi yang terakhir benar-benar merupakan pertanyaan untuk Stack Overflow (lihat FAQ kami ). Bisakah Anda menjelaskan apa pertanyaan utama Anda? Apakah Anda lebih suka Q Anda tinggal di sini, atau dimigrasi ke SO? (Jika yang terakhir, tandai Q Anda & kami akan memigrasikannya untuk Anda, tolong jangan posting-silang .)
gung - Reinstate Monica

Jawaban:

9

Sepertinya Anda mungkin tidak menghitung kemungkinan dengan benar.

Ketika semua yang Anda tahu tentang suatu nilaix Apakah itu

  1. Itu diperoleh secara independen dari suatu distribusi Fθ dan

  2. Itu terletak di antara a dan b>a inklusif (di mana b dan a independen dari x),

maka (menurut definisi) kemungkinannya adalah

PrFθ(axb)=Fθ(b)Fθ(a).
Kemungkinan satu set pengamatan independen karena itu adalah produk dari ekspresi seperti itu, satu per observasi. Kemungkinan log, seperti biasa, akan menjadi jumlah logaritma dari ekspresi tersebut.

Sebagai contoh, berikut ini adalah Rimplementasi di mana nilai berada di vektor , nilai di vektor , dan adalah Lognormal. (Ini bukan solusi untuk tujuan umum; khususnya, ini mengasumsikan bahwa dan untuk semua data.)aleftbrightFθb>aba

#
# Lognormal log-likelihood for interval data.
#
lambda <- function(mu, sigma, left, right) {
  sum(log(pnorm(log(right), mu, sigma) - pnorm(log(left), mu, sigma)))
}

Untuk menemukan kemungkinan log maksimum, kita memerlukan sekumpulan nilai awal yang masuk akal untuk log mean dan mencatat standar deviasi . Estimasi ini menggantikan setiap interval dengan rata-rata geometris dari titik akhir:μσ

#
# Create an initial estimate of lognormal parameters for interval data.
#
lambda.init <- function(left, right) {
  mid <- log(left * right)/2
  c(mean(mid), sd(mid))
}

Mari kita buat beberapa data acak yang didistribusikan secara lognormal dan masukkan ke dalam interval:

set.seed(17)
n <- 12                     # Number of data
z <- exp(rnorm(n, 6, .5))   # Mean = 6, SD = 0.5
left <- 100 * floor(z/100)  # Bin into multiples of 100
right <- left + 100

Pemasangan dapat dilakukan oleh pengoptimal multivariat tujuan umum. (Yang ini adalah minimizer secara default, jadi itu harus diterapkan pada kemungkinan log yang negatif.)

fit <- optim(lambda.init(left,right), 
             fn=function(theta) -lambda(theta[1], theta[2], left, right))
fit$par

6.1188785 0.3957045

Estimasi adalah , tidak jauh dari nilai yang dimaksudkan , dan estimasi adalah , tidak jauh dari nilai yang dimaksudkan : tidak buruk untuk hanya nilai. Untuk melihat seberapa baik kecocokannya, mari kita plot fungsi distribusi kumulatif empiris dan fungsi distribusi yang pas. Untuk membangun ECDF, saya hanya menyisipkan secara linear melalui setiap interval:μ6.126σ0.400.512

#
# ECDF of the data.
#
F <- function(x) (1 + mean((abs(x - left) - abs(x - right)) / (right - left)))/2

y <- sapply(x <- seq(min(left) * 0.8, max(right) / 0.8, 1), F)
plot(x, y, type="l", lwd=2, lty=2, ylab="Cumulative probability")
curve(pnorm(log(x), fit$par[1], fit$par[2]), from=min(x), to=max(x), col="Red", lwd=2, 
  add=TRUE)

Plot

Karena penyimpangan vertikal secara konsisten kecil dan bervariasi naik dan turun, sepertinya cocok.

whuber
sumber
Terima kasih banyak atas masukan Anda @whuber. Saya telah menciptakan kembali contoh Anda dan semuanya masuk akal. Namun, saya tidak dapat membuat ulang pada data saya sendiri dari n = 56 di mana kepala kiri <- c (860, 516, 430, 1118, 860, 602) dan kanan <- c (946, 602, 516 , 1204, 946, 688). Saya mendapatkan pesan peringatan ini: "1: Dalam pnorm (log (kanan), mu, sigma): NaNs diproduksi 2: In pnorm (log (kiri), mu, sigma): NaNs diproduksi" ketika pas dengan optimizer untuk mengekstrak estimasi saya. Itu membawa saya kembali ke masalah saya sebelumnya yang memiliki probabilitas negatif ketika calc. kemungkinan langkah demi langkah dan mengurangi.
Elio Druml
Ini adalah pesan peringatan yang sama yang diberikan oleh fungsi mledist dari paket fitdistrplus. Namun, seperti yang Anda lihat di atas, itu memang memberi saya output untuk estimasi mle yang terlihat relatif baik. Haruskah saya mempercayainya dan / atau apa masalahnya di sini? Terima kasih untuk umpan baliknya.
Elio Druml
Mengapa Anda tidak memposting data Anda, Elio, sehingga kami dapat mendiagnosis masalahnya? Meski begitu, saya tidak yakin ini adalah kesalahan kritis. Anda mungkin mengalami masalah yang sama yang dilaporkan oleh pengguna lain ketika secara numerik meminimalkan fungsi di Mathematica ; penjelasan yang sama mungkin berlaku dalam kasus Anda.
whuber