3.5 Maximum Likelihood Estimation

Hingga saat ini, bab ini berfokus pada distribusi parametrik yang biasa digunakan dalam aplikasi asuransi. Namun, agar berguna dalam pekerjaan terapan, distribusi ini harus menggunakan nilai “realistis” untuk parameter dan untuk ini kita beralih ke data. Pada tingkat dasar, kami berasumsi bahwa analis telah menyediakan sampel acak \(X_1,... , X_n\) dari distribusi dengan fungsi distribusi \(F_x\) (untuk singkatnya, kita terkadang menjatuhkan subskrip \(X\)). Vektor \(θ\) untuk menunjukkan kumpulan parameter untuk \(F\).

Sebelum menggambar dari distribusi, kami mempertimbangkan hasil potensial yang dirangkum oleh variabel acak \(X_i\) (\(i\) = 1,2,..,n). Dari Bab sebelumnya yang membahas perkiraan distribusi frekuensi, dimana \(Pr(X_1 = x_1,...,X_n = x_n)\) untuk mengukur “kemungkinan” menggambar sampel \({x_1,...,x_n}\) Dengan data kontinu, kami menggunakan fungsi kerapatan probabilitas bersama alih-alih probabilitas bersama. Dengan asumsi independensi, pdf bersama dapat ditulis sebagai produk pdf. Dengan demikian, kami mendefinisikan kemungkinannya menjadi

\[\begin{equation} L(\boldsymbol \theta) = \prod_{i=1}^n f(x_i) . \end{equation}\]

Dari notasi, perhatikan bahwa kami menganggap ini sebagai fungsi dari parameter di \(θ\), dengan data \({x_1,...,x_n}\) iadakan tetap. Penaksir kemungkinan maksimum adalah nilai parameter di \(θ\) yang memaksimalkan \(L(θ)\).

Dari kalkulus, kita tahu bahwa memaksimalkan fungsi menghasilkan hasil yang sama dengan memaksimalkan logaritma suatu fungsi (ini karena logaritma adalah fungsi monoton). Karena kita mendapatkan hasil yang sama, untuk memudahkan pertimbangan komputasi, adalah umum untuk mempertimbangkan kemungkinan logaritmik, yang dilambangkan sebagai

\[\begin{equation} l(\boldsymbol \theta) = \log L(\boldsymbol \theta) = \sum_{i=1}^n \log f(x_i) . \end{equation}\]

Contoh 3.5.1. Soal Ujian Aktuaria. Anda diberikan lima pengamatan berikut: 521, 658, 702, 819, 1217. Anda menggunakan Pareto parameter tunggal dengan fungsi distribusi: \[F(x) = 1- \left(\frac{500}{x}\right)^{\alpha}, ~~~~ x>500 .\]

Dengan \(n = 5\), fungsi log-likelihood adalah

\[l(\alpha) = \sum_{i=1}^5 \log f(x_i;\alpha ) = 5 \alpha \log 500 + 5 \log \alpha-(\alpha+1) \sum_{i=1}^5 \log x_i.\]

Gambar dibawah menunjukkan kemungkinan logaritmik sebagai fungsi parameter α.

Kita dapat menentukan nilai maksimum kemungkinan logaritmik dengan mengambil turunan dan mengaturnya sama dengan nol. Ini menghasilkan

Secara alami, ada banyak masalah di mana tidak praktis menggunakan perhitungan tangan untuk optimasi. Untungnya ada banyak rutinitas statistik yang tersedia seperti fungsi .R optim

c1 <- log(521)+log(658)+log(702)+log(819)+log(1217)
nloglike <- function(alpha){-(5*alpha*log(500)+5*log(alpha)-(alpha+1)*c1)}
MLE <- optim(par=1, fn=nloglike)$par

## Warning in optim(par = 1, fn = nloglike): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly

MLE

## [1] 2.453125

2.453125 mengkonfirmasi hasil perhitungan tangan kita di mana penaksir kemungkinan maksimum

Contoh 3.5.4. Dana Properti Wisconsin. Untuk melihat bagaimana penaksir kemungkinan maksimum bekerja dengan data nyata, kami kembali ke data klaim 2010

