8  Distribusi Asimtotik MLE (Bagian II)

8.1.1 Selang Kepercayaan Normal Asimtotik untuk Data Berdistribusi t

Sebagai ilustrasi, perhatikan data berikut dari distribusi t. Di sini kita akan menyimulasikan data:

y=c(-0.71,-1.30,-0.13,-2.03, 1.62, 2.38, 0.48, 0.51,-0.69,
    -2.32,-2.02,1.23,-0.25,0.76,0.65,-0.08,-1.20,0.99,
    2.58,0.73,-0.09,-0.03,0.56,-1.44,0.13)
hist(y)

Data ini saling bebas dan berasal dari \(y_i \sim t(4)\), yaitu distribusi t dengan parameter derajat kebebasan sebesar 4. Jika kita tidak mengetahui nilai sebenarnya dari parameter derajat kebebasan (\(\text{df}\)), mula-mula kita dapat mencoba mencari nilai dugaan kemungkinan maksimumnya secara langsung. Namun, PDF distribusi t tidak menghasilkan MLE yang dapat diselesaikan secara analitik. Kita dapat menemukannya dari data dengan pendugaan kemungkinan maksimum numerik. Inferensi numerik di R dapat dilakukan dengan menulis fungsi negatif log-kemungkinan, lalu mengoptimisasikannya menggunakan optim.

## fungsi negatif log-kemungkinan dengan df dibatasi positif
nll.t=function(df,y){
  if(length(df)!=1L || !is.finite(df) || df<=0) return(Inf)
  value=-sum(dt(y,df=df,log=TRUE))
  if(is.finite(value)) value else Inf
}
## MLE numerik berbatas; simpan dan periksa diagnostik pengoptimal
out=optim(4,fn=nll.t,y=y,method="L-BFGS-B",
          lower=sqrt(.Machine$double.eps),hessian=TRUE)
if(out$convergence!=0L || !is.finite(out$value) ||
   any(!is.finite(out$par))) stop("Pendugaan MLE t gagal")
df.hat=unname(out$par)
df.hat

[1] 6.63125

Terlihat bahwa MLE kita untuk \(\text{df}\) adalah \(\hat{\text{df}}=6.63125\).

Selanjutnya, kita dapat membentuk selang kepercayaan dengan menggunakan hampiran Normal asimtotik untuk MLE sebagai berikut:

## informasi teramati: Hessian negatif log-kemungkinan pada nilai dugaan
J.n=drop(out$hessian)
if(length(J.n)!=1L || !is.finite(J.n) || J.n<=0)
  stop("Informasi teramati tidak positif dan berhingga")
## selang Wald asimtotik dengan informasi teramati plug-in
c(df.hat-1.96*sqrt(1/J.n),df.hat+1.96*sqrt(1/J.n))

[1] -6.6074 19.8699

Perhatikan bahwa selang kepercayaan ini cukup lebar dan bahkan memuat nilai negatif, yang mustahil bagi parameter derajat kebebasan dalam distribusi t!

Selang kepercayaan ini mengandalkan hampiran distribusi Normal yang diketahui menjadi hampiran akurat ketika ukuran sampel \(n\) menuju tak hingga. Namun, untuk setiap ukuran sampel tetap (\(n=25\) dalam contoh ini), kita tidak mengetahui seberapa baik hampiran tersebut. Karena himpunan dukungan distribusi Normal adalah seluruh garis bilangan real, nilai negatif wajar muncul dalam distribusi hampiran itu. Karena parameter \(\text{df}\) harus positif, kita mengetahui bahwa distribusi hampiran ini mungkin kurang sesuai untuk situasi sekarang.

Pada bagian berikutnya, kami menyajikan cara lain untuk membentuk selang kepercayaan hampiran menggunakan sampel bootstrap. Kita akan membahas bootstrap parametrik dan bootstrap nonparametrik.

8.1.2 Bootstrap Parametrik

Jenis bootstrap pertama yang akan kita pelajari adalah bootstrap parametrik. Dalam setiap metode bootstrap, kita memerlukan cara untuk menyimulasikan data dari hampiran terhadap distribusi yang menghasilkan data tersebut.

Dalam bootstrap parametrik, pertama-tama kita menduga parameter model data, lalu menggunakan nilai dugaan parameter itu sebagai parameter dalam model data untuk menyimulasikan data.

Misalkan \(x_i\sim f_X(x_i|\theta)\), \(i=1, \ldots, n\), merupakan titik-titik data yang saling bebas dan diasumsikan berasal dari distribusi \(f_X(x|\theta)\) yang bergantung pada parameter tak diketahui \(\theta\).

