Contoh 1

Mencari nilai maksimum dari \(sin(x) + sin(2 * x) + cos(3 * x)\) dengan interval \([0,2\pi]\)**

f1 <- function(x) sin(x) + sin(2 * x) + cos(3 * x)
curve(f1, from = 0, to = 2 * pi) 

Permasalahan tersebut adalah model optimasi satu peubah, sebaiknya diselesaikan dengan fungsi optimize yang sudah tertanam dalam R untuk mendapatkan nilai minimum dari suatu fungsi dengan satu peubah:

optimize(f1, interval = c(0, 2*pi))              # minimum lokal

## $minimum
## [1] 3.033129
## 
## $objective
## [1] -1.054505

optimize(f1, interval = c(4, 2*pi))              # minimum global

## $minimum
## [1] 5.273383
## 
## $objective
## [1] -2.741405

optimize(f1, interval = c(0, 2*pi), maximum = T) # maksimum lokal

## $maximum
## [1] 4.0598
## 
## $objective
## [1] 1.096473

optimize(f1, interval = c(0, 1.5), maximum = T)  # maksimum global

## $maximum
## [1] 0.3323289
## 
## $objective
## [1] 1.485871

Permasalahan optimasi dengan fungsi satu variabel sudah terselesaikan dengan baik. Berikut ini akan dibahas bagaimana kita mencari nilai minimum dari fungsi lebih dari satu peubah.

Contoh 2:

Mencari nilai \(x_1\) dan \(x2\) yang pada fungsi \[f(x_1,x_2)=100(x_2−x_1^2)^2+(1−x_1)^2\]

f2 <- function (x) {
  x1 <- x[1]
  x2 <- x[2]
  100*(x2 - x1^2)^2 + (1 - x1)^2
}

optim(c( -3 ,3) , f2)

## $par
## [1] 1.006766 1.013557
## 
## $value
## [1] 4.582775e-05
## 
## $counts
## function gradient 
##      193       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Algoritma Steepest Descent

• Metode ini dimulai pada tebakan awal \(x\).
• Tebakan berikutnya dibuat dengan bergerak ke arah gradien negatif. Lokasi minimum sepanjang garis ini kemudian dapat ditemukan dengan menggunakan algoritma pencarian satu dimensi seperti metode Golden Section.
• Pembaruan ke-n kemudian

