diff --git a/R/boot.R b/R/boot.R
index 8b3cd2552..c2d3f9a9e 100644
--- a/R/boot.R
+++ b/R/boot.R
@@ -44,12 +44,12 @@ generate_data <- function(dist, data, args, weighted, censoring, parametric) {
   sample_parametric(dist, args = args, weighted = weighted, censoring = censoring)
 }
 
-boot_filename <- function(i, dist) {
-  paste0("boot_", stringr::str_pad(i, width = 9, pad = "0"), "_", dist, ".csv")
+boot_filename <- function(i, dist, prefix, ext) {
+  paste0(prefix, "_", stringr::str_pad(i, width = 9, pad = "0"), "_", dist, ext)
 }
 
-boot_filepath <- function(i, dist, save_to) {
-  file.path(save_to, boot_filename(i, dist))
+boot_filepath <- function(i, dist, save_to, prefix = "data", ext = ".csv") {
+  file.path(save_to, boot_filename(i, dist, prefix = prefix, ext = ext))
 }
 
 sample_parameters <- function(i, dist, fun, data, args, pars, weighted, censoring, min_pmix, range_shape1, range_shape2, parametric, control, save_to, wts = NULL) {
@@ -66,10 +66,12 @@ sample_parameters <- function(i, dist, fun, data, args, pars, weighted, censorin
     pars <- slnorm_lnorm(new_data)
   }
   if(dist == "multi") {
-    dist <- names(pars)
+    dist2 <- names(pars)
+  } else {
+    dist2 <- dist
   }
 
-  fit <- fun(new_data, dist,
+  fit <- fun(new_data, dist2,
     min_pmix = min_pmix, range_shape1 = range_shape1,
     range_shape2 = range_shape2, control = control, pars = pars, hessian = FALSE,
     censoring = censoring, weighted = weighted
@@ -79,6 +81,11 @@ sample_parameters <- function(i, dist, fun, data, args, pars, weighted, censorin
     return(NULL)
   }
   est <- estimates(fit, multi = TRUE)
+  
+  if(!is.null(save_to)) {
+    saveRDS(est, boot_filepath(i, dist, save_to, prefix = "estimates", ext = ".rds"))
+  }
+  
   if(!is.null(wts)) {
     est[names(wts)] <- unname(wts)
   }
@@ -93,13 +100,6 @@ boot_estimates <- function(fun, dist, estimates, pars, nboot, data, weighted, ce
 
   data <- data[c("left", "right", "weight")]
   
-  if(!is.null(save_to)) {
-    if(!requireNamespace("readr", quietly = TRUE)) {
-      err("Package 'readr' must be installed.")
-    }
-    readr::write_csv(data, boot_filepath(0, dist, save_to))
-  }
-  
   seeds <- seed_streams(nboot)
   
   if(fix_weights) {
@@ -108,6 +108,14 @@ boot_estimates <- function(fun, dist, estimates, pars, nboot, data, weighted, ce
     wts <- NULL
   }
   
+  if(!is.null(save_to)) {
+    if(!requireNamespace("readr", quietly = TRUE)) {
+      err("Package 'readr' must be installed.")
+    }
+    readr::write_csv(data, boot_filepath(0, dist, save_to))
+    saveRDS(estimates, boot_filepath(0, dist, save_to, prefix = "estimates", ext = ".rds"))
+  }
+  
   estimates <- future_map(1:nboot, sample_parameters,
     dist = dist, fun = sfun,
     data = data, args = args, pars = pars,
@@ -117,6 +125,6 @@ boot_estimates <- function(fun, dist, estimates, pars, nboot, data, weighted, ce
     wts = wts,
     .options = furrr::furrr_options(seed = seeds)
   )
-
+  
   estimates[!vapply(estimates, is.null, TRUE)]
 }
diff --git a/tests/testthat/_snaps/hc/hc.csv b/tests/testthat/_snaps/hc/hc.csv
index b3c5c25ea..c43c0e4b5 100644
--- a/tests/testthat/_snaps/hc/hc.csv
+++ b/tests/testthat/_snaps/hc/hc.csv
@@ -1,7 +1,7 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-gamma,5,1.07428,0.983875,0.709111,3.19472,0.356574,parametric,10
-lgumbel,5,1.76939,0.401763,1.39032,2.63032,0.0134466,parametric,10
-llogis,5,1.56226,0.834776,0.851639,3.37274,0.0656452,parametric,10
-lnorm,5,1.68117,0.823601,0.8894,3.02109,0.177236,parametric,10
-lnorm_lnorm,5,1.54141,0.329119,0.952154,1.82195,0.0296268,parametric,10
-weibull,5,1.08673,0.895134,0.819219,3.40858,0.357472,parametric,10
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+gamma,5,1.07428,0.983875,0.709111,3.19472,0.356574,parametric,10,1
+lgumbel,5,1.76939,0.401763,1.39032,2.63032,0.0134466,parametric,10,1
+llogis,5,1.56226,0.834776,0.851639,3.37274,0.0656452,parametric,10,1
+lnorm,5,1.68117,0.823601,0.8894,3.02109,0.177236,parametric,10,1
+lnorm_lnorm,5,1.54141,0.329119,0.952154,1.82195,0.0296268,parametric,10,1
+weibull,5,1.08673,0.895134,0.819219,3.40858,0.357472,parametric,10,1
diff --git a/tests/testthat/_snaps/hc/hc_30.csv b/tests/testthat/_snaps/hc/hc_30.csv
index 621e3096b..4577719d9 100644
--- a/tests/testthat/_snaps/hc/hc_30.csv
+++ b/tests/testthat/_snaps/hc/hc_30.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,0.9294,0.0164012,0.886164,0.940088,1,parametric,100
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,0.9294,0.0164012,0.886164,0.940088,1,parametric,100,0.93
diff --git a/tests/testthat/_snaps/hc/hc_boron.csv b/tests/testthat/_snaps/hc/hc_boron.csv
index 16668a516..8962f99f1 100644
--- a/tests/testthat/_snaps/hc/hc_boron.csv
+++ b/tests/testthat/_snaps/hc/hc_boron.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-invpareto,5,0.386944,0.347873,0.195882,1.10643,1,non-parametric,10
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+invpareto,5,0.386944,0.347873,0.195882,1.10643,1,non-parametric,10,1
diff --git a/tests/testthat/_snaps/hc/hc_burrIII3.csv b/tests/testthat/_snaps/hc/hc_burrIII3.csv
index 2b77a40b0..210258da0 100644
--- a/tests/testthat/_snaps/hc/hc_burrIII3.csv
+++ b/tests/testthat/_snaps/hc/hc_burrIII3.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-burrIII3,5,0.0180122,0.0330222,0.00257335,0.0960996,1,non-parametric,10
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+burrIII3,5,0.0180122,0.0330222,0.00257335,0.0960996,1,non-parametric,10,1
diff --git a/tests/testthat/_snaps/hc/hc_burrIII3_parametric.csv b/tests/testthat/_snaps/hc/hc_burrIII3_parametric.csv
index a6e8c5772..db4d23f56 100644
--- a/tests/testthat/_snaps/hc/hc_burrIII3_parametric.csv
+++ b/tests/testthat/_snaps/hc/hc_burrIII3_parametric.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-burrIII3,5,0.0180122,0.0149874,0.00621681,0.0503455,1,parametric,10
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+burrIII3,5,0.0180122,0.0149874,0.00621681,0.0503455,1,parametric,10,1
diff --git a/tests/testthat/_snaps/hc/hc_fix.csv b/tests/testthat/_snaps/hc/hc_fix.csv
index 6a27743f5..e3271124d 100644
--- a/tests/testthat/_snaps/hc/hc_fix.csv
+++ b/tests/testthat/_snaps/hc/hc_fix.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68947,0.528952,0.921861,2.95419,1,parametric,100
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68947,0.528952,0.921861,2.95419,1,parametric,100,1
diff --git a/tests/testthat/_snaps/hc/hc_fixmulti.csv b/tests/testthat/_snaps/hc/hc_fixmulti.csv
new file mode 100644
index 000000000..47c9b026d
--- /dev/null
+++ b/tests/testthat/_snaps/hc/hc_fixmulti.csv
@@ -0,0 +1,3 @@
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68947,0.528952,0.921861,2.95419,1,parametric,100,1
+average,10,2.62057,0.716268,1.55146,4.31158,1,parametric,100,1
diff --git a/tests/testthat/_snaps/hc/hc_notallestimates.csv b/tests/testthat/_snaps/hc/hc_notallestimates.csv
new file mode 100644
index 000000000..479c4ae38
--- /dev/null
+++ b/tests/testthat/_snaps/hc/hc_notallestimates.csv
@@ -0,0 +1,2 @@
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.64901,0.384673,1.36309,2.46218,1,parametric,10,0.9
diff --git a/tests/testthat/_snaps/hc/hc_save_to.csv b/tests/testthat/_snaps/hc/hc_save_to.csv
index 17a8e2be1..2b5c72b72 100644
--- a/tests/testthat/_snaps/hc/hc_save_to.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to1.csv b/tests/testthat/_snaps/hc/hc_save_to1.csv
index 17a8e2be1..2b5c72b72 100644
--- a/tests/testthat/_snaps/hc/hc_save_to1.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to1.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to11.csv b/tests/testthat/_snaps/hc/hc_save_to11.csv
index abaed854c..4835d8d25 100644
--- a/tests/testthat/_snaps/hc/hc_save_to11.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to11.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68117,NA,1.10678,1.10678,1,parametric,1
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68117,NA,1.10678,1.10678,1,parametric,1,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to1_not_multi.csv b/tests/testthat/_snaps/hc/hc_save_to1_not_multi.csv
index bb683fb75..58b45d6dd 100644
--- a/tests/testthat/_snaps/hc/hc_save_to1_not_multi.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to1_not_multi.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68117,0.52617,1.30753,2.28691,1,parametric,3
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68117,0.52617,1.30753,2.28691,1,parametric,3,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to1_not_multi_default.csv b/tests/testthat/_snaps/hc/hc_save_to1_not_multi_default.csv
index da3c7f866..a118151a8 100644
--- a/tests/testthat/_snaps/hc/hc_save_to1_not_multi_default.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to1_not_multi_default.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.24152,NA,1.62613,1.62613,1,parametric,1
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.24152,NA,1.62613,1.62613,1,parametric,1,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to1_rescale.csv b/tests/testthat/_snaps/hc/hc_save_to1_rescale.csv
index 17a8e2be1..2b5c72b72 100644
--- a/tests/testthat/_snaps/hc/hc_save_to1_rescale.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to1_rescale.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to1data.csv b/tests/testthat/_snaps/hc/hc_save_to1data.csv
index 17a8e2be1..2b5c72b72 100644
--- a/tests/testthat/_snaps/hc/hc_save_to1data.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to1data.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to_not_multi.csv b/tests/testthat/_snaps/hc/hc_save_to_not_multi.csv
index bb683fb75..58b45d6dd 100644
--- a/tests/testthat/_snaps/hc/hc_save_to_not_multi.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to_not_multi.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68117,0.52617,1.30753,2.28691,1,parametric,3
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68117,0.52617,1.30753,2.28691,1,parametric,3,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to_not_multi_default.csv b/tests/testthat/_snaps/hc/hc_save_to_not_multi_default.csv
index da3c7f866..a118151a8 100644
--- a/tests/testthat/_snaps/hc/hc_save_to_not_multi_default.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to_not_multi_default.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.24152,NA,1.62613,1.62613,1,parametric,1
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.24152,NA,1.62613,1.62613,1,parametric,1,1
diff --git a/tests/testthat/_snaps/hc/hc_save_to_rescale.csv b/tests/testthat/_snaps/hc/hc_save_to_rescale.csv
index 17a8e2be1..2b5c72b72 100644
--- a/tests/testthat/_snaps/hc/hc_save_to_rescale.csv
+++ b/tests/testthat/_snaps/hc/hc_save_to_rescale.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68117,0.515645,1.13096,2.11007,1,parametric,3,1
diff --git a/tests/testthat/_snaps/hc/hc_unfix.csv b/tests/testthat/_snaps/hc/hc_unfix.csv
index 303dcd2c8..45b365be3 100644
--- a/tests/testthat/_snaps/hc/hc_unfix.csv
+++ b/tests/testthat/_snaps/hc/hc_unfix.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.68947,0.592316,1.02258,3.29178,1,parametric,100
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68947,0.592316,1.02258,3.29178,1,parametric,100,1
diff --git a/tests/testthat/_snaps/hc/hc_unfixmulti.csv b/tests/testthat/_snaps/hc/hc_unfixmulti.csv
new file mode 100644
index 000000000..f99b60d09
--- /dev/null
+++ b/tests/testthat/_snaps/hc/hc_unfixmulti.csv
@@ -0,0 +1,3 @@
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.68947,0.592316,1.02258,3.29178,1,parametric,100,1
+average,10,2.62057,0.74108,1.65497,4.39447,1,parametric,100,1
diff --git a/tests/testthat/_snaps/hp/hc_fix.csv b/tests/testthat/_snaps/hp/hc_fix.csv
index 6698790f7..9e2d626cf 100644
--- a/tests/testthat/_snaps/hp/hc_fix.csv
+++ b/tests/testthat/_snaps/hp/hc_fix.csv
@@ -1,2 +1,2 @@
-dist,conc,est,se,lcl,ucl,wt,method,nboot
-average,1,1.87688,1.5501,0.365892,5.53175,1,parametric,100
+dist,conc,est,se,lcl,ucl,wt,method,nboot,pboot
+average,1,1.87688,1.5501,0.365892,5.53175,1,parametric,100,1
diff --git a/tests/testthat/_snaps/hp/hc_unfix.csv b/tests/testthat/_snaps/hp/hc_unfix.csv
index cbba18078..8a7e38f8a 100644
--- a/tests/testthat/_snaps/hp/hc_unfix.csv
+++ b/tests/testthat/_snaps/hp/hc_unfix.csv
@@ -1,2 +1,2 @@
-dist,conc,est,se,lcl,ucl,wt,method,nboot
-average,1,1.87688,1.37401,0.0959776,4.81697,1,parametric,100
+dist,conc,est,se,lcl,ucl,wt,method,nboot,pboot
+average,1,1.87688,1.37401,0.0959776,4.81697,1,parametric,100,1
diff --git a/tests/testthat/_snaps/hp/hp_30.csv b/tests/testthat/_snaps/hp/hp_30.csv
index 534a34012..7a7a41c56 100644
--- a/tests/testthat/_snaps/hp/hp_30.csv
+++ b/tests/testthat/_snaps/hp/hp_30.csv
@@ -1,2 +1,2 @@
-dist,conc,est,se,lcl,ucl,wt,method,nboot
-average,1,14.0388,3.71053,7.63528,26.3366,1,parametric,100
+dist,conc,est,se,lcl,ucl,wt,method,nboot,pboot
+average,1,14.0388,3.71053,7.63528,26.3366,1,parametric,100,0.93
diff --git a/tests/testthat/_snaps/zzz-unstable/hc_err.csv b/tests/testthat/_snaps/zzz-unstable/hc_err.csv
index d094488b1..e85b66090 100644
--- a/tests/testthat/_snaps/zzz-unstable/hc_err.csv
+++ b/tests/testthat/_snaps/zzz-unstable/hc_err.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,0.934681,0.33094,0.839423,2.35185,1,parametric,100
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,0.934681,0.33094,0.839423,2.35185,1,parametric,100,0.97
diff --git a/tests/testthat/_snaps/zzz-unstable/hc_err_na.csv b/tests/testthat/_snaps/zzz-unstable/hc_err_na.csv
index 4c5ba460d..83998a928 100644
--- a/tests/testthat/_snaps/zzz-unstable/hc_err_na.csv
+++ b/tests/testthat/_snaps/zzz-unstable/hc_err_na.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,0.934681,NA,NA,NA,1,parametric,100
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,0.934681,NA,NA,NA,1,parametric,100,0.94
diff --git a/tests/testthat/_snaps/zzz-unstable/hc_nonpara.csv b/tests/testthat/_snaps/zzz-unstable/hc_nonpara.csv
index 6560dd04e..7c5562f1c 100644
--- a/tests/testthat/_snaps/zzz-unstable/hc_nonpara.csv
+++ b/tests/testthat/_snaps/zzz-unstable/hc_nonpara.csv
@@ -1,2 +1,2 @@
-dist,percent,est,se,lcl,ucl,wt,method,nboot
-average,5,1.42635,0.0298623,1.37466,1.45809,1,non-parametric,10
+dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
+average,5,1.42635,0.0298623,1.37466,1.45809,1,non-parametric,10,1
diff --git a/tests/testthat/_snaps/zzz-unstable/hp_err.csv b/tests/testthat/_snaps/zzz-unstable/hp_err.csv
index 0c8ada7f2..0b6ace1ea 100644
--- a/tests/testthat/_snaps/zzz-unstable/hp_err.csv
+++ b/tests/testthat/_snaps/zzz-unstable/hp_err.csv
@@ -1,2 +1,2 @@
-dist,conc,est,se,lcl,ucl,wt,method,nboot
-average,1,9.18236,5.44367,2.45953e-42,21.8494,1,parametric,100
+dist,conc,est,se,lcl,ucl,wt,method,nboot,pboot
+average,1,9.18236,5.44367,2.45953e-42,21.8494,1,parametric,100,0.97
diff --git a/tests/testthat/_snaps/zzz-unstable/hp_err_na.csv b/tests/testthat/_snaps/zzz-unstable/hp_err_na.csv
index 8f2cccb55..ddaed4c8c 100644
--- a/tests/testthat/_snaps/zzz-unstable/hp_err_na.csv
+++ b/tests/testthat/_snaps/zzz-unstable/hp_err_na.csv
@@ -1,2 +1,2 @@
-dist,conc,est,se,lcl,ucl,wt,method,nboot
-average,1,9.18236,NA,NA,NA,1,parametric,100
+dist,conc,est,se,lcl,ucl,wt,method,nboot,pboot
+average,1,9.18236,NA,NA,NA,1,parametric,100,0.94
diff --git a/tests/testthat/_snaps/zzz-unstable/hp_nonpara.csv b/tests/testthat/_snaps/zzz-unstable/hp_nonpara.csv
index 30593174b..8b1a7e3d9 100644
--- a/tests/testthat/_snaps/zzz-unstable/hp_nonpara.csv
+++ b/tests/testthat/_snaps/zzz-unstable/hp_nonpara.csv
@@ -1,2 +1,2 @@
-dist,conc,est,se,lcl,ucl,wt,method,nboot
-average,1,2.27801,0.129041,2.14227,2.50138,1,non-parametric,10
+dist,conc,est,se,lcl,ucl,wt,method,nboot,pboot
+average,1,2.27801,0.129041,2.14227,2.50138,1,non-parametric,10,1
diff --git a/tests/testthat/_snaps/zzz-unstable/hp_para.csv b/tests/testthat/_snaps/zzz-unstable/hp_para.csv
index bbfc3ebdd..e5267af21 100644
--- a/tests/testthat/_snaps/zzz-unstable/hp_para.csv
+++ b/tests/testthat/_snaps/zzz-unstable/hp_para.csv
@@ -1,2 +1,2 @@
-dist,conc,est,se,lcl,ucl,wt,method,nboot
-average,1,2.27801,0.0639425,2.2218,2.41308,1,parametric,10
+dist,conc,est,se,lcl,ucl,wt,method,nboot,pboot
+average,1,2.27801,0.0639425,2.2218,2.41308,1,parametric,10,1
diff --git a/tests/testthat/test-hc.R b/tests/testthat/test-hc.R
index b337c8431..e6d66ab8a 100644
--- a/tests/testthat/test-hc.R
+++ b/tests/testthat/test-hc.R
@@ -18,7 +18,7 @@ test_that("hc", {
   set.seed(102)
   hc <- ssd_hc(fits, ci = TRUE, nboot = 10, average = FALSE)
   expect_s3_class(hc, "tbl")
-  expect_snapshot_boot_data(hc, "hc")
+  expect_snapshot_data(hc, "hc")
 })
 
 test_that("ssd_hc list must be named", {
@@ -404,7 +404,7 @@ test_that("ssd_hc cis with non-convergence", {
   expect_identical(attr(fit, "min_pmix"), 0.3)
   hc30 <- ssd_hc(fit, ci = TRUE, nboot = 100, min_pboot = 0.9, multi = FALSE)
   expect_s3_class(hc30, "tbl")
-  expect_snapshot_boot_data(hc30, "hc_30")
+  expect_snapshot_data(hc30, "hc_30")
 })
 
 test_that("ssd_hc cis with error and multiple dists", {
@@ -448,7 +448,7 @@ test_that("ssd_hc_burrlioz gets estimates with invpareto", {
   fit <- ssd_fit_burrlioz(ssddata::ccme_boron)
   set.seed(47)
   hc_boron <- ssd_hc(fit, nboot = 10, ci = TRUE, min_pboot = 0)
-  expect_snapshot_boot_data(hc_boron, "hc_boron")
+  expect_snapshot_data(hc_boron, "hc_boron")
 })
 
 test_that("ssd_hc_burrlioz gets estimates with burrIII3", {
@@ -459,7 +459,7 @@ test_that("ssd_hc_burrlioz gets estimates with burrIII3", {
   expect_identical(names(fit), "burrIII3")
   set.seed(49)
   hc_burrIII3 <- ssd_hc(fit, nboot = 10, ci = TRUE, min_pboot = 0)
-  expect_snapshot_boot_data(hc_burrIII3, "hc_burrIII3")
+  expect_snapshot_data(hc_burrIII3, "hc_burrIII3")
 })
 
 test_that("ssd_hc_burrlioz gets estimates with burrIII3 parametric", {
@@ -473,7 +473,7 @@ test_that("ssd_hc_burrlioz gets estimates with burrIII3 parametric", {
                         nboot = 10, ci = TRUE, min_pboot = 0,
                         parametric = TRUE
   )
-  expect_snapshot_boot_data(hc_burrIII3, "hc_burrIII3_parametric")
+  expect_snapshot_data(hc_burrIII3, "hc_burrIII3_parametric")
 })
 
 test_that("ssd_hc passing all boots ccme_chloride lnorm_lnorm", {
@@ -495,12 +495,30 @@ test_that("ssd_hc save_to", {
   fits <- ssd_fit_dists(ssddata::ccme_boron, dist = "lnorm")
   set.seed(102)
   hc <- ssd_hc(fits, nboot = 3, ci = TRUE, save_to = dir)
-  expect_snapshot_boot_data(hc, "hc_save_to")
-  expect_identical(list.files(dir), c("boot_000000000_multi.csv", "boot_000000001_multi.csv", "boot_000000002_multi.csv", "boot_000000003_multi.csv"))
-  data <- read.csv(file.path(dir, "boot_000000000_multi.csv"))
-  expect_snapshot_boot_data(hc, "hc_save_to1data")
-  boot1 <- read.csv(file.path(dir, "boot_000000001_multi.csv"))
-  expect_snapshot_boot_data(hc, "hc_save_to1")
+  expect_snapshot_data(hc, "hc_save_to")
+  expect_identical(list.files(dir), c("data_000000000_multi.csv", "data_000000001_multi.csv", "data_000000002_multi.csv", 
+                                      "data_000000003_multi.csv", "estimates_000000000_multi.rds", 
+                                      "estimates_000000001_multi.rds", "estimates_000000002_multi.rds", 
+                                      "estimates_000000003_multi.rds"))
+  data <- read.csv(file.path(dir, "data_000000000_multi.csv"))
+  expect_snapshot_data(hc, "hc_save_to1data")
+  boot1 <- read.csv(file.path(dir, "data_000000001_multi.csv"))
+  expect_snapshot_data(hc, "hc_save_to1")
+  ests <- readRDS(file.path(dir, "estimates_000000000_multi.rds"))
+  ests1 <- readRDS(file.path(dir, "estimates_000000001_multi.rds"))
+  
+  expect_identical(names(ests), names(ests1))
+  expect_identical(names(ests), c("burrIII3.weight", "burrIII3.shape1", "burrIII3.shape2", "burrIII3.scale", 
+                                  "gamma.weight", "gamma.shape", "gamma.scale", "gompertz.weight", 
+                                  "gompertz.location", "gompertz.shape", "invpareto.weight", "invpareto.shape", 
+                                  "invpareto.scale", "lgumbel.weight", "lgumbel.locationlog", "lgumbel.scalelog", 
+                                  "llogis.weight", "llogis.locationlog", "llogis.scalelog", "llogis_llogis.weight", 
+                                  "llogis_llogis.locationlog1", "llogis_llogis.scalelog1", "llogis_llogis.locationlog2", 
+                                  "llogis_llogis.scalelog2", "llogis_llogis.pmix", "lnorm.weight", 
+                                  "lnorm.meanlog", "lnorm.sdlog", "lnorm_lnorm.weight", "lnorm_lnorm.meanlog1", 
+                                  "lnorm_lnorm.sdlog1", "lnorm_lnorm.meanlog2", "lnorm_lnorm.sdlog2", 
+                                  "lnorm_lnorm.pmix", "weibull.weight", "weibull.shape", "weibull.scale"
+  ))
 })
 
 test_that("ssd_hc save_to multi = FALSE", {
@@ -509,10 +527,13 @@ test_that("ssd_hc save_to multi = FALSE", {
   fits <- ssd_fit_dists(ssddata::ccme_boron, dist = "lnorm")
   set.seed(102)
   hc <- ssd_hc(fits, nboot = 3, ci = TRUE, save_to = dir, multi = FALSE)
-  expect_snapshot_boot_data(hc, "hc_save_to_not_multi")
-  expect_identical(list.files(dir), c("boot_000000000_lnorm.csv", "boot_000000001_lnorm.csv", "boot_000000002_lnorm.csv", "boot_000000003_lnorm.csv"))
-  boot1 <- read.csv(file.path(dir, "boot_000000001_lnorm.csv"))
-  expect_snapshot_boot_data(hc, "hc_save_to1_not_multi")
+  expect_snapshot_data(hc, "hc_save_to_not_multi")
+  expect_identical(list.files(dir), c("data_000000000_lnorm.csv", "data_000000001_lnorm.csv", "data_000000002_lnorm.csv", 
+                                      "data_000000003_lnorm.csv", "estimates_000000000_lnorm.rds", 
+                                      "estimates_000000001_lnorm.rds", "estimates_000000002_lnorm.rds", 
+                                      "estimates_000000003_lnorm.rds"))
+  data1 <- read.csv(file.path(dir, "data_000000001_lnorm.csv"))
+  expect_snapshot_data(hc, "hc_save_to1_not_multi")
 })
 
 test_that("ssd_hc save_to multi = FALSE default", {
@@ -521,15 +542,21 @@ test_that("ssd_hc save_to multi = FALSE default", {
   fits <- ssd_fit_dists(ssddata::ccme_boron)
   set.seed(102)
   hc <- ssd_hc(fits, nboot = 1, ci = TRUE, save_to = dir, multi = FALSE)
-  expect_snapshot_boot_data(hc, "hc_save_to_not_multi_default")
+  expect_snapshot_data(hc, "hc_save_to_not_multi_default")
   expect_identical(sort(list.files(dir)), 
-                   sort(c("boot_000000000_gamma.csv", "boot_000000000_lgumbel.csv", "boot_000000000_llogis.csv", 
-                          "boot_000000000_lnorm_lnorm.csv", "boot_000000000_lnorm.csv", 
-                          "boot_000000000_weibull.csv", "boot_000000001_gamma.csv", "boot_000000001_lgumbel.csv", 
-                          "boot_000000001_llogis.csv", "boot_000000001_lnorm_lnorm.csv", 
-                          "boot_000000001_lnorm.csv", "boot_000000001_weibull.csv")))
-  boot1 <- read.csv(file.path(dir, "boot_000000001_lnorm.csv"))
-  expect_snapshot_boot_data(hc, "hc_save_to1_not_multi_default")
+                   sort(c("data_000000000_gamma.csv", "data_000000000_lgumbel.csv", "data_000000000_llogis.csv", 
+                          "data_000000000_lnorm_lnorm.csv", "data_000000000_lnorm.csv", 
+                          "data_000000000_weibull.csv", "data_000000001_gamma.csv", "data_000000001_lgumbel.csv", 
+                          "data_000000001_llogis.csv", "data_000000001_lnorm_lnorm.csv", 
+                          "data_000000001_lnorm.csv", "data_000000001_weibull.csv", "estimates_000000000_gamma.rds", 
+                          "estimates_000000000_lgumbel.rds", "estimates_000000000_llogis.rds", 
+                          "estimates_000000000_lnorm_lnorm.rds", "estimates_000000000_lnorm.rds", 
+                          "estimates_000000000_weibull.rds", "estimates_000000001_gamma.rds", 
+                          "estimates_000000001_lgumbel.rds", "estimates_000000001_llogis.rds", 
+                          "estimates_000000001_lnorm_lnorm.rds", "estimates_000000001_lnorm.rds", 
+                          "estimates_000000001_weibull.rds")))
+  boot1 <- read.csv(file.path(dir, "data_000000001_lnorm.csv"))
+  expect_snapshot_data(hc, "hc_save_to1_not_multi_default")
 })
 
 test_that("ssd_hc save_to rescale", {
@@ -538,10 +565,13 @@ test_that("ssd_hc save_to rescale", {
   fits <- ssd_fit_dists(ssddata::ccme_boron, dist = "lnorm", rescale = TRUE)
   set.seed(102)
   hc <- ssd_hc(fits, nboot = 3, ci = TRUE, save_to = dir)
-  expect_snapshot_boot_data(hc, "hc_save_to_rescale")
-  expect_identical(list.files(dir), c("boot_000000000_multi.csv", "boot_000000001_multi.csv", "boot_000000002_multi.csv", "boot_000000003_multi.csv"))
-  boot1 <- read.csv(file.path(dir, "boot_000000001_multi.csv"))
-  expect_snapshot_boot_data(hc, "hc_save_to1_rescale")
+  expect_snapshot_data(hc, "hc_save_to_rescale")
+  expect_identical(list.files(dir), c("data_000000000_multi.csv", "data_000000001_multi.csv", "data_000000002_multi.csv", 
+                                      "data_000000003_multi.csv", "estimates_000000000_multi.rds", 
+                                      "estimates_000000001_multi.rds", "estimates_000000002_multi.rds", 
+                                      "estimates_000000003_multi.rds"))
+  boot1 <- read.csv(file.path(dir, "data_000000001_multi.csv"))
+  expect_snapshot_data(hc, "hc_save_to1_rescale")
 })
 
 test_that("ssd_hc save_to lnorm 1", {
@@ -550,9 +580,10 @@ test_that("ssd_hc save_to lnorm 1", {
   fits <- ssd_fit_dists(ssddata::ccme_boron, dist = "lnorm")
   set.seed(102)
   hc <- ssd_hc(fits, nboot = 1, ci = TRUE, save_to = dir)
-  expect_snapshot_boot_data(hc, "hc_save_to11")
-  expect_identical(list.files(dir), c("boot_000000000_multi.csv", "boot_000000001_multi.csv"))
-  boot1 <- read.csv(file.path(dir, "boot_000000001_multi.csv"))
+  expect_snapshot_data(hc, "hc_save_to11")
+  expect_identical(list.files(dir), c("data_000000000_multi.csv", "data_000000001_multi.csv", "estimates_000000000_multi.rds", 
+                                      "estimates_000000001_multi.rds"))
+  boot1 <- read.csv(file.path(dir, "data_000000001_multi.csv"))
   fit1 <- ssd_fit_dists(boot1, dist = "lnorm", left = "left", right = "right", weight = "weight")
   est <- ssd_hc(fit1)$est
   expect_identical(hc$lcl, est)
@@ -565,13 +596,15 @@ test_that("ssd_hc save_to replaces", {
   fits <- ssd_fit_dists(ssddata::ccme_boron, dist = "lnorm")
   set.seed(102)
   hc <- ssd_hc(fits, nboot = 1, ci = TRUE, save_to = dir)
-  expect_identical(list.files(dir), c("boot_000000000_multi.csv", "boot_000000001_multi.csv"))
-  boot <- read.csv(file.path(dir, "boot_000000001_multi.csv"))
+  expect_identical(list.files(dir), c("data_000000000_multi.csv", "data_000000001_multi.csv", "estimates_000000000_multi.rds", 
+                                      "estimates_000000001_multi.rds"))
+  boot <- read.csv(file.path(dir, "data_000000001_multi.csv"))
   hc2 <- ssd_hc(fits, nboot = 1, ci = TRUE, save_to = dir)
-  expect_identical(list.files(dir), c("boot_000000000_multi.csv", "boot_000000001_multi.csv"))
-  boot2 <- read.csv(file.path(dir, "boot_000000001_multi.csv"))
-  expect_snapshot_boot_data(boot, "hc_boot1_replace")
-  expect_snapshot_boot_data(boot2, "hc_boot2_replace")
+  expect_identical(list.files(dir), c("data_000000000_multi.csv", "data_000000001_multi.csv", "estimates_000000000_multi.rds", 
+                                      "estimates_000000001_multi.rds"))
+  boot2 <- read.csv(file.path(dir, "data_000000001_multi.csv"))
+  expect_snapshot_data(boot, "hc_boot1_replace")
+  expect_snapshot_data(boot2, "hc_boot2_replace")
 })
 
 test_that("ssd_hc fix_weight", {
@@ -579,9 +612,66 @@ test_that("ssd_hc fix_weight", {
   
   set.seed(102)
   hc_unfix <- ssd_hc(fits, nboot = 100, ci = TRUE, fix_weights = FALSE)
-  expect_snapshot_boot_data(hc_unfix, "hc_unfix")
+  expect_snapshot_data(hc_unfix, "hc_unfix")
   
   set.seed(102)
   hc_fix <- ssd_hc(fits, nboot = 100, ci = TRUE, fix_weights = TRUE)
-  expect_snapshot_boot_data(hc_fix, "hc_fix")
+  expect_snapshot_data(hc_fix, "hc_fix")
+})
+
+test_that("ssd_hc multiple values", {
+  fits <- ssd_fit_dists(ssddata::ccme_boron, dist = c("lnorm", "lgumbel"))
+  
+  set.seed(102)
+  hc_unfix <- ssd_hc(fits, percent = c(5,10), nboot = 100, ci = TRUE, fix_weights = FALSE)
+  expect_snapshot_data(hc_unfix, "hc_unfixmulti")
+  
+  set.seed(102)
+  hc_fix <- ssd_hc(fits, percent = c(5,10), nboot = 100, ci = TRUE, fix_weights = TRUE)
+  expect_snapshot_data(hc_fix, "hc_fixmulti")
+})
+
+test_that("ssd_hc multiple values save_to", {
+  dir <- withr::local_tempdir()
+  
+  fits <- ssd_fit_dists(ssddata::ccme_boron, dist = c("lnorm", "lgumbel"))
+  
+  set.seed(102)
+  hc <- ssd_hc(fits, percent = c(5,10), nboot = 2, save_to = dir, ci = TRUE)
+  expect_identical(list.files(dir), c("data_000000000_multi.csv", "data_000000001_multi.csv", "data_000000002_multi.csv", 
+                                      "estimates_000000000_multi.rds", "estimates_000000001_multi.rds", 
+                                      "estimates_000000002_multi.rds"))
+})
+
+test_that("ssd_hc not multi save_to", {
+  dir <- withr::local_tempdir()
+  
+  fits <- ssd_fit_dists(ssddata::ccme_boron, dist = c("lnorm", "lgumbel"))
+  
+  set.seed(102)
+  hc <- ssd_hc(fits, nboot = 2, multi = FALSE, save_to = dir, ci = TRUE)
+  expect_identical(list.files(dir), c("data_000000000_lgumbel.csv", "data_000000000_lnorm.csv", "data_000000001_lgumbel.csv", 
+                                      "data_000000001_lnorm.csv", "data_000000002_lgumbel.csv", "data_000000002_lnorm.csv", 
+                                      "estimates_000000000_lgumbel.rds", "estimates_000000000_lnorm.rds", 
+                                      "estimates_000000001_lgumbel.rds", "estimates_000000001_lnorm.rds", 
+                                      "estimates_000000002_lgumbel.rds", "estimates_000000002_lnorm.rds"
+  ))
+})
+
+test_that("not all estimates if fail", {
+  dir <- withr::local_tempdir()
+  
+  fit <- ssd_fit_dists(ssddata::ccme_boron, dists = c("lnorm", "lnorm_lnorm"))
+  set.seed(49)
+  hc <- ssd_hc(fit, nboot = 10, ci = TRUE,
+                        parametric = TRUE, save_to = dir, min_pboot = 0.8)
+  expect_snapshot_data(hc, "hc_notallestimates")
+  expect_identical(list.files(dir), c("data_000000000_multi.csv", "data_000000001_multi.csv", "data_000000002_multi.csv", 
+                                      "data_000000003_multi.csv", "data_000000004_multi.csv", "data_000000005_multi.csv", 
+                                      "data_000000006_multi.csv", "data_000000007_multi.csv", "data_000000008_multi.csv", 
+                                      "data_000000009_multi.csv", "data_000000010_multi.csv", "estimates_000000000_multi.rds", 
+                                      "estimates_000000001_multi.rds", "estimates_000000002_multi.rds", 
+                                      "estimates_000000003_multi.rds", "estimates_000000004_multi.rds", 
+                                      "estimates_000000005_multi.rds", "estimates_000000006_multi.rds", 
+                                      "estimates_000000007_multi.rds", "estimates_000000009_multi.rds", "estimates_000000010_multi.rds"))
 })
diff --git a/tests/testthat/test-hp.R b/tests/testthat/test-hp.R
index b32ef5c19..54048a391 100644
--- a/tests/testthat/test-hp.R
+++ b/tests/testthat/test-hp.R
@@ -286,7 +286,7 @@ test_that("ssd_hp cis with non-convergence", {
   expect_identical(attr(fit, "min_pmix"), 0.3)
   hp30 <- ssd_hp(fit, conc = 1, ci = TRUE, nboot = 100, min_pboot = 0.9, multi = FALSE)
   expect_s3_class(hp30, "tbl")
-  expect_snapshot_boot_data(hp30, "hp_30")
+  expect_snapshot_data(hp30, "hp_30")
 })
 
 test_that("ssd_hp cis with error and multiple dists", {
@@ -323,10 +323,10 @@ test_that("ssd_hp fix_weight", {
   
   set.seed(102)
   hc_unfix <- ssd_hp(fits, nboot = 100, ci = TRUE, fix_weights = FALSE)
-  expect_snapshot_boot_data(hc_unfix, "hc_unfix")
+  expect_snapshot_data(hc_unfix, "hc_unfix")
   
   set.seed(102)
   hc_fix <- ssd_hp(fits, nboot = 100, ci = TRUE, fix_weights = TRUE)
-  expect_snapshot_boot_data(hc_fix, "hc_fix")
+  expect_snapshot_data(hc_fix, "hc_fix")
 })
 
diff --git a/tests/testthat/test-zzz-unstable.R b/tests/testthat/test-zzz-unstable.R
index 25836667a..266d2f5b2 100644
--- a/tests/testthat/test-zzz-unstable.R
+++ b/tests/testthat/test-zzz-unstable.R
@@ -201,10 +201,10 @@ test_that("ssd_hc cis with error", {
   expect_identical(attr(fit, "min_pmix"), 0.1)
   expect_warning(hc_err <- ssd_hc(fit, ci = TRUE, nboot = 100))
   expect_s3_class(hc_err, "tbl")
-  expect_snapshot_boot_data(hc_err, "hc_err_na")
+  expect_snapshot_data(hc_err, "hc_err_na")
   hc_err <- ssd_hc(fit, ci = TRUE, nboot = 100, min_pboot = 0.92, multi = FALSE)
   expect_s3_class(hc_err, "tbl")
-  expect_snapshot_boot_data(hc_err, "hc_err")
+  expect_snapshot_data(hc_err, "hc_err")
 })
 
 test_that("ssd_hc comparable parametric and non-parametric big sample size", {
@@ -219,7 +219,7 @@ test_that("ssd_hc comparable parametric and non-parametric big sample size", {
   expect_snapshot_data(hc_para, "hc_para")
   set.seed(10)
   hc_nonpara <- ssd_hc(fit, ci = TRUE, nboot = 10, parametric = FALSE, multi = FALSE)
-  expect_snapshot_boot_data(hc_nonpara, "hc_nonpara")
+  expect_snapshot_data(hc_nonpara, "hc_nonpara")
 })
 
 test_that("ssd_hp cis with error", {
@@ -233,10 +233,10 @@ test_that("ssd_hp cis with error", {
   expect_identical(attr(fit, "min_pmix"), 0.1)
   expect_warning(hp_err <- ssd_hp(fit, conc = 1, ci = TRUE, nboot = 100))
   expect_s3_class(hp_err, "tbl")
-  expect_snapshot_boot_data(hp_err, "hp_err_na")
+  expect_snapshot_data(hp_err, "hp_err_na")
   hp_err <- ssd_hp(fit, conc = 1, ci = TRUE, nboot = 100, min_pboot = 0.92, multi = FALSE)
   expect_s3_class(hp_err, "tbl")
-  expect_snapshot_boot_data(hp_err, "hp_err")
+  expect_snapshot_data(hp_err, "hp_err")
 })
 
 test_that("ssd_hp comparable parametric and non-parametric big sample size", {
@@ -248,10 +248,10 @@ test_that("ssd_hp comparable parametric and non-parametric big sample size", {
   fit <- ssd_fit_dists(data, dists = "lnorm")
   set.seed(10)
   hp_para <- ssd_hp(fit, 1, ci = TRUE, nboot = 10, multi = FALSE)
-  expect_snapshot_boot_data(hp_para, "hp_para")
+  expect_snapshot_data(hp_para, "hp_para")
   set.seed(10)
   hp_nonpara <- ssd_hp(fit, 1, ci = TRUE, nboot = 10, parametric = FALSE, multi = FALSE)
-  expect_snapshot_boot_data(hp_nonpara, "hp_nonpara")
+  expect_snapshot_data(hp_nonpara, "hp_nonpara")
 })
 
 test_that("plot geoms", {