Optimization

Quasi Newton Methods

library(readr)
library(tidyr)
library(dplyr)
library(ggplot2)

Metode Quasi-Newton di R dapat diakses melalui fungsi optim(), yang merupakan fungsi optimasi tujuan umum. Fungsi optim() mengimplementasikan berbagai metode tetapi di bagian ini kita akan fokus pada metode "BFGS" dan "L-BFGS-B". Untuk contoh ini kita akan menggunakan konsentrasi rata-rata harian nitrogen dioksida (\(NO_2\)) untuk tahun 2016 yang ditemukan dalam file ini. Secara khusus, kita akan fokus pada data untuk monitor yang berlokasi di Negara Bagian Washington.

library(readr)
library(tidyr)
dat0 <- read_csv("daily_42602_2016.csv.bz2")
names(dat0) <- make.names(names(dat0))
dat <- filter(dat0, State.Name == "Washington") %>%
        unite(site, State.Code, County.Code, Site.Num) %>%
        rename(no2 = Arithmetic.Mean, date = Date.Local) %>%
        select(site, date, no2)

Estimasi kepadatan kernel dari data \(NO_2\) menunjukkan distribusi berikut.

library(ggplot2)
ggplot(dat, aes(x = no2)) + 
        geom_density()

Sebagai langkah awal dalam mengkarakterisasi distribusi nilai NO2 (dan untuk mendemonstrasikan penggunaan optim() untuk model yang pas), kita akan mencoba menyesuaikan model Normal yang terpotong ke data. Normal yang terpotong dapat masuk akal untuk jenis data ini karena mereka benar-benar positif, membuat distribusi Normal standar tidak sesuai.

Untuk normal terpotong, terpotong dari bawah di 0, densitas data adalah \[ f(x)=\frac{\frac{1}{\sigma}\varphi\left(\frac{x-\mu}{\delta}\right)}{\int^\infty_0 \frac{1}{\sigma}\varphi\left(\frac{x-\mu}{\delta}\right)dx} \]

Parameter yang tidak diketahui adalah \(\mu\) dan \(\sigma\). Mengingat kepadatan, kita dapat mencoba untuk memperkirakan dan dengan kemungkinan maksimum. Dalam hal ini, kita akan meminimalkan negative log-likehood dari data.

Kita dapat menggunakan fungsi deriv() untuk menghitung kemungkinan log negatif dan gradiennya secara otomatis. Karena kita menggunakan metode quasi-Newton di sini kita tidak memerlukan matriks Hessian.

nll_one <- deriv(~ -log(dnorm((x - mu)/s) / s) + log(0.5),
                 c("mu", "s"), 
                 function.arg = TRUE)

Fungsi optim() bekerja sedikit berbeda dari nlm() karena alih-alih memiliki gradien sebagai atribut kemungkinan log negatif, gradien harus menjadi fungsi yang terpisah.

Pertama kemungkinan log negatif.

nll <- function(p) {
        v <- nll_one(p[1], p[2])
        sum(v)
}

Fungsi gradien

nll_grad <- function(p) {
        v <- nll_one(p[1], p[2])
        colSums(attr(v, "gradient"))
}

Sekarang kita dapat meneruskan fungsi nll() dan nll_grad() ke optim() untuk mendapatkan estimasi dan . kita akan menggunakan nilai awal \(\mu=1\) dan \(\sigma=5\). Untuk menggunakan metode quasi-Newton "BFGS", Anda perlu menentukannya dalam argumen metode. Metode default untuk optim() adalah metode simpleks Nelder-Mead. kita juga menentukan hessian = TRUE untuk memberi tahu optim() untuk menghitung secara numerik matriks Hessian pada titik optimal.

x <- dat$no2
res <- optim(c(1, 5), nll, gr = nll_grad, 
             method = "BFGS", hessian = TRUE)
res

## $par
## [1] 13.23731  8.26315
## 
## $value
## [1] 4043.641
## 
## $counts
## function gradient 
##       35       19 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
##              [,1]         [,2]
## [1,] 2.087005e+01 5.659674e-04
## [2,] 5.659674e-04 4.174582e+01

