R 语言数据分析实战 - 27 优化问题

[giscus[bot]]

R 语言数据分析实战 - 27  优化问题

https://bookdown.org/xiangyun/data-analysis-in-action/optimization-problems.html

[XiangyunHuang]

高斯过程回归模型的参数估计，REML 限制极大似然估计

library(MASS)
data(topo)

log_lik <- function(x) {
  n <- nrow(topo)
  D <- t(t(rep(1, n)))
  Sigma <- x[1]^2 * exp(-as.matrix(dist(topo[, c("x", "y")])) / x[2])
  inv_Sigma <- solve(Sigma)
  P <- diag(1, n) - D %*% solve(t(D) %*% solve(Sigma, D), t(D)) %*% inv_Sigma
  as.vector(-1 / 2 * log(det(Sigma)) - 1 / 2 * log(det(t(D) %*% solve(Sigma, D))) - 1 / 2 * t(topo[, "z"]) %*% inv_Sigma %*% P %*% topo[, "z"])
}

log_lik(x = c(65, 2))

Sys.setenv(ROI_LOAD_PLUGINS = FALSE)
library(ROI)
library(ROI.plugin.nloptr)

op <- OP(
  objective = F_objective(log_lik, n = 2L),
  bounds = V_bound(lb = c(100, 10), ub = c(150, 50)),
  maximum = TRUE
)

nlp <- ROI_solve(op, solver = "nloptr.directL")
# sigma phi
nlp$solution
# 目标函数值
nlp$objval
# REML 对数似然
nlp$objval - 52 / 2 * log(2 * pi)

library(lattice)

dat <- expand.grid(
  sigma = seq(from = 100, to = 150, length.out = 41),
  phi = seq(from = 10, to = 50, length.out = 31)
)
dat$fn <- apply(dat, 1, log_lik)

levelplot(fn ~ sigma * phi,
  data = dat, aspect = 1,
  xlim = c(100, 150), ylim = c(10, 50),
  xlab = expression(sigma), ylab = expression(phi),
  col.regions = cm.colors, contour = TRUE,
  scales = list(
    x = list(alternating = 1, tck = c(1, 0)),
    y = list(alternating = 1, tck = c(1, 0))
  )
)