Algoritme bootstrap parametrik adalah sebagai berikut:

1. Peroleh Nilai Dugaan \(\hat{\theta}\)
   Dengan menggunakan data \(\{x_1,x_2,\ldots,x_n\}\), dugalah parameter tersebut dengan penduga \(\hat{\theta}=h(\{x_1,x_2,\ldots,x_n\})\).

Pendugaan ini dapat dilakukan dengan MLE analitik, MLE numerik, atau bahkan pendekatan lain seperti metode momen.

2. Simulasikan Sampel Bootstrap
   Dengan menggunakan nilai dugaan parameter \(\hat{\theta}\) dari langkah 1, simulasikan \(M\) himpunan data yang saling bebas.

Setiap himpunan data berukuran sama, yaitu \(n\) seperti data asli, dan masing-masing berasal dari \(x_i^{(m)}\sim f_X(x|\hat{\theta})\).

3. Peroleh Nilai Dugaan \(\{\hat{\theta}^{(m)}\}_{m=1}^{M}\) dengan Data Bootstrap
   Sekarang gunakan masing-masing dari \(M\) himpunan data tersebut untuk memperoleh nilai dugaan bagi \(\theta\), dengan metode yang sama seperti pada langkah 1.

Hasilnya adalah \(M\) nilai dugaan parameter yang berbeda: \(\{\hat{\theta}^{(1)},\hat{\theta}^{(2)},\ldots,\hat{\theta}^{(M)}\}\), masing-masing berasal dari himpunan data simulasi yang berbeda \[ \{x_1^{(1)},x_2^{(1)},\ldots,x_n^{(1)}\} \rightarrow \hat{\theta}^{(1)}\] \[ \{x_1^{(2)},x_2^{(2)},\ldots,x_n^{(2)}\} \rightarrow \hat{\theta}^{(2)}\] \[ \vdots \] \[ \{x_1^{(M)},x_2^{(M)},\ldots,x_n^{(M)}\} \rightarrow \hat{\theta}^{(M)} \] Kumpulan \(M\) nilai dugaan parameter tersebut kira-kira berdistribusi seperti sampel-sampel dari distribusi nilai dugaan kita \(\hat{\theta}\); hampiran itu timbul karena kita menyimulasikan data dari \(f_X(x|\hat{\theta})\) —distribusi pembangkit data dengan mengasumsikan bahwa parameter sebenarnya sama dengan nilai dugaan kita \(\hat{\theta}\), bukan dari distribusi pembangkit yang sebenarnya \(f_X(x|\theta)\).

4. Inferensi terhadap \(\hat{\theta}\)
   Selanjutnya, kumpulan nilai dugaan parameter ini dapat digunakan untuk menghampiri secara numerik sifat-sifat distribusi nilai dugaan kita \(\hat{\theta}\).

Sebagai contoh, jika kita tertarik pada selang kepercayaan 95% untuk \(\theta\), kita dapat menghampirinya dengan selang probabilitas 95% dari distribusi sampling ini, yakni dengan mencari kuantil 0,025 dan 0,975 dari sampel-sampel \(\{\hat{\theta}^{(1)},\hat{\theta}^{(2)},\ldots,\hat{\theta}^{(M)}\}\). \[\text{SK }95\%\text{ untuk }\theta\approx\left(\hat q_{.025}\{\hat\theta^{(1)},\ldots,\hat\theta^{(M)}\},\ \hat q_{.975}\{\hat\theta^{(1)},\ldots,\hat\theta^{(M)}\}\right)\] Untuk menghampiri kuantil 0,025 dan 0,975 dengan baik, kita memerlukan banyak \(M\) sampel bootstrap.

x=c(-0.71,-1.30,-0.13,-2.03, 1.62, 2.38, 0.48, 0.51,-0.69,
    -2.32,-2.02,1.23,-0.25,0.76,0.65,-0.08,-1.20,0.99,
    2.58,0.73,-0.09,-0.03,0.56,-1.44,0.13)

8.1.3 Selang Kepercayaan Bootstrap Nonparametrik

Bootstrap parametrik yang diperkenalkan di atas menyimulasikan himpunan data baru dari distribusi yang (secara hampiran) sama dengan distribusi data teramati kita, dengan menggunakan model yang diasumsikan (distribusi t dalam contoh di atas) dan parameter yang diduga dari data teramati (MLE untuk \(\text{df}\) dalam contoh di atas).