Fungsi optim() mengembalikan daftar dengan 5 elemen (ditambah matriks Hessian jika hessian = TRUE diatur). Elemen pertama yang harus Anda periksa adalah kode konvergensi. Jika konvergensi adalah 0, itu bagus. Apa pun selain 0 dapat menunjukkan masalah, yang sifatnya tergantung pada algoritme yang Anda gunakan (lihat halaman bantuan untuk optim() untuk detail lebih lanjut). Kali ini kita juga memiliki optim() menghitung Hessian (secara numerik) pada titik optimal sehingga kita dapat memperoleh kesalahan standar asimtotik jika kita mau.

Catatan pertama bahwa ada beberapa pesan yang dicetak ke konsol saat algoritme sedang berjalan yang menunjukkan bahwa NaN dihasilkan oleh fungsi target. Ini kemungkinan karena fungsi tersebut mencoba mengambil log angka negatif. Karena kami menggunakan algoritme "BFGS", kami melakukan pengoptimalan tanpa kendala. Oleh karena itu, kemungkinan pencarian algoritme menghasilkan nilai negatif untuk \(\sigma\), yang tidak masuk akal dalam konteks ini. Untuk membatasi pencarian, kita dapat menggunakan metode "L-BFGS-B" yang merupakan algoritma BFGS “memori terbatas” dengan “kendala kotak”. Ini memungkinkan Anda untuk menempatkan batas bawah dan atas pada setiap parameter dalam model.

Perhatikan bahwa optim() memungkinkan fungsi target Anda untuk menghasilkan nilai NA atau NaN, dan memang dari output tampaknya algoritme akhirnya bertemu pada jawabannya. Tetapi karena kita tahu bahwa parameter dalam model ini dibatasi, kita dapat melanjutkan dan menggunakan pendekatan alternatif.

Di sini kita menetapkan batas bawah untuk semua parameter menjadi 0 tetapi membiarkan batas atas menjadi tak terhingga (Inf), yang merupakan default.

res <- optim(c(1, 5), nll, gr = nll_grad, 
             method = "L-BFGS-B", hessian = TRUE,
             lower = 0)
res

## $par
## [1] 13.237470  8.263546
## 
## $value
## [1] 4043.641
## 
## $counts
## function gradient 
##       14       14 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
## 
## $hessian
##              [,1]         [,2]
## [1,] 20.868057205 -0.000250771
## [2,] -0.000250771 41.735838073

Kita dapat melihat sekarang bahwa pesan peringatan hilang, tetapi solusinya identik dengan yang dihasilkan oleh metode "BFGS" asli.

Estimasi kemungkinan maksimum dari \(\mu\) adalah 13,24 dan taksiran dari \(\sigma\) adalah 8.26. Jika kita ingin mendapatkan kesalahan standar asimtotik untuk parameter ini, kita dapat melihat matriks Hessian.

solve(res$hessian) %>%
        diag %>%
        sqrt

## [1] 0.2189067 0.1547909

Namun dalam hal ini, kami tidak terlalu peduli dengan kesalahan standar sehingga kami akan melanjutkan.

Kita dapat memplot kerapatan asli dari data versus model Normal terpotong yang dipasang untuk melihat seberapa baik kita mengkarakterisasi distribusi. Pertama kita akan mengevaluasi model yang dipasang pada 100 titik antara 0 dan 50.

xpts <- seq(0, 50, len = 100)
dens <- data.frame(xpts = xpts,
                   ypts = dnorm(xpts, res$par[1], res$par[2]))

Kemudian kita dapat melapisi model yang dipasang di atas kepadatan menggunakan geom_line().

ggplot(dat, aes(x = no2)) + 
        geom_density() + 
        geom_line(aes(x = xpts, y = ypts), data = dens, col = "steelblue",
                  lty = 2)

Ini tidak cocok. Melihat kepadatan halus data, jelas bahwa ada dua mode pada data, menunjukkan bahwa Normal yang terpotong mungkin tidak cukup untuk mengkarakterisasi data.

Salah satu alternatif dalam kasus ini adalah campuran dari dua Normal, yang mungkin menangkap dua mode. Untuk campuran dua komponen, densitas datanya adalah

\[ \lambda\frac{1}{\sigma}\varphi\left(\frac{x-\mu_1}{\sigma_1}\right)+(1-\lambda)\frac{1}{\sigma}\varphi\left(\frac{x-\mu_2}{\sigma_2}\right) \]

