From 6efea35873a52f296919f2425115693a5a786c23 Mon Sep 17 00:00:00 2001 From: markean Date: Sat, 20 Jan 2024 20:42:30 -0500 Subject: [PATCH] application retel_f added --- DESCRIPTION | 2 +- retel-paper/application/retel_f.txt | 38 +++++ retel-paper/code/application/retel_f.R | 197 +++++++++++++++++++++++++ retel-paper/renv.lock | 12 +- 4 files changed, 242 insertions(+), 7 deletions(-) create mode 100644 retel-paper/application/retel_f.txt create mode 100644 retel-paper/code/application/retel_f.R diff --git a/DESCRIPTION b/DESCRIPTION index 492392e..d93d4a9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: retel Type: Package Title: Regularized Exponentially Tilted Empirical Likelihood -Version: 0.0.0.9703 +Version: 0.0.0.9704 Authors@R: c( person("Eunseop", "Kim", email = "markean@pm.me", role = c("aut", "cre")), diff --git a/retel-paper/application/retel_f.txt b/retel-paper/application/retel_f.txt new file mode 100644 index 0000000..139ce79 --- /dev/null +++ b/retel-paper/application/retel_f.txt @@ -0,0 +1,38 @@ +< Method = RETEL_f> +< MCMC runs = 1000000 > + +theta psrf: + [1] 1.004675 1.005043 1.005131 1.004328 1.001039 1.001766 1.001821 1.003741 + [9] 1.002047 1.003690 1.002943 1.002106 1.001036 1.002151 1.004675 1.000847 +[17] 1.003991 1.003880 1.004994 1.001463 1.002753 1.002078 1.002623 1.002519 +[25] 1.004075 1.007076 1.006724 1.002586 1.002243 1.001770 1.001479 1.001820 +[33] 1.002741 1.004970 1.005725 1.001699 1.001418 1.003418 1.003746 1.001834 +[41] 1.001401 1.001166 1.001079 1.000889 1.000756 1.002257 1.003136 1.005162 +[49] 1.000747 1.001715 1.001792 +theta mpsrf: 1.0416 +beta psrf: 1.000159 1.000238 +beta mpsrf: 1.000156 +s2 psrf: 1.008876 +Acceptance rate: 0.130911 + +< beta & s2 summary > +beta1 mean: 0.01503831 +beta1 95% ci: -0.1047506 0.1350408 +beta2 mean: 0.9332734 +beta2 95% ci: 0.8127158 1.053873 +s2 mean: 0.9406872 +s2 95% ci: 0.4964409 1.608398 + +< theta summary > +theta 95% ci al: 3.7542 +theta aad: 0.2781389 +theta aard: 0.9111071 +theta asd: 0.1108935 +theta asrd: 3.982455 + +< theta_os summary > +theta_os 95% ci al: 20728.25 +theta_os aad: 1535.702 +theta_os aard: 0.03957153 +theta_os asd: 3380624 +theta_os asrd: 0.002336617 diff --git a/retel-paper/code/application/retel_f.R b/retel-paper/code/application/retel_f.R new file mode 100644 index 0000000..60225f7 --- /dev/null +++ b/retel-paper/code/application/retel_f.R @@ -0,0 +1,197 @@ +## 1. Load packages +options(warn = -1) +options(scipen = 999) +suppressMessages(library(actuar)) +library(coda) +suppressMessages(library(here)) +suppressMessages(i_am("code/application/retel_f.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) { + g <- f(y, theta) + log_pd(theta, b, s2) + + retel(f, y, theta, colMeans(g), var(g), log(m)) +} +# 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 + 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] + ) + # 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_f>\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") diff --git a/retel-paper/renv.lock b/retel-paper/renv.lock index fda060d..63da10e 100644 --- a/retel-paper/renv.lock +++ b/retel-paper/renv.lock @@ -1161,14 +1161,14 @@ }, "ps": { "Package": "ps", - "Version": "1.7.5", + "Version": "1.7.6", "Source": "Repository", - "Repository": "CRAN", + "Repository": "RSPM", "Requirements": [ "R", "utils" ], - "Hash": "709d852d33178db54b17c722e5b1e594" + "Hash": "dd2b9319ee0656c8acf45c7f40c59de7" }, "purrr": { "Package": "purrr", @@ -1235,13 +1235,13 @@ }, "retel": { "Package": "retel", - "Version": "0.0.0.9401", + "Version": "0.0.0.9703", "Source": "GitHub", "RemoteType": "github", "RemoteUsername": "markean", "RemoteRepo": "retel", "RemoteRef": "HEAD", - "RemoteSha": "a9ca276f368e9e4df3cb5ec8e612b7242227109f", + "RemoteSha": "c6c0586dbc6eb8908725d5ca7bdd6b4695ade631", "RemoteHost": "api.github.com", "Requirements": [ "Matrix", @@ -1250,7 +1250,7 @@ "matrixcalc", "nloptr" ], - "Hash": "aa14eac22f1e11105fce40f1047150b3" + "Hash": "673b1135be8b6847071373f0c71e9848" }, "rlang": { "Package": "rlang",