\[x_n=x_{n−1}−\alpha f'(x_{n−1})\]

di mana \(\alpha\) dipilih untuk meminimalkan fungsi satu dimensi:

\[g(\alpha)=f(x_{n−1}−\alpha f'(x_{n−1}))\]

Sehingga dalam hal ini, kita perlu mengasumsikan bahwa \(\alpha\) berada dalam suatu interval. Jadi dalam hal ini, kita ambil intervalnya menjadi \([0,h]\) dimana \(h\) adalah nilai yang harus kita pilih. Berikut ini adalah Algoritma Steepest Descent di dalam R.

steepestdescent <- function(f, fprime, start, h,
 tol=1e-7, maxiter=100) {
 x <- start
 g <- function(alpha) { f(x - alpha*fpx) }
 niter <- 0
 while(niter < maxiter & sum(abs(fprime(x))) > tol) {
 fpx <- fprime(x)
 alpha <- golden(g, 0, h)
 x <- x - alpha*fpx
 niter <- niter + 1
 }
 if (niter == maxiter) {
 print("Warning: Maximum number of iterations reached")
 }
 c("Minimizer" = x)
 }

Proses optimasi Steepest Descent secara visualisasi dapat dilihat pada grafik berikut, atau klik di sini untuk melihat contoh lain dari metode Steepest Descent ini, khususnya untuk fungsi kuadrat. selengkapnya

Steepest Descent Dua Dimensi

Berikut ini adalah contoh penerapan Steepest Descent untuk melakukan proses optimasi minimal dengan fungsi 2 variabel.

\[f(x)=f(x_1,x_2)={(2−x_1)^2\over 2x_2^2} + {(3−x_1)^2 \over 2x_2^2}+log(x_2)\]

Perhatikan bahwa \(x_2\) harus positif, jadi kita mungkin perlu melindungi dari nilai negatif \(x_2\).

Pertama: Kita membutuhkan fungsi untuk fungsi dan gradien

f1 <- function(x) {
 (2-x[1])^2/(2*x[2]^2) +(3-x[1])^2/(2*x[2]^2) + log(x[2])
 }
 f1prime <- function(x) {
 c(-(2-x[1])/x[2]^2 - (3-x[1])/x[2]^2,-(2-x[1])^2/x[2]^3 -(3-x[1])^2/x[2]^3 + 1/x[2])
 }

Kedua: Lokasi minimum di sepanjang garis ini kemudian dapat ditemukan dengan menggunakan algoritma pencarian satu dimensi seperti metode Golden Section.

golden <- function (f, a, b, tol = 0.0000001)
 {
 ratio <- 2 / (sqrt(5) + 1)
 x1 <- b - ratio * (b - a)
 x2 <- a + ratio * (b - a)
 f1 <- f(x1)
 f2 <- f(x2)

 while(abs(b - a) > tol) {
 if (f2 > f1) {
 b <- x2
 x2 <- x1
 f2 <- f1
 x1 <- b - ratio * (b - a)
 f1 <- f(x1)
 } else {
 a <- x1
 x1 <- x2
 f1 <- f2
 x2 <- a + ratio * (b - a)
 f2 <- f(x2)
 }
 }
 return((a + b) / 2)
 }

Ketiga: Mari kita coba nilai awal x=(.1,.1).

steepestdescent(f1, f1prime, start=c(.1,.1), h=.1)

## [1] "Warning: Maximum number of iterations reached"

## Minimizer1 Minimizer2 
##  2.4992967  0.7123675

Pada posisi ini kita belum menemukan nilai konvergennya, sehingga dilajutkan pada langkah berikutnya.

Keempat: Dalam hal ini kita kembali ke iterasa pertama dengan menggunakan hasil terbaru sebagai tebakan awal.

steepestdescent(f1, f1prime,start=c(2.4992967, 0.7123675), h=.1)

## Minimizer1 Minimizer2 
##  2.5000000  0.7071068

steepestdescent(f1, f1prime, start=c(.1, .1), h=.1,maxiter=200)

## Minimizer1 Minimizer2 
##  2.5000000  0.7071068

Selesai. Nilainya sudah konvergen.

Multivariat Normal

Kita dapat menggunakan Steepest Descent untuk menghitung perkiraan kemungkinan maksimum rata-rata dalam distibusi multivariate normal, yang diberikan sampel data. Namun, ketika data sangat berkorelasi, seperti pada contoh simulasi di bawah, permukaan kemungkinan log bisa menjadi sulit untuk dioptimalkan. Dalam kasus seperti itu, ridge yang sangat sempit berkembang dalam kemungkinan log yang dapat menyulitkan algoritma Steepest Descent untuk dinavigasi.

Pada contoh di bawah, kita sebenarnya menghitung kemungkinan log negatif karena algoritme dirancang untuk meminimalkan fungsi.

set.seed(2020-09-30)
mu <- c(1, 2)
S <- rbind(c(1, .9), c(.9, 1))
x <- MASS::mvrnorm(500, mu, S)
nloglike <- function(mu1, mu2) {
        dmv <- mvtnorm::dmvnorm(x, c(mu1, mu2), S, log = TRUE)
        -sum(dmv)
}
nloglike <- Vectorize(nloglike, c("mu1", "mu2"))
nx <- 40
ny <- 40
xg <- seq(-5, 5, len = nx)
yg <- seq(-5, 6, len = ny)
g <- expand.grid(xg, yg)
nLL <- nloglike(g[, 1], g[, 2])
z <- matrix(nLL, nx, ny)
par(mar = c(4.5, 4.5, 1, 1))
contour(xg, yg, z, nlevels = 40, xlab = expression(mu[1]), 
        ylab = expression(mu[2]))
abline(h = 0, v = 0, lty = 2)

Perhatikan bahwa pada gambar di atas permukaan sangat membentang dan minimum (1,2) terletak di tengah lembah sempit. Untuk algoritma Steepest Descent kita akan mulai pada titik (−5,−2) dan melacak jalur algoritma.

library(dplyr, warn.conflicts = FALSE)
norm <- function(x) x / sqrt(sum(x^2))
Sinv <- solve(S)                          
step1 <- function(mu, alpha = 1) {
        D <- sweep(x, 2, mu, "-")
        score <- colSums(D) %>% norm
        mu + alpha * drop(Sinv %*% score)
}
steep <- function(mu, n = 10, ...) {
        results <- vector("list", length = n)
        for(i in seq_len(n)) {
                results[[i]] <- step1(mu, ...)
                mu <- results[[i]]
        }
        results
}
m <- do.call("rbind", steep(c(-5, -2), 8))
m <- rbind(c(-5, -2), m)

par(mar = c(4.5, 4.5, 1, 1))
contour(xg, yg, z, nlevels = 40, xlab = expression(mu[1]), 
        ylab = expression(mu[2]))
abline(h = 0, v = 0, lty = 2)
points(m, pch = 20, type = "b")

Kita dapat melihat bahwa jalur algoritma agak berliku karena melintasi lembah yang sempit. Dalam hal ini, kita telah memperbaiki panjang langkah dalam kasus ini, yang mungkin saja tidak optimal. Namun, kita sudah dapat melihat bahwa algoritma memiliki beberapa kesulitan dalam menavigasi permukaan karena arah penurunan paling curam tidak selalu mengarah ke arah minimum.

Model Linier Umum

Model linier umum (Generalized Linear Model) adalah perluasan dari model linier sederhana yang memungkinkan obeservasi berdistribusi tidak ormal. Distribusi yang digunakan biasanya berasal dari keluarga eksponensial yang fungsi densitasnya memiliki beberapa karakteristik yang sama. Dengan GLM, kita biasanya menyajikannya sebagai \(y_i \sim p(y_i|\mu_i)\), di mana \(p\) adalah distribusi keluarga eksponensial, \(\mathbb{E}[y_i]=\mu_i\),

\[g(\mu_i)=x'_i\beta,\]

di mana \(g\) adalah fungsi nonlinier, dan \(Var(y_i)=V(\mu)\) dengan \(V\) adalah fungsi varians yang diketahui.

Berbeda dengan model linier sederhana, estimasi kemungkinan maksimum dari parameter vektor \(\beta\) tidak dapat diperoleh dalam bentuk tertutup, sehingga algoritma iteratif harus digunakan untuk mendapatkan estimasi. Algoritma klasik yang digunakan adalah algoritma penilaian Fisher. Algoritma ini menggunakan pendekatan linier untuk fungsi nonlinier \(g\), yang dapat ditulis sebagai:

\[g(y_i)\approx g(\mu_i)+(y_i-\mu_i)g'(\mu_i).\]

Notasi khas GLM mengacu pada \(z_i=g(\mu_i)+(y_i-\mu_i)g'(\mu_i)\) sebagai variabel respon. Algoritma penilaian Fisher kemudian bekerja sebagai berikut:

1. Mulailah dengan \(\hat{\mu_i}\), beberapa nilai awal.
2. Hitung \(z_i=g(\mu_i)+(y_i-\mu_i)g'(\mu_i)\).
3. Mengingat \(n \times 1\) adalah vektor variabel respon \(z\) dan matriks prediktor \(n \times p\) \(X\) kita dapat menghitung regresi berbobot \(z\) pada \(X\) untuk mendapatkan \[\beta_n=(X'WX)^{-1}X'Wz\] di mana \(W\) adalah matriks diagonal dengan elemen diagonal \[w_{ii}=[g'(\mu_i)^2V(\mu_i)]^{-1}\].
4. Mengingat \(\beta_n,\) kita dapat menghitung ulang \(\hat{\mu_i}=g^{-1}(x'_i\beta_n)\) dan pergi ke 2.

Perhatikan bahwa pada Langkah 3 di atas, bobot hanyalah kebalikan dari varians \(z_i\), yaitu.

\[\begin {align*} Var(z_i) & = Var(g(\mu)+(y_i-\mu_i)g'(\mu_i))\\ & =Var((y_i-\mu_i)g'(\mu_i))\\ & =V(\mu_i)g'(\mu_i)^2 \end{align*}\]

Secara alami, ketika melakukan regresi berbobot, kita akan menggunakan fungsi kebalikan dari variansinya.

Contoh: Regresi Poisson

Untuk regresi Poisson, kita memiliki \(y_i\sim Poisson (\mu_i)\) dimana \(g(\mu)=log \space \mu_i=x'_i\beta\) karena log adalah fungsi tautan kanonik untuk distribusi Poisson, kita juga memiliki \(g'(\mu_i)=\frac{1}{\mu_i}\) dan \(V(\mu_i)=\mu_i.\) Oleh karena itu, algoritma penskoran Fisher adalah 1. Inisialisasi \(\hat{\mu_i},\) mungkin menggunakan \(y_i+1\) (untuk menghindari nol). 2. Biarkan \(z_i=log \space \hat{\mu_i}+(y_i-\hat{\mu_i})\frac{1}{\hat{\mu_i}}\) 3. Regresi \(z\) pada \(X\) menggunakan bobot \[w_{ii}=[\frac{1}{\hat{\mu_i^2}}\hat{\mu_i}]^{-1}=\hat{\mu_i}.\]

Dengan menggunakan contoh regresi Poisson, kita dapat menarik hubungan antara algoritma penilaian Fisher yang biasa untuk pemasangan GLM dan metode Newton. Ingat bahwa jika \(l(\beta)\) adalah kemungkinan log sebagai fungsi dari parameter regresi \(\beta\), maka skema pembaruan Newton adalah \[\beta_{n+1}=\beta_n+l''(\beta_n)^{-1}[-l'(\beta_n)].\]

Kemungkinan log untuk model regresi Poisson dapat ditulis dalam bentuk vektor/matriks sebagai \[l(\beta)=y'X\beta- exp(X\beta)'1\] di mana eksponensial diambil berdasarkan komponen pada vektor \(X\beta\). Fungsi gradiennya adalah \[l'(\beta)=X'y-X'x(X\beta)=X'(y-\mu)\] dan Hessian adalah \[l''(\beta)=-X'WX\] di mana \(W\) adalah matriks diagonal dengan nilai \(w_{ii}=exp(x'_i\beta)\) pada diagonalnya. Iterasi Newton adalah

\[\begin {eqnarray*} \beta_{n+1}=\beta_n+(-X'WX)^{-1}(-X'(y-\mu))\\ =\beta_n+(X'WX)^{-1}XW(z-X\beta_n)\\ =(X'WX)^{-1}X'Wz+\beta_n-(X'WX)^{-1}X'WX\beta_n\\ =(X'WX)^{-1}X'Wz \end {eqnarray*}\]

Oleh karena itu iterasinya persis sama dengan algoritma penilaian Fisher dalam kasus ini. Secara umum, metode Newton dan penskoran Fisher akan bertepatan dengan model linier umum apa pun yang menggunakan keluarga eksponensial dengan fungsi tautan kanonik.

Arah Newton dengan R

Fungsi nlm() dalam R mengimplementasikan metode Newton untuk meminimalkan fungsi yang diberikan vektor yang menyatakan nilai. Secara default, seseorang tidak perlu menyediakan fungsi gradien atau Hessian; mereka akan diestimasi secara numerik oleh algoritma. Namun, untuk tujuan meningkatkan akurasi algoritme, gradien dan Hessian dapat diberikan sebagai atribut dari fungsi target.

Sebagai contoh, kita akan menggunakan fungsi nlm() agar sesuai dengan model regresi logistik sederhana untuk data biner. Model ini menetapkan bahwa \(y_i \sim \space Bernoulli(p_i)\) where \[log \frac{p_i}{1-p_i}=\beta_0+x_i\beta_1\]

dan tujuannya adalah untuk memperkirakan \(\beta\) melalui kemungkinan maksimum. Mengingat distribusi Bernoulli yang diasumsikan, kita dapat menulis kemungkinan log untuk satu pengamatan sebagai:

\[\begin{eqnarray*} \log L(\beta) & = & \log\left\{\prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i}\right\}\\ & = & \sum_{i=1}^n y_i\log p_i + (1-y_i)\log(1-p_i)\\ & = & \sum_{i=1}^n y_i\log\frac{p_i}{1-p_i}+\log(1-p_i)\\ & = & \sum_{i=1}^n y_i(\beta_0 + x_i\beta_1) + \log\left(\frac{1}{1+e^{(\beta_0 + x_i\beta_1)}}\right)\\ & = & \sum_{i=1}^n y_i(\beta_0 + x_i\beta_1) -\log\left(1+e^{(\beta_0 + x_i\beta_1)}\right) \end{eqnarray*}\]

Jika kita mengambil baris terakhir dari turunan di atas dan mengambil satu elemen di dalam jumlah, kita memiliki

\[\ell_i(\beta)=y_i(\beta_0 + x_i\beta_1) -\log\left(1+e^{(\beta_0 + x_i\beta_1)}\right)\]

Kita akan membutuhkan gradien dan Hessian untuk \(\beta\). Karena jumlah dan turunannya dapat dipertukarkan, kita kemudian dapat menjumlahkan masing-masing gradien individu dan Hessian untuk mendapatkan gradien penuh dan Hessian untuk seluruh sampel, sehingga \[\ell^\prime(\beta) = \sum_{i=1}^n\ell_i^\prime(\beta)\] dan \[\ell^\prime(\beta) = \sum_{i=1}^n\ell_i^{\prime\prime}(\beta)\] Sekarang, mengambil gradien dan Hessian dari ekspresi di atas mungkin agak merepotkan. Namun, R menyediakan cara otomatis untuk melakukan diferensiasi simbolik sehingga pekerjaan manual dapat dihindari. Fungsi deriv() menghitung gradien dan Hessian dari ekspresi secara simbolis sehingga dapat digunakan dalam rutinitas minimalisasi. Itu tidak dapat menghitung gradien ekspresi arbitrer, tetapi mendukung berbagai fungsi statistik umum.

Contoh: Tren nilai-p Sepanjang Waktu

Paket tidypvals yang ditulis oleh Jeff Leek berisi kumpulan data yang diambil dari literatur yang mengumpulkan nilai-p yang terkait dengan berbagai publikasi bersama dengan beberapa informasi tentang publikasi tersebut (yaitu jurnal, tahun, DOI). Satu pertanyaan yang muncul adalah apakah ada tren dari waktu ke waktu dalam signifikansi statistik yang diklaim dari publikasi, di mana “signifikansi statistik” didefinisikan memiliki nilai \(p\) kurang dari 0,05.

Paket tidypvals tersedia dari GitHub dan dapat diinstal menggunakan fungsi install_github() dalam paket remotes.

remotes::install_github("jtleek/tidypvals")

Setelah terinstal, kami akan menggunakan kumpulan data jager2014. Secara khusus, kita tertarik untuk membuat indikator apakah nilai \(p\) kurang dari 0,05 dan melakukan regresi pada variabel tahun.

# devtools::install_github('jtleek/tidypvals')
library(tidypvals)
library(dplyr)
jager <- mutate(tidypvals::jager2014, 
pvalue = as.numeric(as.character(pvalue)),
     y = ifelse(pvalue < 0.05 |(pvalue == 0.05 & operator == "lessthan"), 1, 0),
     x = year - 2000) %>% tbl_df
jager

Perhatikan di sini bahwa kita telah mengurangi tahun 2000 dari variabel tahun sehingga \(x=0\) sesuai dengan tahun==2000.

Selanjutnya kita menghitung gradien dan Hessian dari kemungkinan log negatif sehubungan dengan \(\beta_0\) dan \(\beta_1\) menggunakan fungsi deriv(). Di bawah ini, kami menetapkan function.arg=TRUE dalam panggilan ke deriv() karena kami ingin deriv() mengembalikan fungsi yang argumennya adalah b0 dan b1.

nll_one <- deriv(~ -(y * (b0 + x * b1) - log(1 + exp(b0 + b1 * x))),
             c("b0", "b1"), function.arg = TRUE, hessian = TRUE)

Berikut tampilan fungsinya.

nll_one

## function (b0, b1) 
## {
##     .expr6 <- exp(b0 + b1 * x)
##     .expr7 <- 1 + .expr6
##     .expr11 <- .expr6/.expr7
##     .expr15 <- .expr7^2
##     .expr18 <- .expr6 * x
##     .expr19 <- .expr18/.expr7
##     .value <- -(y * (b0 + x * b1) - log(.expr7))
##     .grad <- array(0, c(length(.value), 2L), list(NULL, c("b0", 
##         "b1")))
##     .hessian <- array(0, c(length(.value), 2L, 2L), list(NULL, 
##         c("b0", "b1"), c("b0", "b1")))
##     .grad[, "b0"] <- -(y - .expr11)
##     .hessian[, "b0", "b0"] <- .expr11 - .expr6 * .expr6/.expr15
##     .hessian[, "b0", "b1"] <- .hessian[, "b1", "b0"] <- .expr19 - 
##         .expr6 * .expr18/.expr15
##     .grad[, "b1"] <- -(y * x - .expr19)
##     .hessian[, "b1", "b1"] <- .expr18 * x/.expr7 - .expr18 * 
##         .expr18/.expr15
##     attr(.value, "gradient") <- .grad
##     attr(.value, "hessian") <- .hessian
##     .value
## }

Fungsi nll_one() yang dihasilkan oleh deriv() mengevaluasi kemungkinan log negatif untuk setiap titik data. Output dari nll_one() akan memiliki atribut "gradient" dan "hessian" yang masing-masing mewakili gradien dan Hessian. Misalnya, dengan menggunakan data dari kumpulan data jager, kita dapat mengevaluasi kemungkinan log negatif pada \(\beta_0=0,\beta_1=0.\)

x <- jager$x
y <- jager$y
str(nll_one(0, 0))

##  num [1:15653] 0.693 0.693 0.693 0.693 0.693 ...
##  - attr(*, "gradient")= num [1:15653, 1:2] -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "b0" "b1"
##  - attr(*, "hessian")= num [1:15653, 1:2, 1:2] 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "b0" "b1"
##   .. ..$ : chr [1:2] "b0" "b1"

Fungsi nll_one() mengevaluasi kemungkinan log negatif pada setiap titik data, tetapi tidak menjumlahkan poin seperti yang diperlukan untuk mengevaluasi kemungkinan log negatif penuh. Oleh karena itu, akan dituliskan fungsi terpisah yang dapat melakukan itu untuk kemungkinan log negatif, gradien, dan Hessian.

nll <- function(b) {
        v <- nll_one(b[1], b[2])
        f <- sum(v)                                     ## log-likelihood
        gr <- colSums(attr(v, "gradient"))              ## gradient vector
        hess <- apply(attr(v, "hessian"), c(2, 3), sum) ## Hessian matrix
        attributes(f) <- list(gradient = gr, 
                              hessian = hess)
        f
}

Sekarang, kita dapat mengevaluasi kemungkinan log negatif penuh dengan fungsi nll(). Perhatikan bahwa nll() mengambil satu vektor numerik sebagai input karena inilah yang diharapkan oleh fungsi nlm().

nll(c(0, 0))

## [1] 10849.83
## attr(,"gradient")
##       b0       b1 
##  -4586.5 -21854.5 
## attr(,"hessian")
##          b0        b1
## b0  3913.25  19618.25
## b1 19618.25 137733.75

Menggunakan \(\beta_0=0,\beta_1=0\) sebagai nilai awal, kita dapat memanggil nlm() untuk meminimalkan kemungkinan log negatif.

res <- nlm(nll, c(0, 0))
res

## $minimum
## [1] 7956.976
## 
## $estimate
## [1]  1.57032807 -0.04416515
## 
## $gradient
## [1] -1.451746e-06 -2.782241e-06
## 
## $code
## [1] 1
## 
## $iterations
## [1] 4

Perhatikan terlebih dahulu pada output bahwa ada kode dengan nilai 4 dan jumlah iterasinya adalah 100. Setiap kali jumlah iterasi dalam algoritme pengoptimalan adalah bilangan bulat yang bagus, kemungkinan besar ia memiliki beberapa batas iterasi yang telah ditetapkan. Ini pada gilirannya biasanya berarti algoritme tidak konvergen.

Dalam bantuan untuk nlm() kita juga mempelajari bahwa nilai kode dari 4 berarti “batas iterasi terlampaui”, yang umumnya tidak baik. Untungnya, solusinya sederhana: kita dapat meningkatkan batas iterasi dan membiarkan algoritme berjalan lebih lama.

res <- nlm(nll, c(0, 0), iterlim = 1000)
res

## $minimum
## [1] 7956.976
## 
## $estimate
## [1]  1.57032807 -0.04416515
## 
## $gradient
## [1] -1.451746e-06 -2.782241e-06
## 
## $code
## [1] 1
## 
## $iterations
## [1] 4

Dalam hal ini kita melihat bahwa jumlah iterasi yang digunakan adalah 260, yang jauh di bawah batas iterasi. Sekarang kita mendapatkan kode sama dengan 2 yang berarti bahwa “iterasi berturut-turut dalam toleransi, iterasi saat ini mungkin solusi”. Itu terdengar seperti kabar baik!

Terakhir, sebagian besar algoritme pengoptimalan memiliki opsi untuk menskalakan nilai parameter Anda sehingga nilainya bervariasi secara kasar pada skala yang sama. Jika fungsi target Anda memiliki parameter yang bervariasi pada skala yang sangat berbeda, ini dapat menyebabkan masalah praktis bagi komputer (ini bukan masalah bagi teori). Cara untuk mengatasi hal ini di nlm() adalah dengan menggunakan argumen typsize, yang merupakan vektor yang panjangnya sama dengan vektor parameter yang menyediakan ukuran relatif parameter.

Di sini, saya memberikan typsize = c(1, 0.1), yang menunjukkan kepada nlm() bahwa parameter pertama, \(\beta_0\), kira-kira 10 kali lebih besar dari parameter kedua, \(\beta_1\) ketika fungsi target minimum.

res <- nlm(nll, c(0, 0), iterlim = 1000,
           typsize = c(1, 0.1))
res

## $minimum
## [1] 7956.976
## 
## $estimate
## [1]  1.57032807 -0.04416515
## 
## $gradient
## [1] -1.451745e-06 -2.782238e-06
## 
## $code
## [1] 1
## 
## $iterations
## [1] 4

Dengan menjalankan nlm() anda akan melihat bahwa solusinya sama tetapi jumlah iterasi sebenarnya jauh lebih sedikit dari sebelumnya (4 iterasi) yang berarti algoritme berjalan lebih cepat. Secara umum, penskalaan vektor parameter dengan tepat (jika mungkin) meningkatkan kinerja semua algoritme pengoptimalan dan menurut pengalaman saya, hampir selalu merupakan ide yang bagus. Nilai spesifik yang diberikan pada argumen typsize tidak penting; melainkan hubungan mereka satu sama lain (yaitu urutan besarnya) yang penting.