Umumnya, kita melihat bahwa model ini cocok menggunakan algoritma yang lebih kompleks seperti algoritma EM atau metode rantai Markov Monte Carlo. Meskipun metode tersebut memberikan stabilitas yang lebih besar dalam proses estimasi (seperti yang akan kita lihat nanti), kita sebenarnya dapat menggunakan metode tipe Newton untuk memaksimalkan kemungkinan secara langsung dengan sedikit perhatian.

Pertama kita dapat menuliskan kemungkinan log negatif secara simbolis dan mengizinkan fungsi R deriv() untuk menghitung fungsi gradien.

nll_one <- deriv(~ -log(lambda * dnorm((x-mu1)/s1)/s1 + (1-lambda)*dnorm((x-mu2)/s2)/s2), 
                 c("mu1", "mu2", "s1", "s2", "lambda"), 
                 function.arg = TRUE)

Kemudian, seperti sebelumnya, kita dapat menentukan fungsi negative log-likelihood (nll) dan gradien R (nll_grad) secara terpisah.

nll <- function(p) {
        p <- as.list(p)
        v <- do.call("nll_one", p)
        sum(v)
}
nll_grad <- function(p) {
        v <- do.call("nll_one", as.list(p))
        colSums(attr(v, "gradient"))
}

Terakhir, kita dapat meneruskan fungsi-fungsi tersebut ke optim() dengan vektor parameter awal. Di sini, kami berhati-hati untuk menentukan

• Kami menggunakan metode "L-BFGS-B" sehingga kami menentukan batas bawah 0 untuk semua parameter dan batas atas 1 untuk parameter \(\lambda\)
• Kami mengatur opsi parscale dalam daftar parameter kontrol, yang mirip dengan argumen typsize ke nlm(). Tujuannya di sini adalah untuk memberikan optim() skala untuk setiap parameter di sekitar titik optimal.

x <- dat$no2
pstart <- c(5, 10, 2, 3, 0.5)
res <- optim(pstart, nll, gr = nll_grad, method = "L-BFGS-B",
             control = list(parscale = c(2, 2, 1, 1, 0.1)),
             lower = 0, upper = c(Inf, Inf, Inf, Inf, 1))
res

## $par
## [1]  3.7606598 16.1469811  1.6419640  7.2378153  0.2348927
## 
## $value
## [1] 4879.924
## 
## $counts
## function gradient 
##       17       17 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Kode konvergensi 0 adalah pertanda baik dan estimasi parameter dalam vektor par semuanya tampak masuk akal. Kita dapat melapisi model yang dipasang ke kepadatan yang halus untuk melihat bagaimana modelnya.

xpts <- seq(0, 50, len = 100)
dens <- with(res, {
        data.frame(xpts = xpts, 
                   ypts = par[5]*dnorm(xpts, par[1], par[3]) + (1-par[5])*dnorm(xpts, par[2], par[4]))
})
ggplot(dat, aes(x = no2)) + 
        geom_density() + 
        geom_line(aes(x = xpts, y = ypts), data = dens, col = "steelblue",
                  lty = 2)

Kesesuaiannya masih kurang bagus, tetapi setidaknya model ini menangkap secara kasar lokasi dua mode dalam kepadatan. Juga, tampaknya model menangkap ekor kerapatan dengan cukup baik, meskipun ini perlu diperiksa lebih hati-hati dengan melihat kuantil.

Terakhir, seperti kebanyakan model dan skema pengoptimalan, biasanya ide yang baik untuk memvariasikan titik awal untuk melihat apakah perkiraan kami saat ini adalah mode lokal.

pstart <- c(1, 20, 5, 2, 0.1)
res <- optim(pstart, nll, gr = nll_grad, method = "L-BFGS-B",
             control = list(parscale = c(2, 2, 1, 1, 0.1)),
             lower = 0, upper = c(Inf, Inf, Inf, Inf, 1))
res

## $par
## [1]  3.760571 16.146834  1.641961  7.237776  0.234892
## 
## $value
## [1] 4879.924
## 
## $counts
## function gradient 
##       22       22 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Di sini kita melihat bahwa dengan titik awal yang sedikit berbeda, kita mendapatkan nilai yang sama dan kemungkinan log negatif minimum yang sama.