Bagaimana cara menyesuaikan distribusi Weibull untuk memasukkan data yang mengandung nol?

14

Saya mencoba mereproduksi algoritma prediksi yang ada, yang diturunkan oleh seorang pensiunan peneliti. Langkah pertama adalah mencocokkan beberapa data yang diamati dengan distribusi Weibull, untuk mendapatkan bentuk dan skala yang akan digunakan untuk memprediksi nilai masa depan. Saya menggunakan R untuk melakukan ini. Ini contoh kode saya:

x<-c(23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,51,77,78,144,34,29,45,16,15,37,218,170,44,121)
f<-fitdistr(x, 'weibull')

Ini berfungsi dengan baik kecuali jika ada nol di array input, yang menyebabkannya gagal sepenuhnya. Hal yang sama terjadi di SAS. Seperti yang saya pahami, ini karena salah satu langkah dalam menghitung distribusi Weibull adalah mengambil log natural, yang tidak ditentukan untuk 0. Apakah ada cara yang masuk akal untuk mengatasi ini?

Sejauh ini yang terbaik yang saya temukan adalah menambahkan 1 ke semua nilai input saya, pas kurva, dan kemudian mengurangi satu dari nilai prediksi saya ("menggeser" kurva ke atas dan kemudian mundur ke bawah dengan 1). Ini cocok dengan data yang diprediksi sebelumnya dengan cukup baik, tetapi sepertinya itu cara yang salah untuk melakukannya.

sunting: Nilai-nilai dalam array input diamati, data dunia nyata (jumlah kemunculan sesuatu) untuk rentang tahun. Jadi dalam beberapa tahun jumlah kejadiannya nol. Apakah itu cara terbaik atau tidak (saya setuju mungkin tidak), pembuat algoritma asli mengklaim telah menggunakan distribusi Weibull, dan saya harus mencoba meniru proses mereka.

Ethan Shepherd
sumber
5
Weibull adalah distribusi kontinu sehingga kemungkinan mendapatkan tepat nol memiliki probabilitas nol. Jika Anda mendapatkan banyak nol di data Anda, itu adalah petunjuk langsung bahwa Weibull tidak pantas. Bagaimanapun, data Anda terlihat seperti data jumlah (atau paling tidak, diskrit) sehingga Weibull mungkin bukan pilihan terbaik.
kardinal
Menambahkan beberapa konteks dari mana data berasal akan membantu siapa pun yang mencoba menjawab dengan sangat.
kardinal

Jawaban:

8

(Seperti yang telah ditunjukkan orang lain, distribusi Weibull tidak mungkin menjadi perkiraan yang tepat ketika data hanya bilangan bulat. Berikut ini dimaksudkan hanya untuk membantu Anda menentukan apa yang dilakukan peneliti sebelumnya, benar atau salah.)

Ada beberapa metode alternatif yang tidak terpengaruh oleh nol dalam data, seperti menggunakan berbagai metode penaksir momen. Ini biasanya memerlukan solusi numerik persamaan yang melibatkan fungsi gamma, karena momen-momen distribusi Weibull diberikan dalam hal fungsi ini. Saya tidak terbiasa dengan R, tapi inilah program Sage yang menggambarkan salah satu metode yang lebih sederhana - mungkin bisa disesuaikan dengan R? (Anda dapat membaca tentang ini dan metode lain semacam itu di, misalnya, "Distribusi Weibull: buku pegangan" oleh Horst Rinne, hal. 455ff - namun, ada kesalahan ketik pada eq.12.4b-nya, sebagai '-1' redundan).

"""
Blischke-Scheuer method-of-moments estimation of (a,b)
for the Weibull distribution F(t) = 1 - exp(-(t/a)^b)
""" 

x = [23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,16,15,37,218,170,44,121]
xbar = mean(x)
varx = variance(x)
var("b"); f(b) = gamma(1+2/b)/gamma(1+1/b)^2 - 1 - varx/xbar^2
bhat = find_root(f, 0.01, 100)
ahat = xbar/gamma(1+1/bhat)
print "Estimates: (ahat, bhat) = ", (ahat, bhat)

Ini menghasilkan output

Estimates: (ahat, bhat) =  (81.316784310814455, 1.3811394719075942)


0

x = [23,0,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,0,0,37,218,170,44,121]

maka prosedur yang sama menghasilkan output

Estimates: (ahat, bhat) =  (78.479354097488923, 1.2938352346035282)


EDIT: Saya baru saja menginstal R untuk mencobanya. Dengan risiko membuat jawaban ini terlalu lama, bagi siapa pun yang tertarik inilah kode-R saya untuk metode Blischke-Scheuer:

fit_weibull <- function(x)
{
    xbar <- mean(x)
    varx <- var(x)
    f <- function(b){return(gamma(1+2/b)/gamma(1+1/b)^2 - 1 - varx/xbar^2)}
    bhat <- uniroot(f,c(0.02,50))$root
    ahat <- xbar/gamma(1+1/bhat)
    return(c(ahat,bhat))
}

