+
+Same model but now assuming observations are from Gamma distribution:
+
+``` r
+
+model2 <- bsm_ng(airquality$Ozone,
+ xreg = xreg,
+ beta = normal(rep(0, ncol(xreg)), 0, 1),
+ distribution = "gamma",
+ phi = gamma_prior(1, 2, 0.01),
+ sd_level = gamma_prior(1, 2, 0.1),
+ sd_slope = gamma_prior(1, 2, 0.1))
+
+fit2 <- run_mcmc(model2, iter = 20000, burnin = 5000, particles = 10)
+fit2
+#>
+#> Call:
+#> run_mcmc.nongaussian(model = model2, iter = 20000, particles = 10,
+#> burnin = 5000)
+#>
+#> Iterations = 5001:20000
+#> Thinning interval = 1
+#> Length of the final jump chain = 3858
+#>
+#> Acceptance rate after the burn-in period: 0.257
+#>
+#> Summary for theta:
+#>
+#> variable Mean SE SD 2.5% 97.5% ESS
+#> Temp 0.052808820 0.0002404538 0.008701489 0.0353736458 0.06992423 1310
+#> Wind -0.057351094 0.0004338213 0.015411504 -0.0873384757 -0.02700112 1262
+#> phi 4.006977632 0.0159088062 0.536273508 3.0263444882 5.15527365 1136
+#> sd_level 0.057158663 0.0020138200 0.035366227 0.0083794202 0.14651419 308
+#> sd_slope 0.003894013 0.0001752319 0.003654978 0.0004250207 0.01374575 435
+#> SE_IS ESS_IS
+#> 7.128104e-05 14485
+#> 1.263047e-04 13905
+#> 4.411840e-03 14611
+#> 2.927386e-04 10591
+#> 3.031489e-05 7766
+#>
+#> Summary for alpha_154:
+#>
+#> variable time Mean SE SD 2.5% 97.5% ESS
+#> level 154 -0.200656509 0.0201721601 0.73134471 -1.62501396 1.24522802 1314
+#> slope 154 -0.002689176 0.0005121944 0.02289051 -0.04650504 0.04724173 1997
+#> SE_IS ESS_IS
+#> 0.005987284 9458
+#> 0.000191620 6448
+#>
+#> Run time:
+#> user system elapsed
+#> 7.50 0.01 7.71
+```
+
+Comparison:
+
+``` r
+pred2 <- fitted(fit2, model2)
+
+bind_rows(list(Gaussian = pred, Gamma = pred2), .id = "Model") |>
+ ggplot(aes(x = Time, y = Mean)) +
+ geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`, fill = Model),
+ alpha = 0.25) +
+ geom_line(aes(colour = Model)) +
+ geom_point(data = obs,
+ aes(x = Time, y = Ozone)) +
+ theme_bw()
+```
+
+
+
+Now let’s assume that we also want to use the solar radiation variable
+as predictor for ozone. As it contains few missing values, we cannot use
+it directly. As the number of missing time points is very small, simple
+imputation would likely be acceptable, but let’s consider more another
+approach. For simplicity, the slope terms of the previous models are now
+omitted, and we focus on the Gaussian case. Let $\mu_t$ be the true
+solar radiation at time $t$. Now for ozone $O_t$ we assume following
+model:
+
+$O_t = D_t + \alpha_t + \beta_S \mu_t + \sigma_\epsilon \epsilon_t$
+$\alpha_{t+1} = \alpha_t + \sigma_\eta\eta_t$
+$\alpha_1 \sim N(0, 100^2\textrm{I})$,
+wheere $D_t = \beta X_t$ contains regression terms related to wind and
+temperature, $\alpha_t$ is the time varying intercept term, and
+$\beta_S$ is the effect of solar radiation $\mu_t$.
+
+Now for the observed solar radiation $S_t$ we assume
+
+$S_t = \mu_t$
+$\mu_{t+1} = \mu_t + \sigma_\xi\xi_t,$
+$\mu_1 \sim N(0, 100^2)$,
+i.e. we assume as simple random walk for the $\mu$ which we observe
+without error or not at all (there is no error term in the observation
+equation $S_t=\mu_t$).
+
+We combine these two models as a bivariate Gaussian model with
+`ssm_mlg`:
+
+``` r
+# predictors (not including solar radiation) for ozone
+xreg <- airquality |> select(Wind, Temp) |> as.matrix()
+
+# Function which outputs new model components given the parameter vector theta
+update_fn <- function(theta) {
+ D <- rbind(t(xreg %*% theta[1:2]), 1)
+ Z <- matrix(c(1, 0, theta[3], 1), 2, 2)
+ R <- diag(exp(theta[4:5]))
+ H <- diag(c(exp(theta[6]), 0))
+ # add third dimension so we have p x n x 1, p x m x 1, p x p x 1 arrays
+ dim(Z)[3] <- dim(R)[3] <- dim(H)[3] <- 1
+ list(D = D, Z = Z, R = R, H = H)
+}
+
+# Function for log-prior density
+prior_fn <- function(theta) {
+ sum(dnorm(theta[1:3], 0, 10, log = TRUE)) +
+ sum(dgamma(exp(theta[4:6]), 2, 0.01, log = TRUE)) +
+ sum(theta[4:6]) # log-jacobian
+}
+
+init_theta <- c(0, 0, 0, log(50), log(5), log(20))
+comps <- update_fn(init_theta)
+
+model <- ssm_mlg(y = cbind(Ozone = airquality$Ozone, Solar = airquality$Solar.R),
+ Z = comps$Z, D = comps$D, H = comps$H, T = diag(2), R = comps$R,
+ a1 = rep(0, 2), P1 = diag(100, 2), init_theta = init_theta,
+ state_names = c("alpha", "mu"), update_fn = update_fn,
+ prior_fn = prior_fn)
+
+fit <- run_mcmc(model, iter = 60000, burnin = 10000)
+fit
+#>
+#> Call:
+#> run_mcmc.lineargaussian(model = model, iter = 60000, burnin = 10000)
+#>
+#> Iterations = 10001:60000
+#> Thinning interval = 1
+#> Length of the final jump chain = 12234
+#>
+#> Acceptance rate after the burn-in period: 0.245
+#>
+#> Summary for theta:
+#>
+#> variable Mean SE SD 2.5% 97.5% ESS
+#> theta_1 -3.89121114 0.0233827004 0.58715113 -5.0085134 -2.6915137 631
+#> theta_2 0.98712126 0.0051506907 0.18819758 0.5917823 1.3475147 1335
+#> theta_3 0.06324657 0.0004672314 0.02417334 0.0141425 0.1100184 2677
+#> theta_4 0.82577262 0.0165661049 0.67134723 -0.6840637 1.9160168 1642
+#> theta_5 4.75567621 0.0010887250 0.05858454 4.6446809 4.8704036 2895
+#> theta_6 3.05462451 0.0014803971 0.07640392 2.9032635 3.2028023 2664
+#>
+#> Summary for alpha_154:
+#>
+#> variable time Mean SE SD 2.5% 97.5% ESS
+#> alpha 154 -16.44435 0.3659912 14.99708 -46.321645 13.01863 1679
+#> mu 154 223.60490 1.3409568 116.49063 -6.206301 453.18554 7546
+#>
+#> Run time:
+#> user system elapsed
+#> 7.41 0.11 7.83
+```
+
+Draw predictions:
+
+``` r
+pred <- fitted(fit, model)
+
+obs <- data.frame(Time = 1:nrow(airquality),
+ Solar = airquality$Solar.R) |> filter(!is.na(Solar))
+
+pred |> filter(Variable == "Solar") |>
+ ggplot(aes(x = Time, y = Mean)) +
+ geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`),
+ alpha = 0.5, fill = "steelblue") +
+ geom_line() +
+ geom_point(data = obs,
+ aes(x = Time, y = Solar), colour = "Tomato") +
+ theme_bw()
+```
+
+
+
+``` r
+
+
+obs <- data.frame(Time = 1:nrow(airquality),
+ Ozone = airquality$Ozone) |> filter(!is.na(Ozone))
+
+pred |> filter(Variable == "Ozone") |>
+ ggplot(aes(x = Time, y = Mean)) +
+ geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`),
+ alpha = 0.5, fill = "steelblue") +
+ geom_line() +
+ geom_point(data = obs,
+ aes(x = Time, y = Ozone), colour = "Tomato") +
+ theme_bw()
+```
+
+
+
+See more examples in the paper, vignettes, and in the docs.
diff --git a/benchmarks/replications.Rmd b/benchmarks/replications.Rmd
new file mode 100644
index 00000000..3d785225
--- /dev/null
+++ b/benchmarks/replications.Rmd
@@ -0,0 +1,86 @@
+---
+title: "Replications"
+author: "Jouni Helske"
+date: "11/17/2021"
+output: html_document
+---
+
+```{r srr-tags, eval = FALSE, echo = FALSE}
+#' @srrstats {G5.4, G5.4a, G5.4b, G5.4c, G5.5, G5.6, G5.6a, G5.6b, G5.7} The
+#' algorithms work correctly as per Vihola, Helske, Franks (2020)
+#' (all simulations were implemented with the bssm package) and Helske
+#' and Vihola (2021).
+
+```
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+## Reproducing some of the results of IS-MCMC paper
+
+Full simulation experiments of Vihola, Helske and Franks (2020) takes some time,
+here we only run a single replication to test that the methods work as expected.
+To generate Table 1 in the paper, the code below should be run say 1000 times,
+from which IREs could be computed.
+
+```{r, cache = TRUE}
+library("bssm")
+data(poisson_series)
+s <- sd(log(pmax(0.1, poisson_series)))
+model <- bsm_ng(poisson_series, sd_level = uniform(0.115, 0, 2 * s),
+ sd_slope = uniform(0.004, 0, 2 * s), P1 = diag(0.1, 2), distribution = "poisson")
+
+
+d <- data.frame(mean = NA, se = NA,
+ variable = c("sd_level", "sd_slope", "u_1", "u_100"),
+ mcmc_type = rep(c("approx", "da", "is1", "is2", "pm"),
+ times = 4*c(2, 6, 6, 6, 6)),
+ sampling_method = c(rep("psi", 8),
+ rep(rep(c("bsf", "spdk", "psi"), each = 2 * 4), 4)),
+ local_approx = rep(c(TRUE, FALSE), each = 4),
+ time = NA,
+ acceptance_rate = NA)
+
+iter <- 1e4 # Use less iterations than in the paper for faster experiment
+for(i in seq(1, nrow(d), by = 4)) {
+
+ cat("Testing method '", d$mcmc_type[i], "' with sampling by '",
+ d$sampling_method[i], "' and local_approx '", d$local_approx[i], "'\n",
+ sep = "")
+
+ res <- run_mcmc(model, iter = iter,
+ sampling_method = d$sampling_method[i],
+ particles = switch(d$sampling_method[i],
+ bsf = 200,
+ spdk = 10,
+ psi = 10),
+ mcmc_type = d$mcmc_type[i],
+ local_approx = d$local_approx[i],
+ end_adaptive_phase = TRUE)
+
+ w <- res$counts *
+ if (res$mcmc_type %in% paste0("is", 1:2)) res$weights else 1
+
+ d[((i - 1) + 1):((i - 1) + 4), "mean"] <- c(
+ diagis::weighted_mean(res$theta, w),
+ diagis::weighted_mean(res$alpha[1, 1, ], w),
+ diagis::weighted_mean(res$alpha[100, 1, ], w))
+
+ d[((i - 1) + 1):((i - 1) + 4), "se"] <- c(
+ sqrt(asymptotic_var(res$theta[, 1], w)),
+ sqrt(asymptotic_var(res$theta[, 2], w)),
+ sqrt(asymptotic_var(res$alpha[1, 1, ], w)),
+ sqrt(asymptotic_var(res$alpha[100, 1, ], w)))
+
+ d$time[((i - 1) + 1):((i - 1) + 4)] <- res$time["elapsed"]
+ d$acceptance_rate[((i - 1) + 1):((i - 1) + 4)] <- res$acceptance_rate
+}
+```
+Results:
+```{r}
+library(dplyr)
+d |>
+ arrange(local_approx, variable, mcmc_type, sampling_method)
+```
+
diff --git a/benchmarks/replications.html b/benchmarks/replications.html
new file mode 100644
index 00000000..abe8ed8f
--- /dev/null
+++ b/benchmarks/replications.html
@@ -0,0 +1,530 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/codecov.yml b/codecov.yml
new file mode 100644
index 00000000..e5a1ab93
--- /dev/null
+++ b/codecov.yml
@@ -0,0 +1,14 @@
+comment: false
+
+coverage:
+ status:
+ project:
+ default:
+ target: auto
+ threshold: 1%
+ informational: true
+ patch:
+ default:
+ target: auto
+ threshold: 1%
+ informational: true
diff --git a/codemeta.json b/codemeta.json
new file mode 100644
index 00000000..5ea2691d
--- /dev/null
+++ b/codemeta.json
@@ -0,0 +1,372 @@
+{
+ "@context": "https://doi.org/10.5063/schema/codemeta-2.0",
+ "@type": "SoftwareSourceCode",
+ "identifier": "bssm",
+ "description": "Efficient methods for Bayesian inference of state space models via Markov chain Monte Carlo (MCMC) based on parallel importance sampling type weighted estimators (Vihola, Helske, and Franks, 2020,
+
+
+
+
+
+
+Replications
+Jouni Helske
+11/17/2021
+ +
+
+
+
+
+
+Reproducing some of the results of IS-MCMC paper
+Full simulation experiments of Vihola, Helske and Franks (2020) takes some time, here we only run a single replication to test that the methods work as expected. To generate Table 1 in the paper, the code below should be run say 1000 times, from which IREs could be computed.
+library("bssm")##
+## Attaching package: 'bssm'## The following object is masked from 'package:base':
+##
+## gammadata(poisson_series)
+s <- sd(log(pmax(0.1, poisson_series)))
+model <- bsm_ng(poisson_series, sd_level = uniform(0.115, 0, 2 * s),
+ sd_slope = uniform(0.004, 0, 2 * s), P1 = diag(0.1, 2), distribution = "poisson")
+
+
+d <- data.frame(mean = NA, se = NA,
+ variable = c("sd_level", "sd_slope", "u_1", "u_100"),
+ mcmc_type = rep(c("approx", "da", "is1", "is2", "pm"),
+ times = 4*c(2, 6, 6, 6, 6)),
+ sampling_method = c(rep("psi", 8),
+ rep(rep(c("bsf", "spdk", "psi"), each = 2 * 4), 4)),
+ local_approx = rep(c(TRUE, FALSE), each = 4),
+ time = NA,
+ acceptance_rate = NA)
+
+iter <- 1e4 # Use less iterations than in the paper for faster experiment
+for(i in seq(1, nrow(d), by = 4)) {
+
+ cat("Testing method '", d$mcmc_type[i], "' with sampling by '",
+ d$sampling_method[i], "' and local_approx '", d$local_approx[i], "'\n",
+ sep = "")
+
+ res <- run_mcmc(model, iter = iter,
+ sampling_method = d$sampling_method[i],
+ particles = switch(d$sampling_method[i],
+ bsf = 200,
+ spdk = 10,
+ psi = 10),
+ mcmc_type = d$mcmc_type[i],
+ local_approx = d$local_approx[i],
+ end_adaptive_phase = TRUE)
+
+ w <- res$counts *
+ if (res$mcmc_type %in% paste0("is", 1:2)) res$weights else 1
+
+ d[((i - 1) + 1):((i - 1) + 4), "mean"] <- c(
+ diagis::weighted_mean(res$theta, w),
+ diagis::weighted_mean(res$alpha[1, 1, ], w),
+ diagis::weighted_mean(res$alpha[100, 1, ], w))
+
+ d[((i - 1) + 1):((i - 1) + 4), "se"] <- c(
+ sqrt(asymptotic_var(res$theta[, 1], w)),
+ sqrt(asymptotic_var(res$theta[, 2], w)),
+ sqrt(asymptotic_var(res$alpha[1, 1, ], w)),
+ sqrt(asymptotic_var(res$alpha[100, 1, ], w)))
+
+ d$time[((i - 1) + 1):((i - 1) + 4)] <- res$time["elapsed"]
+ d$acceptance_rate[((i - 1) + 1):((i - 1) + 4)] <- res$acceptance_rate
+}## Testing method 'approx' with sampling by 'psi' and local_approx 'TRUE'
+## Testing method 'approx' with sampling by 'psi' and local_approx 'FALSE'
+## Testing method 'da' with sampling by 'bsf' and local_approx 'TRUE'
+## Testing method 'da' with sampling by 'bsf' and local_approx 'FALSE'
+## Testing method 'da' with sampling by 'spdk' and local_approx 'TRUE'
+## Testing method 'da' with sampling by 'spdk' and local_approx 'FALSE'
+## Testing method 'da' with sampling by 'psi' and local_approx 'TRUE'
+## Testing method 'da' with sampling by 'psi' and local_approx 'FALSE'
+## Testing method 'is1' with sampling by 'bsf' and local_approx 'TRUE'
+## Testing method 'is1' with sampling by 'bsf' and local_approx 'FALSE'
+## Testing method 'is1' with sampling by 'spdk' and local_approx 'TRUE'
+## Testing method 'is1' with sampling by 'spdk' and local_approx 'FALSE'
+## Testing method 'is1' with sampling by 'psi' and local_approx 'TRUE'
+## Testing method 'is1' with sampling by 'psi' and local_approx 'FALSE'
+## Testing method 'is2' with sampling by 'bsf' and local_approx 'TRUE'
+## Testing method 'is2' with sampling by 'bsf' and local_approx 'FALSE'
+## Testing method 'is2' with sampling by 'spdk' and local_approx 'TRUE'
+## Testing method 'is2' with sampling by 'spdk' and local_approx 'FALSE'
+## Testing method 'is2' with sampling by 'psi' and local_approx 'TRUE'
+## Testing method 'is2' with sampling by 'psi' and local_approx 'FALSE'
+## Testing method 'pm' with sampling by 'bsf' and local_approx 'TRUE'
+## Testing method 'pm' with sampling by 'bsf' and local_approx 'FALSE'
+## Testing method 'pm' with sampling by 'spdk' and local_approx 'TRUE'
+## Testing method 'pm' with sampling by 'spdk' and local_approx 'FALSE'
+## Testing method 'pm' with sampling by 'psi' and local_approx 'TRUE'
+## Testing method 'pm' with sampling by 'psi' and local_approx 'FALSE'Results:
+library(dplyr)##
+## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
+##
+## filter, lag## The following objects are masked from 'package:base':
+##
+## intersect, setdiff, setequal, uniond %>%
+ arrange(local_approx, variable, mcmc_type, sampling_method)## mean se variable mcmc_type sampling_method local_approx
+## 1 0.09578772 0.0022141806 sd_level approx psi FALSE
+## 2 0.09183384 0.0025847529 sd_level da bsf FALSE
+## 3 0.09056407 0.0034324313 sd_level da psi FALSE
+## 4 0.09693717 0.0019153340 sd_level da spdk FALSE
+## 5 0.09562614 0.0025997452 sd_level is1 bsf FALSE
+## 6 0.09303018 0.0027318562 sd_level is1 psi FALSE
+## 7 0.09758895 0.0025748789 sd_level is1 spdk FALSE
+## 8 0.10105081 0.0025722112 sd_level is2 bsf FALSE
+## 9 0.09536685 0.0027073170 sd_level is2 psi FALSE
+## 10 0.09590957 0.0023060056 sd_level is2 spdk FALSE
+## 11 0.09698759 0.0038890507 sd_level pm bsf FALSE
+## 12 0.09430755 0.0024183804 sd_level pm psi FALSE
+## 13 0.08704774 0.0035083078 sd_level pm spdk FALSE
+## 14 0.01506997 0.0004887884 sd_slope approx psi FALSE
+## 15 0.01622993 0.0007867216 sd_slope da bsf FALSE
+## 16 0.01736272 0.0007476045 sd_slope da psi FALSE
+## 17 0.01542487 0.0004711415 sd_slope da spdk FALSE
+## 18 0.01533332 0.0005808372 sd_slope is1 bsf FALSE
+## 19 0.01653758 0.0006161702 sd_slope is1 psi FALSE
+## 20 0.01555065 0.0006237624 sd_slope is1 spdk FALSE
+## 21 0.01450235 0.0007518923 sd_slope is2 bsf FALSE
+## 22 0.01555210 0.0005458160 sd_slope is2 psi FALSE
+## 23 0.01522803 0.0005224722 sd_slope is2 spdk FALSE
+## 24 0.01526484 0.0008580442 sd_slope pm bsf FALSE
+## 25 0.01726606 0.0005463481 sd_slope pm psi FALSE
+## 26 0.01719711 0.0005954383 sd_slope pm spdk FALSE
+## 27 -0.07869788 0.0105559757 u_1 approx psi FALSE
+## 28 -0.06847499 0.0092537409 u_1 da bsf FALSE
+## 29 -0.07676795 0.0124369314 u_1 da psi FALSE
+## 30 -0.08390600 0.0098000700 u_1 da spdk FALSE
+## 31 -0.08702685 0.0125916067 u_1 is1 bsf FALSE
+## 32 -0.07410507 0.0095352137 u_1 is1 psi FALSE
+## 33 -0.08393639 0.0114658112 u_1 is1 spdk FALSE
+## 34 -0.08050907 0.0179799047 u_1 is2 bsf FALSE
+## 35 -0.07701605 0.0105821573 u_1 is2 psi FALSE
+## 36 -0.07846330 0.0121624205 u_1 is2 spdk FALSE
+## 37 -0.09249907 0.0106879249 u_1 pm bsf FALSE
+## 38 -0.08991847 0.0111100364 u_1 pm psi FALSE
+## 39 -0.06743786 0.0106694815 u_1 pm spdk FALSE
+## 40 2.63760608 0.0079707360 u_100 approx psi FALSE
+## 41 2.63182502 0.0080169407 u_100 da bsf FALSE
+## 42 2.62174704 0.0073447869 u_100 da psi FALSE
+## 43 2.61369859 0.0072760694 u_100 da spdk FALSE
+## 44 2.62253096 0.0086131210 u_100 is1 bsf FALSE
+## 45 2.62697081 0.0073913295 u_100 is1 psi FALSE
+## 46 2.61395948 0.0065218848 u_100 is1 spdk FALSE
+## 47 2.62980275 0.0117036903 u_100 is2 bsf FALSE
+## 48 2.62820881 0.0061786231 u_100 is2 psi FALSE
+## 49 2.61647003 0.0066978454 u_100 is2 spdk FALSE
+## 50 2.62452635 0.0071852884 u_100 pm bsf FALSE
+## 51 2.62372637 0.0064406343 u_100 pm psi FALSE
+## 52 2.61946902 0.0053898263 u_100 pm spdk FALSE
+## 53 0.09332580 0.0026097129 sd_level approx psi TRUE
+## 54 0.08850704 0.0034589396 sd_level da bsf TRUE
+## 55 0.09495943 0.0026074340 sd_level da psi TRUE
+## 56 0.09321799 0.0025949141 sd_level da spdk TRUE
+## 57 0.08969266 0.0024885812 sd_level is1 bsf TRUE
+## 58 0.09524581 0.0026432277 sd_level is1 psi TRUE
+## 59 0.09200344 0.0025662615 sd_level is1 spdk TRUE
+## 60 0.09164685 0.0029147370 sd_level is2 bsf TRUE
+## 61 0.09085042 0.0026602853 sd_level is2 psi TRUE
+## 62 0.09401393 0.0023809428 sd_level is2 spdk TRUE
+## 63 0.09946007 0.0029061242 sd_level pm bsf TRUE
+## 64 0.09425689 0.0026954951 sd_level pm psi TRUE
+## 65 0.09196438 0.0034327278 sd_level pm spdk TRUE
+## 66 0.01649875 0.0006084374 sd_slope approx psi TRUE
+## 67 0.01629020 0.0007204475 sd_slope da bsf TRUE
+## 68 0.01675024 0.0005625769 sd_slope da psi TRUE
+## 69 0.01596388 0.0005902601 sd_slope da spdk TRUE
+## 70 0.01684928 0.0005048439 sd_slope is1 bsf TRUE
+## 71 0.01619417 0.0006082412 sd_slope is1 psi TRUE
+## 72 0.01549983 0.0005977432 sd_slope is1 spdk TRUE
+## 73 0.01704707 0.0007564945 sd_slope is2 bsf TRUE
+## 74 0.01713414 0.0005378904 sd_slope is2 psi TRUE
+## 75 0.01581035 0.0006098587 sd_slope is2 spdk TRUE
+## 76 0.01493558 0.0007570220 sd_slope pm bsf TRUE
+## 77 0.01681022 0.0005775678 sd_slope pm psi TRUE
+## 78 0.01669934 0.0006237532 sd_slope pm spdk TRUE
+## 79 -0.04763190 0.0098184938 u_1 approx psi TRUE
+## 80 -0.06303314 0.0108945672 u_1 da bsf TRUE
+## 81 -0.07412832 0.0118281401 u_1 da psi TRUE
+## 82 -0.06864390 0.0107736811 u_1 da spdk TRUE
+## 83 -0.06413780 0.0102697767 u_1 is1 bsf TRUE
+## 84 -0.07305205 0.0087975786 u_1 is1 psi TRUE
+## 85 -0.05281387 0.0087550795 u_1 is1 spdk TRUE
+## 86 -0.06758727 0.0181045578 u_1 is2 bsf TRUE
+## 87 -0.08541492 0.0092005051 u_1 is2 psi TRUE
+## 88 -0.06533748 0.0107268023 u_1 is2 spdk TRUE
+## 89 -0.07094869 0.0127049274 u_1 pm bsf TRUE
+## 90 -0.08475790 0.0099826212 u_1 pm psi TRUE
+## 91 -0.07116972 0.0117171069 u_1 pm spdk TRUE
+## 92 2.62666325 0.0064783980 u_100 approx psi TRUE
+## 93 2.63482193 0.0078166866 u_100 da bsf TRUE
+## 94 2.61414177 0.0064441268 u_100 da psi TRUE
+## 95 2.61282885 0.0065739149 u_100 da spdk TRUE
+## 96 2.61354285 0.0064132143 u_100 is1 bsf TRUE
+## 97 2.61608214 0.0071682674 u_100 is1 psi TRUE
+## 98 2.62069831 0.0067225577 u_100 is1 spdk TRUE
+## 99 2.60570014 0.0103065461 u_100 is2 bsf TRUE
+## 100 2.62148516 0.0074049226 u_100 is2 psi TRUE
+## 101 2.61732719 0.0067673451 u_100 is2 spdk TRUE
+## 102 2.61310072 0.0074142079 u_100 pm bsf TRUE
+## 103 2.60707729 0.0068890297 u_100 pm psi TRUE
+## 104 2.60921701 0.0070487136 u_100 pm spdk TRUE
+## time acceptance_rate
+## 1 1.04 0.2336
+## 2 16.59 0.2466
+## 3 2.31 0.2068
+## 4 1.65 0.2202
+## 5 143.79 0.2348
+## 6 4.13 0.2238
+## 7 2.34 0.2262
+## 8 8.51 0.2300
+## 9 1.63 0.2470
+## 10 1.26 0.2396
+## 11 65.08 0.2388
+## 12 7.94 0.2354
+## 13 4.34 0.2410
+## 14 1.04 0.2336
+## 15 16.59 0.2466
+## 16 2.31 0.2068
+## 17 1.65 0.2202
+## 18 143.79 0.2348
+## 19 4.13 0.2238
+## 20 2.34 0.2262
+## 21 8.51 0.2300
+## 22 1.63 0.2470
+## 23 1.26 0.2396
+## 24 65.08 0.2388
+## 25 7.94 0.2354
+## 26 4.34 0.2410
+## 27 1.04 0.2336
+## 28 16.59 0.2466
+## 29 2.31 0.2068
+## 30 1.65 0.2202
+## 31 143.79 0.2348
+## 32 4.13 0.2238
+## 33 2.34 0.2262
+## 34 8.51 0.2300
+## 35 1.63 0.2470
+## 36 1.26 0.2396
+## 37 65.08 0.2388
+## 38 7.94 0.2354
+## 39 4.34 0.2410
+## 40 1.04 0.2336
+## 41 16.59 0.2466
+## 42 2.31 0.2068
+## 43 1.65 0.2202
+## 44 143.79 0.2348
+## 45 4.13 0.2238
+## 46 2.34 0.2262
+## 47 8.51 0.2300
+## 48 1.63 0.2470
+## 49 1.26 0.2396
+## 50 65.08 0.2388
+## 51 7.94 0.2354
+## 52 4.34 0.2410
+## 53 3.74 0.2404
+## 54 20.38 0.2408
+## 55 4.83 0.2250
+## 56 4.02 0.2448
+## 57 147.64 0.2434
+## 58 5.78 0.2432
+## 59 4.56 0.2322
+## 60 11.19 0.2510
+## 61 3.89 0.2384
+## 62 3.59 0.2342
+## 63 61.67 0.2428
+## 64 11.47 0.2358
+## 65 6.64 0.2358
+## 66 3.74 0.2404
+## 67 20.38 0.2408
+## 68 4.83 0.2250
+## 69 4.02 0.2448
+## 70 147.64 0.2434
+## 71 5.78 0.2432
+## 72 4.56 0.2322
+## 73 11.19 0.2510
+## 74 3.89 0.2384
+## 75 3.59 0.2342
+## 76 61.67 0.2428
+## 77 11.47 0.2358
+## 78 6.64 0.2358
+## 79 3.74 0.2404
+## 80 20.38 0.2408
+## 81 4.83 0.2250
+## 82 4.02 0.2448
+## 83 147.64 0.2434
+## 84 5.78 0.2432
+## 85 4.56 0.2322
+## 86 11.19 0.2510
+## 87 3.89 0.2384
+## 88 3.59 0.2342
+## 89 61.67 0.2428
+## 90 11.47 0.2358
+## 91 6.64 0.2358
+## 92 3.74 0.2404
+## 93 20.38 0.2408
+## 94 4.83 0.2250
+## 95 4.02 0.2448
+## 96 147.64 0.2434
+## 97 5.78 0.2432
+## 98 4.56 0.2322
+## 99 11.19 0.2510
+## 100 3.89 0.2384
+## 101 3.59 0.2342
+## 102 61.67 0.2428
+## 103 11.47 0.2358
+## 104 6.64 0.2358