Pendekatan alternatif untuk bootstrap adalah menyimulasikan himpunan data baru langsung dari data teramati melalui pengambilan sampel ulang (resampling). Pendekatan ini disebut bootstrap nonparametrik.

Penalaran di balik bootstrap nonparametrik adalah sebagai berikut. Jika kita telah mengamati data \(x_i \sim f_X(x)\) yang saling bebas dan berdistribusi identik, maka informasi terbaik yang kita miliki tentang distribusi \(f_X\) sebenarnya adalah data itu sendiri. Kita dapat menghampiri satu pengamatan baru dari \(f_X\) hanya dengan memilih secara acak salah satu nilai teramati \(x_1,x_2,\ldots,x_n\). Hal ini disebut “pengambilan sampel dari distribusi empiris data”.

Distribusi empiris dari suatu himpunan data \(x_1,x_2,\ldots,x_n\) adalah distribusi diskret dengan himpunan dukungan berupa nilai-nilai unik dalam data \(x_1,x_2,\ldots,x_n\) dan PMF yang didefinisikan oleh \[P_n(X=x)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}\{x_i=x\}\]

Dengan demikian, distribusi empiris memberikan bobot 1/n kepada setiap pengamatan; suatu nilai yang berulang memperoleh massa sebesar frekuensi relatifnya dalam data.

x=c(-0.71,-1.30,-0.13,-2.03, 1.62, 2.38, 0.48, 0.51,-0.69,
    -2.32,-2.02,1.23,-0.25,0.76,0.65,-0.08,-1.20,0.99,
    2.58,0.73,-0.09,-0.03,0.56,-1.44,0.13)

8.1.4 Kode Templat untuk Selang Kepercayaan Bootstrap

Untuk mempermudah penerapan metode bootstrap, berikut seperangkat kode templat beserta beberapa komentar.

## baca data
x=c(1,2,3)

## (1) tentukan MLE; ganti ... dan batas dengan spesifikasi model
nll=function(theta,x){
  ## kode negatif log-kemungkinan model ditempatkan di sini
}
out=optim(...,nll,x=x,method="L-BFGS-B",lower=...,upper=...)
if(out$convergence!=0L || !is.finite(out$value) ||
   any(!is.finite(out$par))) stop("Pendugaan data asli gagal")
theta.hat=out$par
## siapkan bootstrap yang dapat direproduksi
RNGversion("4.3.0")
RNGkind("Mersenne-Twister","Inversion","Rejection")
set.seed(4150803)
M=1000L
n=length(x)
sim.data=vector("list",M)
theta.hat.vals=rep(NA_real_,M)
convergence=rep(NA_integer_,M)

for(m in seq_len(M)){
  ## (2) ganti rexp(...) dengan pembangkitan dari model yang dipasang
  sim.data[[m]]=rexp(n,theta.hat)
  ## (3) pasang model yang sama pada data bootstrap
  out.sim=tryCatch(
    optim(...,nll,x=sim.data[[m]],method="L-BFGS-B",lower=...,upper=...),
    error=function(e) NULL
  )
  if(is.null(out.sim)) next
  convergence[m]=out.sim$convergence
  if(out.sim$convergence==0L && is.finite(out.sim$value) &&
     all(is.finite(out.sim$par))) theta.hat.vals[m]=out.sim$par
}
failed=is.na(convergence) | convergence!=0L | !is.finite(theta.hat.vals)
if(any(failed)) stop(sprintf("%d pendugaan bootstrap gagal",sum(failed)))
quantile(theta.hat.vals,c(.025,.975))
sessionInfo()

## baca data
x=c(1,2,3)

## (1) tentukan MLE; ganti ... dan batas dengan spesifikasi model
nll=function(theta,x){
  ## kode negatif log-kemungkinan model ditempatkan di sini
}
out=optim(...,nll,x=x,method="L-BFGS-B",lower=...,upper=...)
if(out$convergence!=0L || !is.finite(out$value) ||
   any(!is.finite(out$par))) stop("Pendugaan data asli gagal")
theta.hat=out$par
## siapkan bootstrap yang dapat direproduksi
RNGversion("4.3.0")
RNGkind("Mersenne-Twister","Inversion","Rejection")
set.seed(4150804)
M=1000L
n=length(x)
sim.data=vector("list",M)
theta.hat.vals=rep(NA_real_,M)
convergence=rep(NA_integer_,M)