Ini mereproduksi (hingga lima digit signifikan) dua contoh Sage di atas:

x <- c(23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,
     51,77,78,144,34,29,45,16,15,37,218,170,44,121)
fit_weibull(x)
[1] 81.316840  1.381145

x <- c(23,0,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,0,0,37,218,170,44,121)
fit_weibull(x)
[1] 78.479180  1.293821
res
sumber
4

θfitdistrθθfitdistr

foo <- function(theta, x)
{
  if (theta <= -min(x)) return(Inf);
  f <- fitdistr(x+theta, 'weibull')
  -2*f$loglik
}

Kemudian minimalkan fungsi ini menggunakan optimasi satu dimensi:

bar <- optimize(foo, lower=-min(x)+0.001, upper=-min(x)+10, x=x)

di mana saya baru saja membuat "+10" berdasarkan apa-apa.

Untuk data dengan tiga nilai terkecil yang diganti dengan nol, kita mendapatkan:

> bar
$minimum
[1] 2.878442

$objective
[1] 306.2792

> fitdistr(x+bar$minimum, 'weibull')
     shape        scale   
   1.2836432   81.1678283 
 ( 0.1918654) (12.3101211)
> 

bar$minimum adalah MLE dari θfitdistrθ

Jbowman
sumber
2

Itu harus gagal, Anda harus bersyukur bahwa itu gagal.

Pengamatan Anda menunjukkan bahwa kegagalan terjadi pada saat Anda mulai mengamatinya. Jika ini adalah proses nyata, yang berasal dari data nyata (dan bukan data yang disimulasikan), Anda perlu menjelaskan alasan mengapa Anda mendapatkan angka nol. Saya telah melihat studi bertahan hidup di mana 0 kali muncul sebagai konsekuensi dari salah satu dari beberapa hal:

  1. Data sebenarnya terpotong: objek beresiko dan gagal sebelum penelitian dimulai dan Anda ingin berpura-pura telah mengamati semuanya selama ini.
  2. Instrumen dikalibrasi dengan buruk: Anda tidak memiliki presisi pengukuran yang cukup untuk penelitian ini dan kegagalan yang terjadi di dekat waktu mulai dikodekan sebagai tepat nol.
  3. Sesuatu yang dikodekan sebagai nol bukanlah nol. Mereka adalah orang atau objek yang dikeluarkan dari analisis dengan satu atau lain cara. Angka nol hanya muncul di data sebagai konsekuensi dari penggabungan, pengurutan, atau pengodean ulang nilai yang hilang.

Jadi untuk kasus 1: Anda perlu menggunakan metode sensor yang tepat, bahkan jika itu berarti menarik catatan secara retrospektif. Kasus 2 berarti Anda dapat menggunakan algoritma EM karena Anda memiliki masalah presisi. Metode Bayesian bekerja dengan cara yang sama di sini. Kasus 3 berarti Anda hanya perlu mengecualikan nilai-nilai yang seharusnya hilang.

AdamO
sumber
OP menjelaskan bahwa peneliti sebelumnya memilih agar sesuai dengan distribusi Weibull, meskipun datanya adalah jumlah dunia nyata - jumlah integer non-negatif dari jumlah kemunculan sesuatu. Tidak jelas bagaimana ketiga kasus Anda berhubungan dengan situasi seperti itu.
res
Oh, catatan bagus! Menyesuaikan distribusi Weibull sangat salah. Ini memiliki dukungan terus menerus dan tidak pernah digunakan untuk memodelkan jumlah tetapi waktu bertahan hidup. Distribusi binomial negatif akan menjadi semacam distribusi dua parameter yang ekuivalen untuk jumlah pemodelan, yang tentu saja tergantung pada sifat proses pembuatan data (yang mana kami memiliki 0 informasi, seperti masalahnya dinyatakan). Terima kasih telah menunjukkan itu padaku.
AdamO
1

Saya setuju dengan jawaban kardinal di atas. Namun, juga cukup umum untuk menambahkan konstanta untuk menghindari nol. Nilai lain yang umum digunakan adalah 0,5, tetapi konstanta positif apa pun mungkin telah digunakan. Anda dapat mencoba rentang nilai untuk melihat apakah Anda dapat mengidentifikasi nilai tepat yang digunakan oleh peneliti sebelumnya. Maka Anda bisa yakin bahwa Anda dapat mereproduksi hasilnya, sebelum mencari distribusi yang lebih baik.

John Bauer
sumber
0

[Dengan asumsi Weibull tepat] Buku Johnson Kotz dan Balakrishnan memiliki banyak cara untuk memperkirakan parameter Weibull. Beberapa di antaranya tidak bergantung pada data yang tidak termasuk nol (misalnya menggunakan mean dan standar deviasi, atau menggunakan persentil tertentu).

Johnson, NL, Kotz, S., dan Balakrishnan, N. (1994). Distribusi Univariat Berkelanjutan. New York: Wiley, kira-kira di halaman 632.

zbicyclist
sumber