Mengestimasi Distribusi Frekuensi
Estimasi Parameter dan Distribusi Frekuensi MLE
| Kontak | : \(\downarrow\) |
| garryjuliuspermana@gmail.com | |
| https://www.instagram.com/boring.garr/ | |
| RPubs | https://rpubs.com/garr/ |
2.4 Mengestimasi Distibusi Frekuensi
Pada bagian ini, Anda akan mempelajari bagaimana caranya:
Menentukan kemungkinan atau probability suatu sampel pengamatan dari distribusi diskrit
Menentukan penaksir kemungkinan maksimum untuk sampel acak pengamatan dari distribusi diskrit
Menghitung penaksir kemungkinan maksimum untuk distribusi binomial, distribusi Poisson, dan distribusi binomial negative
2.4.1 Estimasi Parameter
Pada Bagian 2.2 kita telah mempelajari tiga distribusi yang penting dalam memodelkan berbagai jenis data perhitungan yang timbul dari asuransi. Sekarang mari kita anggap bahwa kita memiliki sekumpulan data perhitungan yang ingin kita sesuaikan dengan sebuah distribusi, dan kita telah menentukan bahwa salah satu dari distribusi ini \((a,b,0)\)
distribusi ini lebih tepat daripada yang lain. Karena masing-masing membentuk kelas distribusi jika kita mengizinkan parameter atau nilainya untuk mengambil nilai apa pun yang diizinkan, masih ada tugas untuk menentukan nilai terbaik dari parameter untuk data yang ada. Ini adalah masalah estimasi titik statistik, dan dalam masalah inferensi parametrik, paradigma inferensi statistik kemungkinan maksimum biasanya menghasilkan perkiraan yang efisien. Pada bagian ini kita akan menjelaskan paradigma ini dan menurunkan perkiraan kemungkinan maksimum.
Misalkan kita mengamati bahwa kita mengamati independen dan terdistribusi secara identik, iid, variabel acak \(X_1,X_2,…,X_n\) dari distribusi dengan Prabability Mass Function/pmf \(pθ\), dimana \(θ\) adalah vektor parameter dan nilai yang tidak diketahui dalam ruang parameter \(Θ⊆\mathbb{R}^d\). Sebagai contoh, dalam kasus distribusi Poisson, ada satu parameter sehingga \(d=1\) dan
\[\begin{align*} p_\theta(x)=e^{-\theta}\frac{\theta^x}{x!}, \quad x=0,1,\ldots, \end{align*}\]
dengan \(θ>0\). Dalam kasus distribusi binomial, kita memiliki
\[\begin{align*} p_\theta(x)= {m \choose x} q^x(1-q)^{m-x}, \quad x=0,1,\ldots,m. \end{align*}\]
Untuk beberapa aplikasi, kita dapat melihat m sebagai sebuah parameter dan dengan demikian ambil \(d = 2\) sehingga \(θ = (m,q) ∈{0,1,2,...}× [0,1]\).
Misalkan pengamatannya adalah \(x_1,...,x_n\) nilai pengamatan dari sampel acak \(X_1,X_2,...,X_n\) yang disajikan sebelumnya. Dalam kasus ini, probabilitas pengamatan sampel ini dari \(pθ\) sama dengan
\[\begin{align*} \prod_{i=1}^n p_\theta(x_i). \end{align*}\]
Hal di atas, dilambangkan dengan \(L(θ)\) yang dipandang sebagai fungsi dari \(θ\) disebut dengan kemungkinan. Perhatikan bahwa kita telah menekan ketergantungannya pada data, untuk menekankan bahwa kita melihatnya sebagai fungsi dari vektor parameter. Sebagai contoh, dalam kasus distribusi Poisson, kita memiliki
\[\begin{align*} L(\lambda)=e^{-n\lambda} \lambda^{\sum_{i=1}^n x_i} \left(\prod_{i=1}^n x_i!\right)^{-1}. \end{align*}\]
Dalam kasus distribusi binomial, kita memiliki
\[\begin{align*} L(m,q)=\left(\prod_{i=1}^n {m \choose x_i}\right) q^{\sum_{i=1}^n x_i} (1-q)^{nm-\sum_{i=1}^n x_i} . \end{align*}\]
Penaksir kemungkinan maksimum (mle) untuk \(θ\) adalah pemaksimum dari kemungkinan; dalam artian, mle memilih himpunan nilai parameter yang paling baik menjelaskan pengamatan yang diamati.
Kasus Khusus: Tiga Hasil Bernoulli. Sebagai ilustrasi, pertimbangkan sebuah sampel dengan ukuran \(n=3\) dari distribusi Bernoulli (binomial dengan \(m = 1\)) dengan nilai \(0,1,0\). Peluang dalam kasus ini dengan mudah diperiksa sama dengan
\[\begin{align*} L(q)=q(1-q)^2, \end{align*}\]
dan plot kemungkinan diberikan pada Gambar 2.1. Seperti yang ditunjukkan pada plot, nilai maksimum kemungkinan sama dengan \(4/27\) dan dicapai pada \(q = 1/3\), dan karenanya estimasi kemungkinan maksimum untuk \(q\) adalah \(1/3\) untuk sampel yang diberikan. Dalam hal ini, kita dapat menggunakan aljabar untuk menunjukkan bahwa
\[\begin{align*} q(1-q)^2=\left(q-\frac{1}{3}\right)^2\left(q-\frac{4}{3}\right)+\frac{4}{27}, \end{align*}\]
dan menyimpulkan bahwa nilai maksimumnya sama dengan \(4/27\), dan dicapai pada \(q = 1/3\) (menggunakan fakta bahwa suku pertama tidak positif dalam interval \([0,1]\)).
Tetapi seperti yang terlihat, cara menurunkan mle menggunakan aljabar ini tidak berlaku umum. Secara umum, seseorang menggunakan kalkulus untuk menurunkan mle - perhatikan bahwa untuk beberapa kemungkinan, seseorang mungkin harus menggunakan metode optimasi lainnya, terutama ketika kemungkinan tersebut memiliki banyak ekstrema lokal. Sudah menjadi kebiasaan untuk memaksimalkan secara ekuivalen logaritma dari kemungkinan \(L(⋅)\), dilambangkan dengan \(l(⋅\)), dan melihat himpunan nol dari turunan pertamanya$ l′(⋅)$. Dalam kasus kemungkinan di atas, \(l(q) = log(q) + 2log(1-q)\), dan
\[\begin{align*} l'(q)=\frac{\rm d}{{\rm d}q}l(q)=\frac{1}{q}-\frac{2}{1-q}. \end{align*}\]
Angka nol unik dari \(l′(⋅)\) sama dengan \(1/3\) , dan karena \(l′′(⋅)\) adalah negatif, kita memiliki \(1/3\) adalah pemaksimum unik dari kemungkinan dan karenanya merupakan estimasi kemungkinan maksimum.
library(png)
img <- readPNG("figur.png")
plot(1:3, 1:3, type = "n", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
rasterImage(img, 1, 1, 3, 3)
2.4.2 Distribusi Frekuensi MLE
Berikut ini, kami menurunkan penaksir kemungkinan maksimum, mle(Maximum likelihood estimator), untuk tiga anggota kelas \((a,b,0)\). Kita mulai dengan meringkas pembahasan di atas. Dalam situasi mengamati iid, independen dan berdistribusi identik, variabel acak \(X_1, X_2,...,X_n\) dari sebuah distribusi dengan pmf \(p_θ\), di mana \(θ\) memiliki nilai yang tidak diketahui dalam \(Θ⊆R^d\) , kemungkinan \(L(⋅)\), sebuah fungsi pada \(Θ\) didefinisikan sebagai
\[\begin{align*} L(\theta)=\prod_{i=1}^n p_\theta(x_i), \end{align*}\]
di mana \(x_1,...,x_n\) adalah nilai yang diamati. MLE dari \(θ\), dinotasikan sebagai \(\hat{θ}_{MLE}\), adalah sebuah fungsi yang memetakan observasi ke sebuah elemen dari himpunan pemaksimum \(L(⋅)\), yaitu
\[\begin{align*} \{\theta \vert L(\theta)=\max_{\eta\in\Theta}L(\eta)\}. \end{align*}\]
Perhatikan bahwa himpunan di atas merupakan fungsi dari pengamatan, meskipun ketergantungan ini tidak dibuat secara eksplisit. Dalam kasus tiga distribusi yang kita pelajari, dan secara umum, himpunan di atas adalah himpunan tunggal dengan probabilitas yang cenderung mendekati satu (dengan meningkatnya ukuran sampel). Dengan kata lain, untuk banyak distribusi yang umum digunakan dan ketika ukuran sampel besar, estimasi kemungkinan didefinisikan secara unik dengan probabilitas yang tinggi.
Berikut ini, kita asumsikan bahwa kita telah mengamati n variabel acak ke-i \(X_1,X_2,...,X_n\) dari distribusi yang dipertimbangkan, meskipun nilai parametriknya tidak diketahui. Selain itu, \(x_1,x_2,...,x_n\) akan menunjukkan nilai yang diamati. Kami mencatat bahwa dalam kasus data cacahan, dan data dari distribusi diskrit secara umum, kemungkinannya dapat direpresentasikan sebagai
\[\begin{align*} L(\theta)=\prod_{k\geq 0} \left(p_\theta(k)\right)^{m_k}, \end{align*}\]
dimana \(m_k\) adalah jumlah observasi yang sama dengan \(k\) . Secara matematis, kita memiliki
\[\begin{align*} m_k= \left\vert \{i\vert x_i=k, 1\leq i \leq n\} \right\vert=\sum_{i= 1}^n I(x_i=k), \quad k\geq 0. \end{align*}\]
Perhatikan bahwa transformasi ini mempertahankan semua data, menyusunnya dengan cara yang efisien. Untuk \(n\) yang besar hal ini menyebabkan kompresi data dalam arti kecukupan. Di bawah ini, kami menyajikan ekspresi untuk mle dalam hal \((m_k)_{k≥1}\) juga.
Kasus Khusus: Distribusi Poisson. Dalam kasus ini, seperti yang disebutkan di atas, kemungkinan diberikan oleh
\[\begin{align*} L(\lambda)=\left(\prod_{i=1}^n x_i!\right)^{-1}e^{-n\lambda}\lambda^{\sum_{i=1}^n x_i} . \end{align*}\]
Dengan menggunakan logaritma, log-kemungkinan adalah
\[\begin{align*} l(\lambda)= -\sum_{i=1}^n \log(x_i!) -n\lambda +\log(\lambda) \cdot \sum_{i=1}^n x_i . \end{align*}\]
Mengambil turunannya, kami memiliki
\[\begin{align*} l'(\lambda)= -n +\frac{1}{\lambda}\sum_{i=1}^n x_i. \end{align*}\]
Dalam mengevaluasi \(l′′(λ)\), ketika \(∑^n_{i=1}xi>0\), \(l′′<0\). Akibatnya, maksimum dicapai pada rata-rata sampel, \(\bar{x}\) yang disajikan di bawah ini. Ketika \(∑^n_{i=1}xi=0\) maka kemungkinan adalah fungsi menurun dan karenanya maksimum dicapai pada nilai parameter yang paling kecil; hal ini mengakibatkan estimasi kemungkinan maksimum menjadi nol. Oleh karena itu, kita memiliki
\[\begin{align*} \overline{x} = \hat{\lambda}_{\rm MLE} = \frac{1}{n}\sum_{i=1}^n x_i. \end{align*}\]
Perhatikan bahwa rata-rata sampel juga dapat dihitung sebagai
\[\begin{align*} \overline{x} = \frac{1}{n} \sum_{k\geq 1} k \cdot m_k ~. \end{align*}\]
Patut dicatat bahwa dalam kasus Poisson, distribusi yang tepat dari \(λ_{MLE}\) tersedia dalam bentuk tertutup - ini adalah Poisson berskala - ketika distribusi yang mendasarinya adalah Poisson. Hal ini dikarenakan jumlah variabel acak Poisson independen adalah Poisson juga. Tentu saja, untuk ukuran sampel yang besar, seseorang dapat menggunakan Teorema Batas Tengah/Central Limit Theorem (CLT) biasa untuk mendapatkan perkiraan normal. Perhatikan bahwa perkiraan yang terakhir berlaku bahkan jika distribusi yang mendasari adalah distribusi apa pun dengan momen kedua yang terbatas.
Kasus Khusus: Distribusi Binomial. Tidak seperti kasus distribusi Poisson, ruang parameter dalam kasus binomial adalah \(2\) dimensi. Oleh karena itu, masalah optimasi sedikit lebih menantang. Kita mulai dengan mengamati bahwa kemungkinan diberikan oleh
\[\begin{align*} L(m,q)= \left(\prod_{i=1}^n {m \choose x_i}\right) q^{\sum_{i=1}^n x_i} (1-q)^{nm-\sum_{i=1}^n x_i} . \end{align*}\]
Dengan menggunakan logaritma, log-kemungkinan adalah
\[\begin{align*} \begin{array}{ll} l(m,q) &= \sum_{i=1}^n \log\left({m \choose x_i}\right) + \left({\sum_{i=1}^n x_i}\right)\log(q) \\ & \ \ \ + \left({nm-\sum_{i=1}^n x_i}\right)\log(1-q) \\ &= \sum_{i=1}^n \log\left({m \choose x_i}\right) + n \overline{x}\log(q) + n\left({m- \overline{x}}\right)\log(1-q) , \end{array} \end{align*}\]
di mana \(\bar{x} = n^{-1}∑^n_{i=1}xi\). Perhatikan bahwa karena m hanya mengambil nilai bilangan bulat non-negatif, kita tidak dapat menggunakan kalkulus multivariat untuk menemukan nilai optimal. Namun demikian, kita dapat menggunakan kalkulus variabel tunggal untuk menunjukkan bahwa
\[\begin{equation} \hat{q}_{MLE}\times \hat{m}_{MLE} = \overline{x}. \tag{2.2} \end{equation}\]
verifikasi persamaan 2.2
Terhadap hal ini, kami mencatat bahwa untuk nilai \(m\) yang tetap
\[\begin{align*} \frac{\delta}{\delta q} l(m,q) = \frac{n \overline{x}}{q}- \frac{n\left({m- \overline{x}}\right)}{1-q}, \end{align*}\]
dan itu
\[\begin{align*} \frac{\delta^2}{\delta q^2} l(m,q)= -\frac{n \overline{x}}{q^2}+ \frac{n\left({m- \overline{x}}\right)}{(1-q)^2} \le 0. \end{align*}\]
Hal di atas mengimplikasikan bahwa untuk setiap nilai \(m\) yang tetap yang tetap, nilai maksimum dari \(q\) memenuhi
\[\begin{align*} mq=\overline{x}, \end{align*}\]
dan karenanya kita membuat persamaan (2.2)
Dengan persamaan (2.2), hal di atas mereduksi tugas menjadi pencarian \(\hat{m}_{MLE}\) yang merupakan pemaksimum dari
\[\begin{equation} L\left(m,\frac{\overline{x}}{m} \right). \tag{2.3} \end{equation}\]
Perhatikan bahwa kemungkinan akan menjadi nol untuk nilai \(m\) yang lebih kecil dari \(max_{1≤i≤n}xi\), dan karenanya \(\hat{m}_{MLE}≥max_{1≤i≤n}xi\).
Catatan Teknis tentang Perkiraan Poisson untuk Binomial
Dalam menentukan algoritma untuk menghitung \(\hat{m}_{MLE}\), pertama-tama kami tunjukkan bahwa untuk beberapa set data, \(\hat{m}_{MLE}\) dapat sama dengan \(∞\), yang mengindikasikan bahwa distribusi Poisson akan memberikan kecocokan yang lebih baik dibandingkan distribusi binomial mana pun. Hal ini terjadi karena distribusi binomial dengan parameter \((m, \bar{x}/m)\) mendekati distribusi Poisson dengan parameter \(\bar{x}\) dengan \(m\) mendekati tak terhingga. Fakta bahwa beberapa set data lebih memilih distribusi Poisson seharusnya tidak mengherankan karena dalam pengertian di atas, himpunan distribusi Poisson berada pada batas himpunan distribusi binomial. Menariknya, dalam Olkin, Petkau, dan Zidek (1981) mereka menunjukkan bahwa jika rata-rata sampel kurang dari atau sama dengan varians sampel maka \(\hat{m}_{MLE} = ∞\) jika tidak, maka terdapat sebuah \(m\) berhingga yang memaksimumkan persamaan (2.3).
Pada Gambar 2.2 di bawah ini, kami menampilkan plot \(L(m, \bar{x}/m)\) untuk tiga sampel yang berbeda dengan ukuran \(5\); mereka hanya berbeda dalam nilai maksimum sampel. Sampel pertama dari \((2,2,2,4,5)\) memiliki rasio rata-rata sampel terhadap varians sampel yang lebih besar dari \(1\) \((1.875)\), sampel kedua dari \((2,2,2,4,6)\) memiliki rasio sebesar \(1.25\) yang lebih mendekati \(1\), dan sampel ketiga dari \((2,2,2,4,6)\) memiliki rasio sebesar \(1.25\) yang lebih mendekati \(1\), dan sampel keempat dari \((2,2,2,4,6)\) memiliki rasio sebesar \(1.25\) yang lebih mendekati \(1\). dan sampel ketiga \((2,2,2,4,7)\) memiliki rasio lebih kecil dari \(1\) \((0.885)\). Untuk ketiga sampel ini, seperti yang ditunjukkan pada Gambar 2.2, \(\bar{m}_{MLE}\) sama dengan \(7\), \(18\) dan \(∞\) masing-masing. Perhatikan bahwa nilai pembatas dari \(L(m, \bar{x}/m)\) sebagai \(m\) mendekati tak terhingga sama dengan
\[\begin{equation} \left(\prod_{i=1}^n x_i! \right)^{-1} \exp\left(-n \overline{x}~\right) \left(~\overline{x}~\right)^{n\overline{x}}. \tag{2.4} \end{equation}\]
Juga, perhatikan bahwa Gambar 2.2 menunjukkan bahwa mle dari \(m\) tidak kuat, yaitu perubahan pada sebagian kecil set data dapat menyebabkan perubahan besar pada penaksir.
Diskusi di atas menyarankan algoritma sederhana berikut ini:
- Langkah 1. Jika rata-rata sampel kurang dari atau sama dengan varians sampel, maka tetapkan \(\hat{m}_{MLE}=∞\). Distribusi yang disarankan oleh MLE adalah distribusi Poisson dengan \(\hat{λ}=\bar{x}\).
- Langkah 2. Jika rata-rata sampel lebih besar dari varians sampel, maka hitunglah \(L(m,\bar{x}/m)\) untuk \(m\) yang lebih besar atau sama dengan maksimum sampel sampai \(L(m,\bar{x}/m)\) mendekati nilai kemungkinan Poisson yang diberikan dalam (2.4). Nilai \(m\) yang sesuai dengan nilai maksimum \(L(m,\bar{x}/m)\) di antara yang dihitung sama dengan \(\hat{m}_{MLE}\).
Kami mencatat bahwa jika distribusi yang mendasari adalah distribusi binomial dengan parameter \((m,q)\) (dengan \(q>0\)) maka \(\hat{m}_{MLE}\) sama dengan \(m\) untuk ukuran sampel yang besar. Juga, \(\hat{q}_{MLE}\) akan memiliki distribusi normal asimtotik dan konvergen dengan probabilitas satu untuk \(q\) .
library(png)
img <- readPNG("figur2.png")
plot(2:3, 2:3, type = "n", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
rasterImage(img, 2, 2, 3, 3)
Kasus Khusus: Distribusi Binomial Negatif. Kasus distribusi binomial negatif mirip dengan distribusi binomial dalam arti kita memiliki dua parameter dan mles tidak tersedia dalam bentuk tertutup. Perbedaan di antara keduanya adalah bahwa tidak seperti parameter binomial \(m\) yang mengambil nilai bilangan bulat positif, parameter \(r\) dari binomial negatif dapat mengambil nilai real positif. Hal ini membuat masalah optimasi menjadi sedikit lebih kompleks. Kita mulai dengan mengamati bahwa kemungkinan dapat dinyatakan dalam bentuk berikut:
\[\begin{align*} L(r,\beta)=\left(\prod_{i=1}^n {r+x_i-1 \choose x_i}{r+x_i-1 \choose x_i}\right) (1+\beta)^{-n(r+\overline{x})} \beta^{n\overline{x}}. \end{align*}\]
Hal di atas menyiratkan bahwa log-kemungkinan diberikan oleh
\[\begin{align*} l(r,\beta)=\sum_{i=1}^n \log{r+x_i-1 \choose x_i} -n(r+\overline{x}) \log(1+\beta) +n\overline{x}\log\beta, \end{align*}\]
dan karenanya
\[\begin{align*} \frac{\delta}{\delta\beta} l(r,\beta) = -\frac{n(r+\overline{x})}{1+\beta} + \frac{n\overline{x}}{\beta}. \end{align*}\]
Dengan menyamakan nilai di atas dengan nol, kita mendapatkan
\[\begin{align*} \frac{\delta}{\delta\beta} l(r,\beta) = -\frac{n(r+\overline{x})}{1+\beta} + \frac{n\overline{x}}{\beta}. \end{align*}\]
Hal di atas mereduksi masalah optimasi dua dimensi menjadi masalah satu dimensi - kita perlu memaksimalkan
\[\begin{align*} \hat{r}_{MLE}\times \hat{\beta}_{MLE} = \overline{x}. \end{align*}\]
terhadap \(r\), dengan \(r\) yang memaksimumkan adalah mle dan \(\hat{β}_{MLE}=\bar{x}/\hat{r}_{MLE}\). Dalam Levin, Reeds, dan kawan-kawan (1977) ditunjukkan bahwa jika varians sampel lebih besar daripada rata-rata sampel maka terdapat unik \(r >0\) yang memaksimalkan \(l(r, \bar{x}/r)\) dan karenanya merupakan mle yang unik untuk \(r\) dan \(β\). Selain itu, mereka juga menunjukkan bahwa jika \(\hat{σ}^2≤\bar{x}\) , maka kemungkinan binomial negatif akan didominasi oleh kemungkinan Poisson dengan \(\hat{λ} = \bar{x}\) . Dengan kata lain, distribusi Poisson memberikan kecocokan yang lebih baik terhadap data. Jaminan dalam kasus \(\hat{σ}^2>\hat{μ}\) memungkinkan kita untuk menggunakan beberapa algoritma untuk memaksimalkan \(l(r,\bar{x}/r)\) . Untuk metode alternatif dalam menghitung kemungkinan, kita perhatikan bahwa
\[\begin{array}{ll} l(r,\overline{x}/r)&=\sum_{i=1}^n \sum_{j=1}^{x_i}\log(r-1+j) - \sum_{i=1}^n\log(x_i!) \\ & \ \ \ - n(r+\overline{x}) \log(r+\overline{x}) + nr\log(r) + n\overline{x}\log(\overline{x}), \end{array}\]yang menghasilkan
\[\begin{align*} \left(\frac{1}{n}\right)\frac{\delta}{\delta r}l(r,\overline{x}/r)=\frac{1}{n}\sum_{i=1}^n \sum_{j=1}^{x_i}\frac{1}{r-1+j} - \log(r+\overline{x}) + \log(r). \end{align*}\]
Kita perhatikan bahwa, pada ekspresi di atas untuk suku-suku yang melibatkan penjumlahan ganda, jumlah bagian dalam sama dengan nol jika \(x_i = 0\). Estimasi kemungkinan maksimum untuk \(r\) adalah akar dari ekspresi terakhir dan kita dapat menggunakan algoritma pencarian akar untuk menghitungnya. Selain itu, kita juga memiliki
\[\begin{align*} \left(\frac{1}{n}\right)\frac{\delta^2}{\delta r^2}l(r,\overline{x}/r)=\frac{\overline{x}}{r(r+\overline{x})}-\frac{1}{n}\sum_{i=1}^n \sum_{j=1}^{x_i}\frac{1}{(r-1+j)^2}. \end{align*}\]
Algoritma pencarian akar berulang yang sederhana namun cepat konvergen adalah metode Newton, yang secara kebetulan diyakini telah digunakan oleh orang Babilonia untuk menghitung akar kuadrat. Dalam metode ini, sebuah perkiraan awal dipilih untuk akar dan perkiraan baru untuk akar tersebut dihasilkan secara berurutan sampai konvergen. Menerapkan metode Newton pada masalah kita akan menghasilkan algoritma berikut: Langkah i. Pilih sebuah solusi hampiran, katakanlah \(r_0\) . Tetapkan \(k\) ke \(0\). Langkah ii. Tetapkan \(r_{k+1}\) sebagai
\[\begin{align*} r_{k+1}= r_k - \frac{\frac{1}{n}\sum_{i=1}^n \sum_{j=1}^{x_i}\frac{1}{r_k-1+j} - \log(r_k+\overline{x})+\log(r_k)}{\frac{\overline{x}}{r_k(r_k+\overline{x})}-\frac{1}{n}\sum_{i=1}^n\sum_{j=1}^{x_i}\frac{1}{(r_k-1+j)^2}} \end{align*}\]
Langkah iii. Jika \(r_{k+1} ∼ r_k\) maka laporkan \(r_{k+1}\) sebagai estimasi maksimum kemungkinan; jika tidak, naikkan \(k\) sebesar \(1\) dan ulangi Langkah ii.
Sebagai contoh, kami mensimulasikan \(5\) sampel pengamatan \(41,49,40,27,23\) dari binomial negatif dengan parameter \(r = 10\) dan \(β = 5\). Memilih nilai awal dari \(r\) sedemikian sehingga
\[\begin{align*} r\beta=\hat{\mu} \quad \hbox{and} \quad r\beta(1+\beta)=\hat{\sigma}^2 \end{align*}\]
di mana \(\hat{μ}\) merupakan estimasi rata-rata dan \(\hat{σ}^{2}\) adalah estimasi varians. Hal ini menghasilkan nilai awal untuk \(r\) sebesar \(23,14286\). Iterasi dari \(r\) dari metode Newton adalah
\[\begin{align*} 21.39627, 21.60287, 21.60647, 21.60647; \end{align*}\]
konvergensi cepat yang terlihat di atas adalah tipikal dari metode Newton. Oleh karena itu, dalam contoh ini, \(\hat{r}_{MLE}∼21.60647\) dan \(\hat{β}_{MLE} = 1.66616\).
Newton<-function(x,abserr){
mu<-mean(x);
sigma2<-mean(x^2)-mu^2;
r<-mu^2/(sigma2-mu);
b<-TRUE;
iter<-0;
while (b) {
tr<-r;
m1<-mean(c(x[x==0],sapply(x[x>0],function(z){sum(1/(tr:(tr-1+z)))})));
m2<-mean(c(x[x==0],sapply(x[x>0],function(z){sum(1/(tr:(tr-1+z))^2)})));
r<-tr-(m1-log(1+mu/tr))/(mu/(tr*(tr+mu))-m2);
b<-!(abs(tr-r)<abserr);
iter<-iter+1;
}
c(r,iter)
}Untuk meringkas pembahasan kita tentang MLE untuk kelas distribusi \((a,b,0)\), pada Gambar 2.3 di bawah ini kita memplot nilai maksimum dari peluang Poisson , \(L(m,\bar{x}/m)\) untuk binomial, dan \(L(r,\bar{x}/r)\) untuk binomial negatif, untuk tiga sampel berukuran 5 yang diberikan pada Tabel 2.1. Data tersebut dibuat untuk mencakup tiga urutan dari rata-rata dan varians sampel. Seperti yang ditunjukkan pada Gambar 2.3, dan didukung oleh teori, jika (\(\hat{μ}\)≤\(\hat{σ}^2\) ) maka binomial negatif menghasilkan nilai kemungkinan maksimum yang lebih tinggi; jika \(\hat{μ}\)=\(\hat{σ}^2\) maka Poisson memiliki nilai kemungkinan tertinggi; dan terakhir dalam kasus \(\hat{μ}\)>\(\hat{σ}^2\) binomial memberikan kecocokan yang lebih baik daripada yang lain. Jadi sebelum mencocokkan data frekuensi dengan distribusi \((a,b,0)\) \((a,b,0)\), yang terbaik adalah memulai dengan memeriksa urutan \(\hat{μ}\) dan \(\hat{σ}^2\) . Sekali lagi kami tekankan bahwa Poisson berada pada batas distribusi binomial negatif dan binomial. Jadi, dalam kasus \(\hat{μ}\)≥\(\hat{σ}^2\) (\(\hat{μ}\)≤\(\hat{σ}^2\)), Poisson menghasilkan kecocokan yang lebih baik daripada binomial negatif (binomial, resp.), yang ditunjukkan oleh \(\hat{r}\) = ∞ (\(\hat{m}\)=∞ masing-masing).
\[\begin{align*} \small{ \begin{array}{c|c|c} \hline \text{Data} & \text{Mean }(\hat{\mu}) & \text{Variance }(\hat{\sigma}^2) \\ \hline (2,3,6,8,9) & 5.60 & 7.44 \\ (2,5,6,8,9) & 6 & 6\\ (4,7,8,10,11) & 8 & 6\\\hline \end{array} } \end{align*}\]
library(png)
img <- readPNG("figur3.png")
plot(2:3, 2:3, type = "n", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
rasterImage(img, 2, 2, 3, 3)