for(m in seq_len(M)){
  ## (2) resampel dari distribusi empiris
  sim.data[[m]]=sample(x,size=n,replace=TRUE)
  ## (3) pasang model yang sama pada data bootstrap
  out.sim=tryCatch(
    optim(...,nll,x=sim.data[[m]],method="L-BFGS-B",lower=...,upper=...),
    error=function(e) NULL
  )
  if(is.null(out.sim)) next
  convergence[m]=out.sim$convergence
  if(out.sim$convergence==0L && is.finite(out.sim$value) &&
     all(is.finite(out.sim$par))) theta.hat.vals[m]=out.sim$par
}
failed=is.na(convergence) | convergence!=0L | !is.finite(theta.hat.vals)
if(any(failed)) stop(sprintf("%d pendugaan bootstrap gagal",sum(failed)))
quantile(theta.hat.vals,c(.025,.975))
sessionInfo()

8.1.5 Contoh Dua Parameter dengan Data Berdistribusi Pareto

## load the library needed for "dpareto"
library(EnvStats)
## read in data
x=c(5.92, 10.87,  6.50, 14.62, 20.49,  5.12,  7.28, 15.24,
   7.47,  6.78, 24.07, 6.76,  8.81,  7.65,  5.28, 15.80,
   5.76,  5.11,  6.10, 23.44, 15.04,  9.02,  8.34, 66.05,
   8.52,  9.26,  7.40,  7.85,  5.93,  5.41, 26.00, 16.00,
   8.99, 11.06,  5.06,  6.92, 10.17,  5.65,  6.06,  5.70)
hist(x)

## MLE analitik Pareto-I; domain dan kegagalan diperiksa eksplisit
pareto.mle=function(z){
  if(!length(z) || any(!is.finite(z)) || any(z<=0))
    stop("Data Pareto harus positif dan berhingga")
  L.hat=min(z)
  denominator=sum(log(z/L.hat))
  if(!is.finite(denominator) || denominator<=0)
    stop("Pendugaan bentuk Pareto gagal")
  a.hat=length(z)/denominator
  if(!is.finite(a.hat) || a.hat<=0) stop("Pendugaan bentuk Pareto gagal")
  c(L=L.hat,a=a.hat)
}
theta.hat=pareto.mle(x)
L.hat=unname(theta.hat["L"])
a.hat=unname(theta.hat["a"])

RNGversion("4.3.0")
RNGkind("Mersenne-Twister","Inversion","Rejection")
set.seed(4150805)
M=1000L
n=length(x)
L.hat.vals=rep(NA_real_,M)
a.hat.vals=rep(NA_real_,M)
for(m in seq_len(M)){
  sim.data=sample(x,size=n,replace=TRUE)
  fit=tryCatch(pareto.mle(sim.data),error=function(e) NULL)
  if(is.null(fit)) next
  L.hat.vals[m]=fit["L"]
  a.hat.vals[m]=fit["a"]
}
failed=!is.finite(L.hat.vals) | !is.finite(a.hat.vals)
if(any(failed)) stop(sprintf("%d pendugaan bootstrap gagal",sum(failed)))

## Hanya kontraexample: kuantil ini bukan selang kepercayaan sah untuk L
quantile(L.hat.vals,c(.025,.975))
sessionInfo()

 2.5% 97.5% 
 5.06  5.28 

## (4) find a 95% CI for a using the bootstrap samples
quantile(a.hat.vals,c(.025,.975))

    2.5%    97.5% 
1.277910 2.270643 

8.2.1 Metode Delta

Jika kita menduga parameter \(\theta\) dengan pendugaan kemungkinan maksimum, maka—di bawah konsistensi, identifiabilitas, parameter benar yang interior, diferensiabilitas, pertukaran turunan dan integral yang sah, serta informasi tak singular—distribusi limit MLE \(\hat{\theta}\) adalah \[\sqrt{n}(\hat{\theta}_n-\theta_0)\xrightarrow{d}N\!\left(0,I_1(\theta_0)^{-1}\right)\] dengan \(I_n(\theta)\) merupakan informasi Fisher harapan untuk sampel, yang didefinisikan sebagai \[I_n(\theta)=-\operatorname{E}_{\theta}\!\left[\ell_n''(\theta;X)\right]\] dengan \(\ell(\theta)=log(L(\theta))\) dan nilai harapannya diambil terhadap semua \(x_1,x_2,\ldots,x_n\).