Cuplikan kode berikut menunjukkan cara menyesuaikan eksponensial, gamma, Pareto, lognormal, dan \(GB2\) Model. Untuk konsistensi, kode menggunakan package R dan VGAM. Akronim adalah singkatan dari Vector Generalized Linear and Additive Models; Seperti yang disarankan oleh namanya, paket ini dapat melakukan jauh lebih dari sesuai dengan model-model ini meskipun cukup untuk tujuan kita. Satu-satunya pengecualian adalah GB2 kepadatan yang tidak banyak digunakan di luar aplikasi asuransi; Namun, kita dapat mengkodekan kepadatan ini dan menghitung penaksir kemungkinan maksimum menggunakan optim General Purpose Optimizer.

library(VGAM)

## Loading required package: stats4

## Loading required package: splines

claim_lev <- read.csv("CLAIMLEVEL.csv", header = TRUE) 
claim_data <- subset(claim_lev, Year == 2010); 

# Inference assuming a GB2 Distribution - this is more complicated
# The likelihood function of GB2 distribution (negative for optimization)
lik_gb2 <- function (param) {
  a_1 <- param[1]
  a_2 <- param[2]
  mu <- param[3]
  sigma <- param[4]
  yt <- (log(claim_data$Claim) - mu) / sigma
  logexpyt <- ifelse(yt > 23, yt, log(1 + exp(yt)))
  logdens <- a_1 * yt - log(sigma) - log(beta(a_1,a_2)) - 
    (a_1+a_2) * logexpyt - log(claim_data$Claim) 
  return(-sum(logdens))
}
# "optim" is a general purpose minimization function
gb2_bop <- optim(c(1, 1, 0, 1), lik_gb2, method = c("L-BFGS-B"), 
                 lower = c(0.01, 0.01, -500, 0.01), 
                 upper = c(500, 500, 500, 500), hessian = TRUE)
# Nonparametric Plot
plot(density(log(claim_data$Claim)), main = "", xlab = "Log Expenditures",
     ylim = c(0 ,0.37))
x <- seq(0, 15, by = 0.01)
#Exponential
fit.exp <- vglm(Claim ~ 1, exponential, data = claim_data)

## Warning in vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2, :
## iterations terminated because half-step sizes are very small

## Warning in vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2, : some
## quantities such as z, residuals, SEs may be inaccurate due to convergence at a
## half-step

theta = 1 / exp(coef(fit.exp))
fexp_ex <- dgamma(exp(x), scale = exp(-coef(fit.exp)), shape = 1) * exp(x)
lines(x, fexp_ex, col = "red", lty =2)
# Inference assuming a gamma distribution
fit.gamma <- vglm(Claim ~ 1, family = gamma2, data = claim_data)
theta <- exp(coef(fit.gamma)[1]) / exp(coef(fit.gamma)[2])  # theta = mu / alpha
alpha <- exp(coef(fit.gamma)[2]) 
fgamma_ex <- dgamma(exp(x), shape = alpha, scale = theta) * exp(x)
lines(x, fgamma_ex, col = "blue", lty =3)
#Pareto
fit.pareto <- vglm(Claim ~ 1, paretoII, loc = 0, data = claim_data)
fpareto_ex <- dparetoII(exp(x), loc = 0, shape = exp(coef(fit.pareto)[2]), 
                        scale = exp(coef(fit.pareto)[1])) * exp(x)
lines(x, fpareto_ex, col = "purple")
# Lognormal
fit.LN <- vglm(Claim ~ 1, family = lognormal, data = claim_data)
flnorm_ex <- dlnorm(exp(x), mean = coef(fit.LN)[1],
                    sd = exp(coef(fit.LN)[2])) * exp(x)
lines(x, flnorm_ex, col = "lightblue")
# Density for GB II
gb2_density <- function (x) {
  a_1 <- gb2_bop$par[1]
  a_2 <- gb2_bop$par[2]
  mu <- gb2_bop$par[3]
  sigma <- gb2_bop$par[4]
  xt <- (log(x) - mu) / sigma
  logexpxt <- ifelse (xt > 23, yt, log(1 + exp(xt)))
  logdens <- a_1 * xt - log(sigma) - log(beta(a_1, a_2)) - 
    (a_1+a_2) * logexpxt -log(x) 
  exp(logdens)
  }
fGB2_ex = gb2_density(exp(x)) * exp(x)
lines(x, fGB2_ex, col="green")
legend("topleft", c("log(Expend)", "Exponential", "Gamma", "Pareto", 
                    "Lognormal", "GB2"), cex=0.8,
       lty = c(4,2,3,1,1,1), #4 is "longdash"
       col = c("black","red","blue","purple","lightblue","green"))