diff --git a/DESCRIPTION b/DESCRIPTION
index 5956330..f6460ca 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: retel
 Type: Package
 Title: Regularized Exponentially Tilted Empirical Likelihood
-Version: 0.0.0.9607
+Version: 0.0.0.9700
 Authors@R: c(
     person("Eunseop", "Kim", email = "markean@pm.me", 
             role = c("aut", "cre")),
diff --git a/retel-paper/code/figures/figure4b.R b/retel-paper/code/figures/figure4b.R
index f296c68..59e590a 100644
--- a/retel-paper/code/figures/figure4b.R
+++ b/retel-paper/code/figures/figure4b.R
@@ -1,139 +1,94 @@
-## 1. Load packages ----
+## 1. Load packages
+options(warn = -1)
+options(scipen = 999)
 library(coda)
 library(ggplot2)
 library(melt)
 library(mvtnorm)
 library(nloptr)
+library(retel)
 
 
-## 2. Functions ----
-d <- function(l, x, theta) {
-  mean(exp(l * (x - theta)))
-}
-# ETEL
-obj <- function(l, x, theta) {
-  d(l, x, theta)
-}
-gr_obj <- function(l, x, theta) {
-  mean(exp(l * (x - theta)) * (x - theta))
-}
-ETEL <- function(x, theta) {
-  if (theta > max(x) || theta < min(x)) {
-    return(-10000)
-  }
-  n <- length(x)
-  out <- nloptr(
-    x0 = 0, eval_f = obj, eval_grad_f = gr_obj, opts = opts, theta = theta,
-    x = x,
-  )
-  lambda <- out$solution
-  -n * log(d(lambda, x, theta)) + lambda * sum(x - theta)
-}
-# RETEL
-d2 <- function(x, l, theta, tau) {
-  n <- length(x)
-  sum(exp(l * (x - theta))) / (n + tau)
-}
-penalty2 <- function(l, tau, mu = 0, sigma = 1) {
-  (tau / (n + tau)) * exp(l * mu + 0.5 * l^2 * sigma^2)
+## 2. Constants
+# Sample size for each group
+n <- 2L
+# Standard deviation of prior
+sd_mu <- 10
+# MCMC parameters
+B <- 5000L
+# Tau
+tau <- 1
+# Optimization
+ctrl <- el_control(th = 1000, maxit_l = 50L, nthreads = 1L)
+
+
+## 3. Functions
+f <- function(x, par) {
+  x - par
 }
-continuous_obj2 <- function(l, x, theta, tau, mu = 0, sigma = 1) {
-  d2(x, l, theta, tau) + penalty2(l, tau, mu, sigma)
+# Prior and posterior computation
+log_prior <- function(mu, theta) {
+  dnorm(mu, mean = 0, sd = sd_mu, log = TRUE) -
+    sum(log(pi * (1 + (theta - mu)^2L)))
 }
-gr_continuous_obj2 <- function(l, x, theta, tau, mu = 0, sigma = 1) {
-  n <- length(x)
-  sum(exp(l * (x - theta)) * (x - theta)) / (n + tau) +
-    (tau / (n + tau)) * exp(l * mu + 0.5 * l^2 * sigma^2) * (mu + sigma^2 * l)
+loglr_retel_f <- function(theta, x1, x2, tau = 1) {
+  retel(f, x1, theta[1L], mean(x1) - theta[1L], 1, tau) +
+    retel(f, x2, theta[2L], mean(x2) - theta[2L], 1, tau)
 }
-RETEL1 <- function(x, theta, tau = 1, mu = 0, sigma = 1) {
-  n <- length(x)
-  out <- nloptr(
-    x0 = 0, eval_f = continuous_obj2,
-    eval_grad_f = gr_continuous_obj2, opts = opts, theta = theta, tau = tau,
-    x = x, mu = mu, sigma = sigma
-  )
-  lambda <- out$solution
-  -(n + 1) * log(continuous_obj2(lambda, x, theta, tau, mu, sigma)) +
-    lambda * sum(x - theta) + (lambda * mu + 0.5 * lambda^2 * sigma^2)
+log_posterior_unnormalized_retel_f <- function(mu_p, theta, x1, x2, tau = 1) {
+  log_prior(mu_p, theta) + loglr_retel_f(theta, x1, x2, tau)
 }
-RETEL2 <- function(x, theta, tau = 1, mu = 0, sigma = 1) {
-  n <- length(x)
-  out <- nloptr(
-    x0 = 0, eval_f = continuous_obj2,
-    eval_grad_f = gr_continuous_obj2, opts = opts, theta = theta, tau = tau,
-    x = x, mu = mu, sigma = sigma
-  )
-  l <- out$solution
-  -n * log(continuous_obj2(l, x, theta, tau, mu, sigma)) + l * sum(x - theta)
+loglr_retel_r <- function(theta, x1, x2, tau = 1) {
+  retel(f, x1, theta[1L], mean(x1) - theta[1L], 1, tau, "reduced") +
+    retel(f, x2, theta[2L], mean(x2) - theta[2L], 1, tau, "reduced")
 }
-# Prior and posterior computation
-log_prior <- function(mu, theta) {
-  dnorm(mu, mean = 0, sd = sd_mu, log = TRUE) -
-    sum(log(pi * (1 + (theta - mu)^2)))
+log_posterior_unnormalized_retel_r <- function(mu_p, theta, x1, x2, tau = 1) {
+  log_prior(mu_p, theta) + loglr_retel_r(theta, x1, x2, tau)
 }
-logLR_EL <- function(theta, x1, x2) {
+loglr_el <- function(theta, x1, x2) {
   fit1 <- el_mean(x1, par = theta[1L], control = ctrl)
   fit2 <- el_mean(x2, par = theta[2L], control = ctrl)
   lr1 <- if (conv(fit1)) logLR(fit1) else -10000
   lr2 <- if (conv(fit2)) logLR(fit2) else -10000
   lr1 + lr2
 }
-log_posterior_unnormalized_EL <- function(mu_p, theta, x1, x2) {
-  log_prior(mu_p, theta) + logLR_EL(theta, x1, x2)
-}
-
-logLR_ETEL <- function(theta, x1, x2) {
-  ETEL(x1, theta[1L]) + ETEL(x2, theta[2L])
-}
-log_posterior_unnormalized_ETEL <- function(mu_p, theta, x1, x2) {
-  log_prior(mu_p, theta) + logLR_ETEL(theta, x1, x2)
+log_posterior_unnormalized_el <- function(mu_p, theta, x1, x2) {
+  log_prior(mu_p, theta) + loglr_el(theta, x1, x2)
 }
-logLR_RETEL1 <- function(theta, x1, x2, tau = 1) {
-  RETEL1(x1, theta[1L], tau, mu = mean(x1) - theta[1L], sigma = 1) +
-    RETEL1(x2, theta[2L], tau, mu = mean(x2) - theta[2L], sigma = 1)
-}
-log_posterior_unnormalized_RETEL1 <- function(mu_p, theta, x1, x2, tau = 1) {
-  log_prior(mu_p, theta) + logLR_RETEL1(theta, x1, x2, tau)
-}
-logLR_RETEL2 <- function(theta, x1, x2, tau = 1) {
-  RETEL2(x1, theta[1L], tau, mu = mean(x1) - theta[1L], sigma = 1) +
-    RETEL2(x2, theta[2L], tau, mu = mean(x2) - theta[2L], sigma = 1)
+loglr_etel <- function(theta, x1, x2) {
+  if (theta[1L] > max(x1) || theta[1L] < min(x1)) {
+    loglr1 <- -10000
+  } else {
+    loglr1 <- etel(f, x1, theta[1L])
+  }
+  if (theta[2L] > max(x2) || theta[2L] < min(x2)) {
+    loglr2 <- -10000
+  } else {
+    loglr2 <- etel(f, x2, theta[2L])
+  }
+  loglr1 + loglr2
 }
-log_posterior_unnormalized_RETEL2 <- function(mu_p, theta, x1, x2, tau = 1) {
-  log_prior(mu_p, theta) + logLR_RETEL2(theta, x1, x2, tau)
+log_posterior_unnormalized_etel <- function(mu_p, theta, x1, x2) {
+  log_prior(mu_p, theta) + loglr_etel(theta, x1, x2)
 }
 
 
-## 3. Parameters ----
-# Sample size for each group
-n <- 2
-# Standard deviation of prior
-sd_mu <- 10
-# MCMC parameters
-B <- 5000
-# Tau
-tau <- 1
-# Optimization
-opts <- list("algorithm" = "NLOPT_LD_LBFGS", "xtol_rel" = 1e-04)
-ctrl <- el_control(th = 1000, maxit_l = 50, nthreads = 1L)
-
-
-## 4. Simulation ----
+## 4. Simulations
 set.seed(25424)
 # Generate data
 x1 <- rnorm(n, mean = -3, sd = 1)
 x2 <- rnorm(n, mean = 3, sd = 1)
-# RETEL_f ----
+# RETEL_f
 sd_proposal_mu <- 10
 sd_proposal_theta <- 2
 # MCMC chain 1
 mu_sample1 <- vector("numeric", length = B)
-mu_sample1[1L] <- rnorm(1)
-theta_sample1 <- matrix(nrow = B, ncol = 2)
+mu_sample1[1L] <- rnorm(1L)
+theta_sample1 <- matrix(nrow = B, ncol = 2L)
 theta_sample1[1L, ] <- c(rnorm(1L, mean = -2), rnorm(1L, mean = 2))
 acceptace1 <- vector("logical", length = B)
 acceptace1[1L] <- FALSE
-for (i in seq_len(B)[-1]) {
+for (i in seq_len(B)[-1L]) {
   # Sample proposal value
   prop <- rmvnorm(1L,
     mean = c(mu_sample1[i - 1L], theta_sample1[i - 1L, ]),
@@ -142,12 +97,12 @@ for (i in seq_len(B)[-1]) {
   mu_p <- prop[1L]
   theta_p <- prop[2L:3L]
   # Compute log ratio of unnormailzed posterior densities
-  logr <- log_posterior_unnormalized_RETEL1(mu_p, theta_p, x1, x2, tau) -
-    log_posterior_unnormalized_RETEL1(
+  logr <- log_posterior_unnormalized_retel_f(mu_p, theta_p, x1, x2, tau) -
+    log_posterior_unnormalized_retel_f(
       mu_sample1[i - 1L], theta_sample1[i - 1L, ], x1, x2, tau
     )
   # Sample uniform random variable
-  u <- runif(1)
+  u <- runif(1L)
   # Accept or reject
   if (isTRUE(log(u) < logr)) {
     mu_sample1[i] <- mu_p
@@ -155,18 +110,18 @@ for (i in seq_len(B)[-1]) {
     acceptace1[i] <- TRUE
   } else {
     mu_sample1[i] <- mu_sample1[i - 1L]
-    theta_sample1[i, ] <- theta_sample1[i - 1, ]
+    theta_sample1[i, ] <- theta_sample1[i - 1L, ]
     acceptace1[i] <- FALSE
   }
 }
 # MCMC chain 2
 mu_sample2 <- vector("numeric", length = B)
-mu_sample2[1L] <- rnorm(1)
-theta_sample2 <- matrix(nrow = B, ncol = 2)
+mu_sample2[1L] <- rnorm(1L)
+theta_sample2 <- matrix(nrow = B, ncol = 2L)
 theta_sample2[1L, ] <- c(rnorm(1L, mean = -2), rnorm(1L, mean = 2))
 acceptace2 <- vector("logical", length = B)
 acceptace2[1L] <- FALSE
-for (i in seq_len(B)[-1]) {
+for (i in seq_len(B)[-1L]) {
   # Sample proposal value
   prop <- rmvnorm(1L,
     mean = c(mu_sample2[i - 1L], theta_sample2[i - 1L, ]),
@@ -175,12 +130,12 @@ for (i in seq_len(B)[-1]) {
   mu_p <- prop[1L]
   theta_p <- prop[2L:3L]
   # Compute log ratio of unnormailzed posterior densities
-  logr <- log_posterior_unnormalized_RETEL1(mu_p, theta_p, x1, x2, tau) -
-    log_posterior_unnormalized_RETEL1(
+  logr <- log_posterior_unnormalized_retel_f(mu_p, theta_p, x1, x2, tau) -
+    log_posterior_unnormalized_retel_f(
       mu_sample2[i - 1L], theta_sample2[i - 1L, ], x1, x2, tau
     )
   # Sample uniform random variable
-  u <- runif(1)
+  u <- runif(1L)
   # Accept or reject
   if (isTRUE(log(u) < logr)) {
     mu_sample2[i] <- mu_p
@@ -188,35 +143,35 @@ for (i in seq_len(B)[-1]) {
     acceptace2[i] <- TRUE
   } else {
     mu_sample2[i] <- mu_sample2[i - 1L]
-    theta_sample2[i, ] <- theta_sample2[i - 1, ]
+    theta_sample2[i, ] <- theta_sample2[i - 1L, ]
     acceptace2[i] <- FALSE
   }
 }
 # KL divergence
-df_retel1 <- data.frame(
+df_retel_f <- data.frame(
   index = seq_len(2L * B), mu = c(mu_sample1, mu_sample2),
   theta1 = c(theta_sample1[, 1L], theta_sample2[, 1L]),
   theta2 = c(theta_sample1[, 2L], theta_sample2[, 2L]),
   acceptace = c(acceptace1, acceptace2),
   chain = rep(c(1L, 2L), each = B),
-  type = "retel1"
+  type = "retel_f"
 )
 # Potential scale reduction factors
-mu_c1 <- mcmc(df_retel1$mu[seq_len(B)])
-mu_c2 <- mcmc(df_retel1$mu[seq(B + 1, 2 * B)])
-retel1_mu_psrf <- gelman.diag(mcmc.list(mu_c1, mu_c2))$psrf[1L]
+mu_c1 <- mcmc(df_retel_f$mu[seq_len(B)])
+mu_c2 <- mcmc(df_retel_f$mu[seq(B + 1L, 2L * B)])
+retel_f_mu_psrf <- gelman.diag(mcmc.list(mu_c1, mu_c2))$psrf[1L]
 
-# RETEL_r ----
+# RETEL_r
 sd_proposal_mu <- 10
 sd_proposal_theta <- 1
 # MCMC chain 1
 mu_sample1 <- vector("numeric", length = B)
-mu_sample1[1L] <- rnorm(1)
-theta_sample1 <- matrix(nrow = B, ncol = 2)
+mu_sample1[1L] <- rnorm(1L)
+theta_sample1 <- matrix(nrow = B, ncol = 2L)
 theta_sample1[1L, ] <- c(rnorm(1L, mean = -2), rnorm(1L, mean = 2))
 acceptace1 <- vector("logical", length = B)
 acceptace1[1L] <- FALSE
-for (i in seq_len(B)[-1]) {
+for (i in seq_len(B)[-1L]) {
   # Sample proposal value
   prop <- rmvnorm(1L,
     mean = c(mu_sample1[i - 1L], theta_sample1[i - 1L, ]),
@@ -225,12 +180,12 @@ for (i in seq_len(B)[-1]) {
   mu_p <- prop[1L]
   theta_p <- prop[2L:3L]
   # Compute log ratio of unnormailzed posterior densities
-  logr <- log_posterior_unnormalized_RETEL2(mu_p, theta_p, x1, x2, tau) -
-    log_posterior_unnormalized_RETEL2(
+  logr <- log_posterior_unnormalized_retel_r(mu_p, theta_p, x1, x2, tau) -
+    log_posterior_unnormalized_retel_r(
       mu_sample1[i - 1L], theta_sample1[i - 1L, ], x1, x2, tau
     )
   # Sample uniform random variable
-  u <- runif(1)
+  u <- runif(1L)
   # Accept or reject
   if (isTRUE(log(u) < logr)) {
     mu_sample1[i] <- mu_p
@@ -238,18 +193,18 @@ for (i in seq_len(B)[-1]) {
     acceptace1[i] <- TRUE
   } else {
     mu_sample1[i] <- mu_sample1[i - 1L]
-    theta_sample1[i, ] <- theta_sample1[i - 1, ]
+    theta_sample1[i, ] <- theta_sample1[i - 1L, ]
     acceptace1[i] <- FALSE
   }
 }
 # MCMC chain 2
 mu_sample2 <- vector("numeric", length = B)
-mu_sample2[1L] <- rnorm(1)
-theta_sample2 <- matrix(nrow = B, ncol = 2)
+mu_sample2[1L] <- rnorm(1L)
+theta_sample2 <- matrix(nrow = B, ncol = 2L)
 theta_sample2[1L, ] <- c(rnorm(1L, mean = -2), rnorm(1L, mean = 2))
 acceptace2 <- vector("logical", length = B)
 acceptace2[1L] <- FALSE
-for (i in seq_len(B)[-1]) {
+for (i in seq_len(B)[-1L]) {
   # Sample proposal value
   prop <- rmvnorm(1L,
     mean = c(mu_sample2[i - 1L], theta_sample2[i - 1L, ]),
@@ -258,12 +213,12 @@ for (i in seq_len(B)[-1]) {
   mu_p <- prop[1L]
   theta_p <- prop[2L:3L]
   # Compute log ratio of unnormailzed posterior densities
-  logr <- log_posterior_unnormalized_RETEL2(mu_p, theta_p, x1, x2, tau) -
-    log_posterior_unnormalized_RETEL2(
+  logr <- log_posterior_unnormalized_retel_r(mu_p, theta_p, x1, x2, tau) -
+    log_posterior_unnormalized_retel_r(
       mu_sample2[i - 1L], theta_sample2[i - 1L, ], x1, x2, tau
     )
   # Sample uniform random variable
-  u <- runif(1)
+  u <- runif(1L)
   # Accept or reject
   if (isTRUE(log(u) < logr)) {
     mu_sample2[i] <- mu_p
@@ -271,38 +226,38 @@ for (i in seq_len(B)[-1]) {
     acceptace2[i] <- TRUE
   } else {
     mu_sample2[i] <- mu_sample2[i - 1L]
-    theta_sample2[i, ] <- theta_sample2[i - 1, ]
+    theta_sample2[i, ] <- theta_sample2[i - 1L, ]
     acceptace2[i] <- FALSE
   }
 }
 # KL divergence
-df_retel2 <- data.frame(
+df_retel_r <- data.frame(
   index = seq_len(2L * B), mu = c(mu_sample1, mu_sample2),
   theta1 = c(theta_sample1[, 1L], theta_sample2[, 1L]),
   theta2 = c(theta_sample1[, 2L], theta_sample2[, 2L]),
   acceptace = c(acceptace1, acceptace2),
   chain = rep(c(1L, 2L), each = B),
-  type = "retel2"
+  type = "retel_r"
 )
 # Potential scale reduction factors
-mu_c1 <- mcmc(df_retel2$mu[seq_len(B)])
-mu_c2 <- mcmc(df_retel2$mu[seq(B + 1, 2 * B)])
-retel2_mu_psrf <- gelman.diag(mcmc.list(mu_c1, mu_c2))$psrf[1L]
+mu_c1 <- mcmc(df_retel_r$mu[seq_len(B)])
+mu_c2 <- mcmc(df_retel_r$mu[seq(B + 1L, 2L * B)])
+retel_r_mu_psrf <- gelman.diag(mcmc.list(mu_c1, mu_c2))$psrf[1L]
 
-# EL ----
+# EL
 sd_proposal_mu <- 10
 sd_proposal_theta <- 0.05
 # MCMC chain 1
 mu_sample1 <- vector("numeric", length = B)
-mu_sample1[1L] <- rnorm(1)
-theta_sample1 <- matrix(nrow = B, ncol = 2)
+mu_sample1[1L] <- rnorm(1L)
+theta_sample1 <- matrix(nrow = B, ncol = 2L)
 theta_sample1[1L, ] <- c(
   runif(1L, min = min(x1), max = max(x1)),
   runif(1L, min = min(x2), max = max(x2))
 )
 acceptace1 <- vector("logical", length = B)
 acceptace1[1L] <- FALSE
-for (i in seq_len(B)[-1]) {
+for (i in seq_len(B)[-1L]) {
   # Sample proposal value
   prop <- rmvnorm(1L,
     mean = c(mu_sample1[i - 1L], theta_sample1[i - 1L, ]),
@@ -311,12 +266,12 @@ for (i in seq_len(B)[-1]) {
   mu_p <- prop[1L]
   theta_p <- prop[2L:3L]
   # Compute log ratio of unnormailzed posterior densities
-  logr <- log_posterior_unnormalized_EL(mu_p, theta_p, x1, x2) -
-    log_posterior_unnormalized_EL(
+  logr <- log_posterior_unnormalized_el(mu_p, theta_p, x1, x2) -
+    log_posterior_unnormalized_el(
       mu_sample1[i - 1L], theta_sample1[i - 1L, ], x1, x2
     )
   # Sample uniform random variable
-  u <- runif(1)
+  u <- runif(1L)
   # Accept or reject
   if (isTRUE(log(u) < logr)) {
     mu_sample1[i] <- mu_p
@@ -324,21 +279,21 @@ for (i in seq_len(B)[-1]) {
     acceptace1[i] <- TRUE
   } else {
     mu_sample1[i] <- mu_sample1[i - 1L]
-    theta_sample1[i, ] <- theta_sample1[i - 1, ]
+    theta_sample1[i, ] <- theta_sample1[i - 1L, ]
     acceptace1[i] <- FALSE
   }
 }
 # MCMC chain 2
 mu_sample2 <- vector("numeric", length = B)
-mu_sample2[1L] <- rnorm(1)
-theta_sample2 <- matrix(nrow = B, ncol = 2)
+mu_sample2[1L] <- rnorm(1L)
+theta_sample2 <- matrix(nrow = B, ncol = 2L)
 theta_sample2[1L, ] <- c(
   runif(1L, min = min(x1), max = max(x1)),
   runif(1L, min = min(x2), max = max(x2))
 )
 acceptace2 <- vector("logical", length = B)
 acceptace2[1L] <- FALSE
-for (i in seq_len(B)[-1]) {
+for (i in seq_len(B)[-1L]) {
   # Sample proposal value
   prop <- rmvnorm(1L,
     mean = c(mu_sample2[i - 1L], theta_sample2[i - 1L, ]),
@@ -347,12 +302,12 @@ for (i in seq_len(B)[-1]) {
   mu_p <- prop[1L]
   theta_p <- prop[2L:3L]
   # Compute log ratio of unnormailzed posterior densities
-  logr <- log_posterior_unnormalized_EL(mu_p, theta_p, x1, x2) -
-    log_posterior_unnormalized_EL(
+  logr <- log_posterior_unnormalized_el(mu_p, theta_p, x1, x2) -
+    log_posterior_unnormalized_el(
       mu_sample2[i - 1L], theta_sample2[i - 1L, ], x1, x2
     )
   # Sample uniform random variable
-  u <- runif(1)
+  u <- runif(1L)
   # Accept or reject
   if (isTRUE(log(u) < logr)) {
     mu_sample2[i] <- mu_p
@@ -360,7 +315,7 @@ for (i in seq_len(B)[-1]) {
     acceptace2[i] <- TRUE
   } else {
     mu_sample2[i] <- mu_sample2[i - 1L]
-    theta_sample2[i, ] <- theta_sample2[i - 1, ]
+    theta_sample2[i, ] <- theta_sample2[i - 1L, ]
     acceptace2[i] <- FALSE
   }
 }
@@ -375,23 +330,23 @@ df_el <- data.frame(
 )
 # Potential scale reduction factors
 mu_c1 <- mcmc(df_el$mu[seq_len(B)])
-mu_c2 <- mcmc(df_el$mu[seq(B + 1, 2 * B)])
+mu_c2 <- mcmc(df_el$mu[seq(B + 1L, 2L * B)])
 el_mu_psrf <- gelman.diag(mcmc.list(mu_c1, mu_c2))$psrf[1L]
 
-# ETEL ----
+# ETEL
 sd_proposal_mu <- 6
 sd_proposal_theta <- 0.2
 # MCMC chain 1
 mu_sample1 <- vector("numeric", length = B)
-mu_sample1[1L] <- runif(1, min = min(x1), max = max(x1))
-theta_sample1 <- matrix(nrow = B, ncol = 2)
+mu_sample1[1L] <- runif(1L, min = min(x1), max = max(x1))
+theta_sample1 <- matrix(nrow = B, ncol = 2L)
 theta_sample1[1L, ] <- c(
   runif(1L, min = min(x1), max = max(x1)),
   runif(1L, min = min(x2), max = max(x2))
 )
 acceptace1 <- vector("logical", length = B)
 acceptace1[1L] <- FALSE
-for (i in seq_len(B)[-1]) {
+for (i in seq_len(B)[-1L]) {
   # Sample proposal value
   prop <- rmvnorm(1L,
     mean = c(mu_sample1[i - 1L], theta_sample1[i - 1L, ]),
@@ -400,12 +355,12 @@ for (i in seq_len(B)[-1]) {
   mu_p <- prop[1L]
   theta_p <- prop[2L:3L]
   # Compute log ratio of unnormailzed posterior densities
-  logr <- log_posterior_unnormalized_ETEL(mu_p, theta_p, x1, x2) -
-    log_posterior_unnormalized_ETEL(
+  logr <- log_posterior_unnormalized_etel(mu_p, theta_p, x1, x2) -
+    log_posterior_unnormalized_etel(
       mu_sample1[i - 1L], theta_sample1[i - 1L, ], x1, x2
     )
   # Sample uniform random variable
-  u <- runif(1)
+  u <- runif(1L)
   # Accept or reject
   if (isTRUE(log(u) < logr)) {
     mu_sample1[i] <- mu_p
@@ -413,21 +368,21 @@ for (i in seq_len(B)[-1]) {
     acceptace1[i] <- TRUE
   } else {
     mu_sample1[i] <- mu_sample1[i - 1L]
-    theta_sample1[i, ] <- theta_sample1[i - 1, ]
+    theta_sample1[i, ] <- theta_sample1[i - 1L, ]
     acceptace1[i] <- FALSE
   }
 }
 # MCMC chain 2
 mu_sample2 <- vector("numeric", length = B)
-mu_sample2[1L] <- runif(1, min = min(x2), max = max(x2))
-theta_sample2 <- matrix(nrow = B, ncol = 2)
+mu_sample2[1L] <- runif(1L, min = min(x2), max = max(x2))
+theta_sample2 <- matrix(nrow = B, ncol = 2L)
 theta_sample2[1L, ] <- c(
   runif(1L, min = min(x1), max = max(x1)),
   runif(1L, min = min(x2), max = max(x2))
 )
 acceptace2 <- vector("logical", length = B)
 acceptace2[1L] <- FALSE
-for (i in seq_len(B)[-1]) {
+for (i in seq_len(B)[-1L]) {
   # Sample proposal value
   prop <- rmvnorm(1L,
     mean = c(mu_sample2[i - 1L], theta_sample2[i - 1L, ]),
@@ -436,12 +391,12 @@ for (i in seq_len(B)[-1]) {
   mu_p <- prop[1L]
   theta_p <- prop[2L:3L]
   # Compute log ratio of unnormailzed posterior densities
-  logr <- log_posterior_unnormalized_ETEL(mu_p, theta_p, x1, x2) -
-    log_posterior_unnormalized_ETEL(
+  logr <- log_posterior_unnormalized_etel(mu_p, theta_p, x1, x2) -
+    log_posterior_unnormalized_etel(
       mu_sample2[i - 1L], theta_sample2[i - 1L, ], x1, x2
     )
   # Sample uniform random variable
-  u <- runif(1)
+  u <- runif(1L)
   # Accept or reject
   if (isTRUE(log(u) < logr)) {
     mu_sample2[i] <- mu_p
@@ -449,7 +404,7 @@ for (i in seq_len(B)[-1]) {
     acceptace2[i] <- TRUE
   } else {
     mu_sample2[i] <- mu_sample2[i - 1L]
-    theta_sample2[i, ] <- theta_sample2[i - 1, ]
+    theta_sample2[i, ] <- theta_sample2[i - 1L, ]
     acceptace2[i] <- FALSE
   }
 }
@@ -464,22 +419,23 @@ df_etel <- data.frame(
 )
 # Potential scale reduction factors
 mu_c1 <- mcmc(df_etel$mu[seq_len(B)])
-mu_c2 <- mcmc(df_etel$mu[seq(B + 1, 2 * B)])
+mu_c2 <- mcmc(df_etel$mu[seq(B + 1L, 2L * B)])
 etel_mu_psrf <- gelman.diag(mcmc.list(mu_c1, mu_c2))$psrf[1L]
 
 
-# Plot ----
+# Plot
+# 4 x 3
 # Potential scale reduction factor
-c(retel1_mu_psrf, retel2_mu_psrf, el_mu_psrf, etel_mu_psrf)
+c(retel_f_mu_psrf, retel_r_mu_psrf, el_mu_psrf, etel_mu_psrf)
 # Acceptance rate
 c(
-  mean(df_retel1$acceptace), mean(df_retel2$acceptace), mean(df_el$acceptace),
+  mean(df_retel_f$acceptace), mean(df_retel_r$acceptace), mean(df_el$acceptace),
   mean(df_etel$acceptace)
 )
 legend_labels <- c(
   expression(RETEL[italic(f)]), expression(RETEL[italic(r)]), "EL", "ETEL"
 )
-df <- rbind(df_retel1, df_retel2, df_el, df_etel)
+df <- rbind(df_retel_f, df_retel_r, df_el, df_etel)
 ggplot(df, aes(x = mu, color = type, linetype = type)) +
   geom_density(adjust = 2) +
   geom_vline(
@@ -508,17 +464,17 @@ ggplot(df, aes(x = mu, color = type, linetype = type)) +
   ) +
   scale_colour_manual(
     labels = legend_labels,
-    breaks = c("retel1", "retel2", "el", "etel"),
-    values = c("black", "red", "blue", "green")
+    values = c("black", "red", "blue", "green"),
+    breaks = c("retel_f", "retel_r", "el", "etel")
   ) +
   scale_linetype_manual(
     labels = legend_labels,
-    breaks = c("retel1", "retel2", "el", "etel"),
-    values = c("solid", "dashed", "dotdash", "twodash")
+    values = c("solid", "dashed", "dotdash", "twodash"),
+    breaks = c("retel_f", "retel_r", "el", "etel")
   ) +
   scale_shape_manual(
     labels = legend_labels,
-    breaks = c("retel1", "retel2", "el", "etel"),
-    values = c(15, 16, 2, 5)
+    values = c(15, 16, 2, 5),
+    breaks = c("retel_f", "retel_r", "el", "etel")
   ) +
   xlim(c(-6, 6))
diff --git a/retel-paper/figures/figure4b.pdf b/retel-paper/figures/figure4b.pdf
index edff70a..4951b20 100644
Binary files a/retel-paper/figures/figure4b.pdf and b/retel-paper/figures/figure4b.pdf differ