Syarat koreksi: hasil limit ini memerlukan model terdominasi dan dapat diidentifikasi, parameter benar yang interior, konsistensi, diferensiabilitas dan pertukaran integral-turunan yang sah, serta informasi Fisher yang positif dan taksingular. I₁ adalah informasi harapan per pengamatan; Jₙ pada nilai dugaan hanya merupakan besaran plug-in untuk galat baku, bukan varians acak dalam hukum limit.

Jika \(g(\theta)\) terdiferensialkan pada parameter benar—tanpa perlu dapat dibalik—maka pendekatan umum untuk inferensi adalah metode delta. Metode ini memakai distribusi limit Gaussian MLE untuk memperoleh hampiran bagi distribusi penduga parameter yang ditransformasi \(\hat{\tau}=g(\hat{\theta})\). Ketika \(n\rightarrow \infty\), \[\sqrt{n}\{g(\hat{\theta}_n)-g(\theta_0)\}\xrightarrow{d}N\!\left(0,[g'(\theta_0)]^2I_1(\theta_0)^{-1}\right)\] dan dengan demikian, selang kepercayaan asimtotik hampiran 95% bagi parameter yang ditransformasi \(\tau=g(\theta)\) adalah \[g(\hat{\theta}_n)\pm1.96\frac{|g'(\hat{\theta}_n)|}{\sqrt{J_n(\hat{\theta}_n)}}\]

Metode delta orde pertama hanya memerlukan g terdiferensialkan pada parameter benar. Turunan dalam varians limit dievaluasi pada nilai benar; turunan pada nilai dugaan dipakai hanya sebagai plug-in. Jika turunannya nol, limit orde pertama degenerat dan metode orde lebih tinggi mungkin diperlukan.

8.2.2 Bootstrap dan Transformasi

Ketika bootstrap digunakan untuk mengambil sampel dari distribusi penarikan sampel bagi nilai dugaan parameter, hasilnya berupa sejumlah besar nilai, yaitu \(M\) nilai dugaan bootstrap \(\hat{\theta}^{(1)},\hat{\theta}^{(2)},\ldots,\hat{\theta}^{(M)}\). Jika kita ingin melakukan inferensi terhadap transformasi \(\theta\), kita cukup mentransformasi nilai-nilai dugaan ini secara langsung lalu melakukan inferensi terhadap sampel yang dihasilkan. Artinya, jika kita ingin melakukan inferensi terhadap \(\tau=g(\theta)\), kita melakukan langkah-langkah berikut:

1. Dapatkan \(M\) nilai dugaan bootstrap sebagaimana ditunjukkan pada Bagian 8.1. Hasilnya adalah \(\hat{\theta}^{(1)},\hat{\theta}^{(2)},\ldots,\hat{\theta}^{(M)}\)
2. Transformasikan masing-masing dari \(M\) nilai dugaan tersebut sehingga diperoleh \(g(\hat{\theta}^{(1)}),g(\hat{\theta}^{(2)}),\ldots,g(\hat{\theta}^{(M)})\)
3. Kemudian, inferensi—misalnya menentukan selang kepercayaan—dapat dilakukan secara langsung pada nilai-nilai dugaan bootstrap yang telah ditransformasi ini. Sebagai contoh, selang kepercayaan bootstrap 95% bagi \(\tau=g(\theta)\) ditentukan oleh kuantil .025 dan .975 dari \(\{g(\hat{\theta}^{(1)}),g(\hat{\theta}^{(2)}),\ldots,g(\hat{\theta}^{(M)})\}\)

Bootstrap dan pendekatan berbasis sampel lainnya sering menjadi pilihan utama untuk melakukan inferensi terhadap transformasi parameter karena alasan berikut:

• Menerapkan transformasi pada setiap replikasi bootstrap tidak menuntut fungsi \(g\) dapat dibalik. Namun, keabsahan inferensi bootstrap tetap bergantung pada statistik, transformasi, skema resampling, desain sampel, dan syarat keteraturan; bootstrap parametrik juga mengasumsikan keluarga distribusi yang dipasang. Untuk metode delta orde pertama, \(g\) harus terdiferensialkan pada parameter benar.

• Jika model memiliki lebih dari satu parameter (yaitu \(\boldsymbol\theta=(\theta_1,\theta_2)\)) dan transformasi \(g\) merupakan fungsi dari beberapa parameter: \(\tau=g(\theta_1,\theta_2)\), pendekatan bootstrap sangat mudah diterapkan karena kita hanya perlu menerapkan fungsi \(g\) pada nilai dugaan bootstrap bagi \(\theta_1\) dan \(\theta_2\).