From 8aa1c701e9b5b4335d39ea791b2fc51d720fdf38 Mon Sep 17 00:00:00 2001 From: markean Date: Sat, 20 Jan 2024 13:11:15 -0500 Subject: [PATCH] application etel added --- DESCRIPTION | 2 +- retel-paper/application/etel.txt | 38 ++++++ retel-paper/code/application/etel.R | 195 ++++++++++++++++++++++++++++ retel-paper/renv.lock | 23 ++++ 4 files changed, 257 insertions(+), 1 deletion(-) create mode 100644 retel-paper/application/etel.txt create mode 100644 retel-paper/code/application/etel.R diff --git a/DESCRIPTION b/DESCRIPTION index b6dd852..66d6fb1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: retel Type: Package Title: Regularized Exponentially Tilted Empirical Likelihood -Version: 0.0.0.9701 +Version: 0.0.0.9702 Authors@R: c( person("Eunseop", "Kim", email = "markean@pm.me", role = c("aut", "cre")), diff --git a/retel-paper/application/etel.txt b/retel-paper/application/etel.txt new file mode 100644 index 0000000..b8bcb83 --- /dev/null +++ b/retel-paper/application/etel.txt @@ -0,0 +1,38 @@ +< Method = ETEL> +< MCMC runs = 1000000 > + +theta psrf: + [1] 1.000572 1.001384 1.000966 1.006470 1.003515 1.000916 1.002042 1.000081 + [9] 1.004781 1.003241 1.002479 1.005392 1.000518 1.000366 1.003687 1.002037 +[17] 1.004757 1.004614 1.002384 1.001342 1.003717 1.001176 1.000437 1.002256 +[25] 1.000432 1.007704 1.004953 1.001959 1.001242 1.001574 1.001423 1.003690 +[33] 1.001718 1.005683 1.000767 1.000986 1.000950 1.001274 1.001884 1.006432 +[41] 1.000426 1.001921 1.001409 1.007172 1.004278 1.003832 1.003266 1.000522 +[49] 1.000887 1.008399 1.001018 +theta mpsrf: 1.037999 +beta psrf: 1.00013 1.000173 +beta mpsrf: 1.000106 +s2 psrf: 1.005759 +Acceptance rate: 0.127195 + +< beta & s2 summary > +beta1 mean: 0.01515873 +beta1 95% ci: -0.1031014 0.1335065 +beta2 mean: 0.9325279 +beta2 95% ci: 0.8142174 1.051342 +s2 mean: 0.928918 +s2 95% ci: 0.4898741 1.573498 + +< theta summary > +theta 95% ci al: 3.728802 +theta aad: 0.2792567 +theta aard: 0.9421465 +theta asd: 0.1149375 +theta asrd: 4.696289 + +< theta_os summary > +theta_os 95% ci al: 20588.02 +theta_os aad: 1541.874 +theta_os aard: 0.03974086 +theta_os asd: 3503906 +theta_os asrd: 0.002438499 diff --git a/retel-paper/code/application/etel.R b/retel-paper/code/application/etel.R new file mode 100644 index 0000000..7454675 --- /dev/null +++ b/retel-paper/code/application/etel.R @@ -0,0 +1,195 @@ +## 1. Load packages +options(warn = -1) +options(scipen = 999) +suppressMessages(library(actuar)) +library(coda) +suppressMessages(library(here)) +suppressMessages(i_am("code/application/etel.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) { + log_pd(theta, b, s2) + etel(f, y, theta) +} +# 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 = ETEL>\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 87a19ba..fda060d 100644 --- a/retel-paper/renv.lock +++ b/retel-paper/renv.lock @@ -131,6 +131,19 @@ ], "Hash": "4f57884290cc75ab22f4af9e9d4ca862" }, + "actuar": { + "Package": "actuar", + "Version": "3.3-4", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "R", + "expint", + "graphics", + "stats" + ], + "Hash": "31c925f608675e450657c4622bd68b68" + }, "backports": { "Package": "backports", "Version": "1.4.1", @@ -503,6 +516,16 @@ ], "Hash": "daf4a1246be12c1fa8c7705a0935c1a0" }, + "expint": { + "Package": "expint", + "Version": "0.1-8", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "R" + ], + "Hash": "57208be351e926869d7121aca0c80a9b" + }, "fansi": { "Package": "fansi", "Version": "1.0.6",