application retel_r added

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: retel
 Type: Package
 Title: Regularized Exponentially Tilted Empirical Likelihood
-Version: 0.0.0.9702
+Version: 0.0.0.9703
 Authors@R: c(
     person("Eunseop", "Kim", email = "markean@pm.me", 
             role = c("aut", "cre")),

retel-paper/application/retel_r.txt

@@ -0,0 +1,38 @@
+< Method = RETEL_r>
+< MCMC runs = 1000000 >
+
+theta psrf:
+ [1] 1.001059 1.000586 1.001367 1.007783 1.004875 1.003851 1.001367 1.001441
+ [9] 1.005139 1.002440 1.000639 1.004450 1.002942 1.003306 1.005548 1.002568
+[17] 1.003495 1.000593 1.005675 1.005214 1.000807 1.002109 1.002791 1.001512
+[25] 1.002791 1.009413 1.001343 1.004359 1.001531 1.001163 1.004050 1.002345
+[33] 1.005026 1.003266 1.004194 1.006318 1.000942 1.004174 1.002949 1.011339
+[41] 1.002836 1.001378 1.003321 1.001866 1.001987 1.002733 1.004184 1.003738
+[49] 1.003413 1.002805 1.002670
+theta mpsrf: 1.045526
+beta psrf: 1.000172 1.000165
+beta mpsrf: 1.00033
+s2 psrf: 1.00668
+Acceptance rate:  0.129885
+
+< beta & s2 summary >
+beta1 mean:    0.01456249
+beta1 95% ci:  -0.1072708 0.1344021
+beta2 mean:    0.9333335
+beta2 95% ci:  0.812334 1.054236
+s2 mean:       0.9604577
+s2 95% ci:     0.5066496 1.649349
+
+< theta summary >
+theta 95% ci al:  3.781057
+theta aad:        0.2724186
+theta aard:       0.9002461
+theta asd:        0.1086403
+theta asrd:       3.831403
+
+< theta_os summary >
+theta_os 95% ci al:  20876.54
+theta_os aad:        1504.119
+theta_os aard:       0.0387572
+theta_os asd:        3311934
+theta_os asrd:       0.002296179

retel-paper/code/application/retel_r.R

@@ -0,0 +1,204 @@
+## 1. Load packages
+options(warn = -1)
+options(scipen = 999)
+suppressMessages(library(actuar))
+library(coda)
+suppressMessages(library(here))
+suppressMessages(i_am("code/application/retel_r.R"))
+library(mvtnorm)
+library(retel)
+
+
+## 2. Load data
+data("income", package = "retel")
+raw_income <- read.csv(paste0(dirname(getwd()), "/data-raw/income.csv"))
+
+
+## 3. Constants
+# Model
+fit <- lm(mi_1989 ~ -1 + mi_1979 + ami, income)
+b0 <- coef(fit)
+x <- model.matrix(fit)
+y <- as.numeric(model.response(fit$model))
+m <- nrow(x)
+p <- ncol(x)
+# MCMC parameters
+B <- 250000L
+
+
+## 4. Functions
+f <- function(x, par) {
+  cbind(x - par, (x - par)^2L - 1)
+}
+# Log prior density functions
+log_pd_theta <- function(theta, b, s2) {
+  dmvnorm(theta, mean = x %*% b, sigma = s2 * diag(m), log = TRUE)
+}
+log_pd_beta <- function(b, s2) {
+  dmvnorm(b, mean = b0, sigma = 0.1 * s2 * solve(crossprod(x)), log = TRUE)
+}
+log_pd_s2 <- function(s2) {
+  dinvgamma(s2, shape = 4, scale = 1, log = TRUE)
+}
+log_pd <- function(theta, b, s2) {
+  log_pd_theta(theta, b, s2) + log_pd_beta(b, s2) - log(s2)
+}
+log_posterior_unnormalized <- function(theta, b, s2, mu, sigma) {
+  g <- f(y, theta)
+  if (missing(mu)) {
+    mu <- colMeans(g)
+  }
+  if (missing(sigma)) {
+    sigma <- var(g)
+  }
+  log_pd(theta, b, s2) +
+    retel(f, y, theta, mu, sigma, log(m), "reduced")
+}
+# MCMC
+mcmc_fn <- function(B) {
+  theta_sample <- matrix(nrow = B, ncol = m)
+  theta_sample[1L, ] <- rnorm(m, mean = mean(y), sd = 0.1)
+  beta_sample <- matrix(nrow = B, ncol = p)
+  beta_sample[1L, ] <- rnorm(p, mean = b0, sd = 0.1)
+  s2_sample <- vector("numeric", length = B)
+  s2_sample[1L] <- rinvgamma(1L, shape = 4, scale = 1)
+  acceptace <- vector("logical", length = B)
+  acceptace[1L] <- FALSE
+  for (i in seq_len(B)[-1]) {
+    # Sample proposal value
+    s2_p <- rnorm(1L, mean = s2_sample[i - 1L], sd = 0.06)
+    if (s2_p < 0) {
+      s2_p <- 0
+    }
+    b_p <- rmvnorm(1L,
+      mean = beta_sample[i - 1L, ], sigma = 0.01 * diag(p)
+    ) |>
+      as.vector()
+    theta_p <- rmvnorm(1L,
+      mean = theta_sample[i - 1L, ], sigma = 0.06 * diag(m)
+    ) |>
+      as.vector()
+    # Compute log ratio of unnormailzed posterior densities
+    g <- cbind(y - theta_p, (y - theta_p)^2L - 1)
+    logr <- log_posterior_unnormalized(theta_p, b_p, s2_p) -
+      log_posterior_unnormalized(
+        theta_sample[i - 1L, ], beta_sample[i - 1L, ],
+        s2_sample[i - 1L], colMeans(g), var(g)
+      )
+    # Sample uniform random variable
+    u <- runif(1L)
+    # Accept or reject
+    if (isTRUE(log(u) < logr)) {
+      theta_sample[i, ] <- theta_p
+      beta_sample[i, ] <- b_p
+      s2_sample[i] <- s2_p
+      acceptace[i] <- TRUE
+    } else {
+      theta_sample[i, ] <- theta_sample[i - 1L, ]
+      beta_sample[i, ] <- beta_sample[i - 1L, ]
+      s2_sample[i] <- s2_sample[i - 1L]
+      acceptace[i] <- FALSE
+    }
+  }
+  list(
+    theta = theta_sample, beta = beta_sample, s2 = s2_sample, rate = acceptace
+  )
+}
+# Metrics
+aad <- function(y, e) {
+  mean(abs(y - e))
+}
+aard <- function(y, e) {
+  mean(abs((y - e) / y))
+}
+asd <- function(y, e) {
+  mean((y - e)^2L)
+}
+asrd <- function(y, e) {
+  mean(((y - e) / y)^2L)
+}
+
+
+## 5. Simulations
+set.seed(2536345)
+cat("< Method = RETEL_r>\n")
+cat("< MCMC runs =", 4L * B, ">\n")
+c1 <- mcmc_fn(B)
+c2 <- mcmc_fn(B)
+c3 <- mcmc_fn(B)
+c4 <- mcmc_fn(B)
+theta_c1 <- mcmc(c1$theta)
+theta_c2 <- mcmc(c2$theta)
+theta_c3 <- mcmc(c3$theta)
+theta_c4 <- mcmc(c4$theta)
+beta_c1 <- mcmc(c1$beta)
+beta_c2 <- mcmc(c2$beta)
+beta_c3 <- mcmc(c3$beta)
+beta_c4 <- mcmc(c4$beta)
+s2_c1 <- mcmc(c1$s2)
+s2_c2 <- mcmc(c2$s2)
+s2_c3 <- mcmc(c3$s2)
+s2_c4 <- mcmc(c4$s2)
+
+
+## 6. Results
+theta <- rbind(c1$theta, c2$theta, c3$theta, c4$theta)
+beta <- rbind(c1$beta, c2$beta, c3$beta, c4$beta)
+s2 <- c(c1$s2, c2$s2, c3$s2, c4$s2)
+accept <- c(c1$rate, c2$rate, c4$rate, c4$rate)
+# Potential scale reduction factors
+cat("\ntheta psrf:")
+cat("\n")
+gelman.diag(mcmc.list(theta_c1, theta_c2, theta_c3, theta_c4))$psrf[, 1L]
+cat(
+  "theta mpsrf:",
+  gelman.diag(mcmc.list(theta_c1, theta_c2, theta_c3, theta_c4))$mpsrf
+)
+cat(
+  "\nbeta psrf:",
+  gelman.diag(mcmc.list(beta_c1, beta_c2, beta_c3, beta_c4))$psrf[, 1L]
+)
+cat(
+  "\nbeta mpsrf:",
+  gelman.diag(mcmc.list(beta_c1, beta_c2, beta_c3, beta_c4))$mpsrf
+)
+cat(
+  "\ns2 psrf:",
+  gelman.diag(mcmc.list(s2_c1, s2_c2, s2_c3, s2_c4))$psrf[1L]
+)
+# Acceptance rate
+cat("\nAcceptance rate: ", mean(mean(c(c1$rate, c2$rate, c4$rate, c4$rate))))
+cat("\n\n")
+# Summary of beta and s2
+cat("< beta & s2 summary >\n")
+cat("beta1 mean:   ", mean(beta[, 1L]))
+cat("\nbeta1 95% ci: ", quantile(beta[, 1L], c(0.025, 0.975)))
+cat("\nbeta2 mean:   ", mean(beta[, 2L]))
+cat("\nbeta2 95% ci: ", quantile(beta[, 2L], c(0.025, 0.975)))
+cat("\ns2 mean:      ", mean(s2))
+cat("\ns2 95% ci:    ", quantile(s2, c(0.025, 0.975)))
+# Summary of theta
+theta_median <- apply(theta, MARGIN = 2L, FUN = median)
+theta_ci <- t(apply(theta, MARGIN = 2L, FUN = quantile, c(0.025, 0.975)))
+theta_al <- mean(theta_ci[, 2L] - theta_ci[, 1L])
+cat("\n\n< theta summary >\n")
+cat("theta 95% ci al: ", theta_al)
+cat("\ntheta aad:       ", aad(y, theta_median))
+cat("\ntheta aard:      ", aard(y, theta_median))
+cat("\ntheta asd:       ", asd(y, theta_median))
+cat("\ntheta asrd:      ", asrd(y, theta_median))
+# Original scale
+y_raw <- raw_income$mi_1989
+theta_os <- t(apply(theta,
+  MARGIN = 1L, FUN = function(x) x * sd(y_raw) + mean(y_raw)
+))
+theta_os_median <- apply(theta_os, MARGIN = 2L, FUN = median)
+theta_os_ci <- t(apply(theta_os, MARGIN = 2L, FUN = quantile, c(0.025, 0.975)))
+theta_os_al <- mean(theta_os_ci[, 2L] - theta_os_ci[, 1L])
+cat("\n\n< theta_os summary >\n")
+cat("theta_os 95% ci al: ", theta_os_al)
+cat("\ntheta_os aad:       ", aad(y_raw, theta_os_median))
+cat("\ntheta_os aard:      ", aard(y_raw, theta_os_median))
+cat("\ntheta_os asd:       ", asd(y_raw, theta_os_median))
+cat("\ntheta_os asrd:      ", asrd(y_raw, theta_os_median))
+cat("\n")