diff --git a/.github/workflows/Test-coverage.yaml b/.github/workflows/Test-coverage.yaml
index fd80e891..7681de61 100644
--- a/.github/workflows/Test-coverage.yaml
+++ b/.github/workflows/Test-coverage.yaml
@@ -63,6 +63,6 @@ jobs:
       - name: Test coverage (Windows)
         if: runner.os == 'Windows'
         run: |
-          options(covr.gcov = 'C:/rtools40/mingw64/bin/gcov.exe');
+          options(covr.gcov = 'C:/rtools44/mingw64/bin/gcov.exe');
           covr::codecov(type = "tests", function_exclusions = "sample_mpi")
         shell: Rscript {0}
diff --git a/.github/workflows/cmdstan-tarball-check.yaml b/.github/workflows/cmdstan-tarball-check.yaml
index d9869f78..c79fae4f 100644
--- a/.github/workflows/cmdstan-tarball-check.yaml
+++ b/.github/workflows/cmdstan-tarball-check.yaml
@@ -23,7 +23,7 @@ jobs:
       matrix:
         config:
           - {os: macOS-latest, r: 'release', rtools: ''}
-          - {os: windows-latest, r: 'release', rtools: '42'}
+          - {os: windows-latest, r: 'release', rtools: '44'}
           - {os: ubuntu-20.04, r: 'release', rtools: ''}
     env:
       R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
diff --git a/DESCRIPTION b/DESCRIPTION
index 192bb489..4aedc431 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -10,6 +10,7 @@ Authors@R:
       person(given = "Andrew", family = "Johnson", role = c("aut", "cre"),
              email = "andrew.johnson@arjohnsonau.com",
              comment = c(ORCID = "0000-0001-7000-8065")),
+      person(given = "Steve", family = "Bronder", role = "aut"),
       person(given = "Ben", family = "Bales", role = "ctb"),
       person(given = "Mitzi", family = "Morris", role = "ctb"),
       person(given = "Mikhail", family = "Popov", role = "ctb"),
@@ -18,7 +19,8 @@ Authors@R:
            email = "will.landau@gmail.com", comment = c(ORCID = "0000-0003-1878-3253")),
       person(given = "Jacob", family = "Socolar", role = "ctb"),
       person(given = "Martin", family = "Modrák", role = "ctb"),
-      person(given = "Steve", family = "Bronder", role = "ctb"))
+      person(given = "Ven", family = "Popov", role = "ctb")
+      )
 Description: A lightweight interface to 'Stan' <https://mc-stan.org>.
     The 'CmdStanR' interface is an alternative to 'RStan' that calls the command
     line interface for compilation and running algorithms instead of interfacing
diff --git a/NAMESPACE b/NAMESPACE
index 559452fb..b157025d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -6,6 +6,15 @@ S3method(as_draws,CmdStanMCMC)
 S3method(as_draws,CmdStanMLE)
 S3method(as_draws,CmdStanPathfinder)
 S3method(as_draws,CmdStanVB)
+S3method(process_init,"function")
+S3method(process_init,CmdStanLaplace)
+S3method(process_init,CmdStanMCMC)
+S3method(process_init,CmdStanMLE)
+S3method(process_init,CmdStanPathfinder)
+S3method(process_init,CmdStanVB)
+S3method(process_init,default)
+S3method(process_init,draws)
+S3method(process_init,list)
 export(as.CmdStanDiagnose)
 export(as.CmdStanGQ)
 export(as.CmdStanLaplace)
@@ -40,3 +49,4 @@ export(write_stan_json)
 export(write_stan_tempfile)
 import(R6)
 importFrom(posterior,as_draws)
+importFrom(stats,aggregate)
diff --git a/R/args.R b/R/args.R
index 7b4a8686..c22ebc7a 100644
--- a/R/args.R
+++ b/R/args.R
@@ -77,11 +77,12 @@ CmdStanArgs <- R6::R6Class(
       }
       self$output_dir <- repair_path(self$output_dir)
       self$output_basename <- output_basename
-      if (is.function(init)) {
-        init <- process_init_function(init, length(self$proc_ids), model_variables)
-      } else if (is.list(init) && !is.data.frame(init)) {
-        init <- process_init_list(init, length(self$proc_ids), model_variables)
+      if (inherits(self$method_args, "PathfinderArgs")) {
+        num_inits <- self$method_args$num_paths
+      } else {
+        num_inits <- length(self$proc_ids)
       }
+      init <- process_init(init, num_inits, model_variables)
       self$init <- init
       self$opencl_ids <- opencl_ids
       self$num_threads = NULL
@@ -691,7 +692,12 @@ validate_cmdstan_args <- function(self) {
     assert_file_exists(self$data_file, access = "r")
   }
   num_procs <- length(self$proc_ids)
-  validate_init(self$init, num_procs)
+  if (inherits(self$method_args, "PathfinderArgs")) {
+    num_inits <- self$method_args$num_paths
+  } else {
+    num_inits <- length(self$proc_ids)
+  }
+  validate_init(self$init, num_inits)
   validate_seed(self$seed, num_procs)
   if (!is.null(self$opencl_ids)) {
     if (cmdstan_version() < "2.26") {
@@ -1018,18 +1024,96 @@ validate_exe_file <- function(exe_file) {
   invisible(TRUE)
 }
 
+
+#' Generic for processing inits
+#' @noRd
+process_init <- function(init, ...) {
+  UseMethod("process_init")
+}
+
+#' Default method
+#' @noRd
+#' @export
+process_init.default <- function(init, ...) {
+  return(init)
+}
+
+#' Write initial values to files if provided as posterior `draws` object
+#' @noRd
+#' @param init A type that inherits the `posterior::draws` class.
+#' @param num_procs Number of inits requested
+#' @param model_variables  A list of all parameters with their types and
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
+#' @param warn_partial Should a warning be thrown if inits are only specified
+#'   for a subset of parameters? Can be controlled by global option
+#'   `cmdstanr_warn_inits`.
+#' @return A character vector of file paths.
+#' @export
+process_init.draws <- function(init, num_procs, model_variables = NULL,
+                               warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                               ...) {
+  if (!is.null(model_variables)) {
+    variable_names = names(model_variables$parameters)
+  } else {
+    variable_names = colnames(draws)[!grepl("__", colnames(draws))]
+  }
+  draws <- posterior::as_draws_df(init)
+  # Since all other process_init functions return `num_proc` inits
+  # This will only happen if a raw draws object is passed
+  if (nrow(draws) < num_procs) {
+    idx <- rep(1:nrow(draws), ceiling(num_procs / nrow(draws)))[1:num_procs]
+    draws <- draws[idx,]
+  } else if (nrow(draws) > num_procs) {
+    draws <- posterior::resample_draws(draws, ndraws = num_procs,
+                                       method ="simple_no_replace")
+  }
+  draws_rvar = posterior::as_draws_rvars(draws)
+  variable_names <- variable_names[variable_names %in% names(draws_rvar)]
+  draws_rvar <- posterior::subset_draws(draws_rvar, variable = variable_names)
+  inits = lapply(1:num_procs, function(draw_iter) {
+    init_i = lapply(variable_names, function(var_name) {
+      x = drop(posterior::draws_of(drop(
+        posterior::subset_draws(draws_rvar[[var_name]], draw=draw_iter))))
+      return(x)
+    })
+    bad_names = unlist(lapply(variable_names, function(var_name) {
+      x = drop(posterior::draws_of(drop(
+        posterior::subset_draws(draws_rvar[[var_name]], draw=draw_iter))))
+      if (any(is.infinite(x)) || any(is.na(x))) {
+        return(var_name)
+      }
+      return("")
+    }))
+    any_na_or_inf = bad_names != ""
+    if (any(any_na_or_inf)) {
+      err_msg = paste0(paste(bad_names[any_na_or_inf], collapse = ", "), " contains NA or Inf values!")
+      if (length(any_na_or_inf) > 1) {
+        err_msg = paste0("Variables: ", err_msg)
+      } else {
+        err_msg = paste0("Variable: ", err_msg)
+      }
+      stop(err_msg)
+    }
+    names(init_i) = variable_names
+    return(init_i)
+  })
+  return(process_init(inits, num_procs, model_variables, warn_partial))
+}
+
 #' Write initial values to files if provided as list of lists
 #' @noRd
 #' @param init List of init lists.
-#' @param num_procs Number of CmdStan processes.
+#' @param num_procs Number of inits needed.
 #' @param model_variables  A list of all parameters with their types and
-#'   number of dimensions. Typically the output of model$variables().
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
 #' @param warn_partial Should a warning be thrown if inits are only specified
 #'   for a subset of parameters? Can be controlled by global option
 #'   `cmdstanr_warn_inits`.
 #' @return A character vector of file paths.
-process_init_list <- function(init, num_procs, model_variables = NULL,
-                              warn_partial = getOption("cmdstanr_warn_inits", TRUE)) {
+#' @export
+process_init.list <- function(init, num_procs, model_variables = NULL,
+                              warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                              ...) {
   if (!all(sapply(init, function(x) is.list(x) && !is.data.frame(x)))) {
     stop("If 'init' is a list it must be a list of lists.", call. = FALSE)
   }
@@ -1083,10 +1167,11 @@ process_init_list <- function(init, num_procs, model_variables = NULL,
   }
   init_paths <-
     tempfile(
-      pattern = paste0("init-", seq_along(init), "-"),
+      pattern = "init-",
       tmpdir = cmdstan_tempdir(),
-      fileext = ".json"
+      fileext = ""
     )
+  init_paths <- paste0(init_paths, "_", seq_along(init), ".json")
   for (i in seq_along(init)) {
     write_stan_json(init[[i]], init_paths[i])
   }
@@ -1096,11 +1181,14 @@ process_init_list <- function(init, num_procs, model_variables = NULL,
 #' Write initial values to files if provided as function
 #' @noRd
 #' @param init Function generating a single list of initial values.
-#' @param num_procs Number of CmdStan processes.
+#' @param num_procs Number of inits needed.
 #' @param model_variables A list of all parameters with their types and
-#'   number of dimensions. Typically the output of model$variables().
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
 #' @return A character vector of file paths.
-process_init_function <- function(init, num_procs, model_variables = NULL) {
+#' @export
+process_init.function <- function(init, num_procs, model_variables = NULL,
+                                  warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                                  ...) {
   args <- formals(init)
   if (is.null(args)) {
     fn_test <- init()
@@ -1116,7 +1204,193 @@ process_init_function <- function(init, num_procs, model_variables = NULL) {
   if (!is.list(fn_test) || is.data.frame(fn_test)) {
     stop("If 'init' is a function it must return a single list.")
   }
-  process_init_list(init_list, num_procs, model_variables)
+  process_init(init_list, num_procs, model_variables)
+}
+
+#' Validate a fit is a valid init
+#' @noRd
+validate_fit_init = function(init, model_variables) {
+  # Convert from data.table to data.frame
+  if (all(init$return_codes() == 1)) {
+    stop("We are unable to create initial values from a model with no samples. Please check the results of the model used for inits before continuing.")
+  } else if (!is.null(model_variables) &&!any(names(model_variables$parameters) %in% init$metadata()$stan_variables)) {
+    stop("None of the names of the parameters for the model used for initial values match the names of parameters from the model currently running.")
+  }
+}
+
+#' Write initial values to files if provided as a `CmdStanMCMC` class
+#' @noRd
+#' @param init A `CmdStanMCMC` class
+#' @param num_procs Number of inits requested
+#' @param model_variables  A list of all parameters with their types and
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
+#' @param warn_partial Should a warning be thrown if inits are only specified
+#'   for a subset of parameters? Can be controlled by global option
+#'   `cmdstanr_warn_inits`.
+#' @return A character vector of file paths.
+#' @export
+process_init.CmdStanMCMC <- function(init, num_procs, model_variables = NULL,
+                                     warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                                     ...) {
+  validate_fit_init(init, model_variables)
+  draws_df = init$draws(format = "df")
+  if (is.null(model_variables)) {
+    model_variables = list(parameters = colnames(draws_df)[2:(length(colnames(draws_df)) - 3)])
+  }
+  init_draws_df = posterior::resample_draws(draws_df, ndraws = num_procs,
+                                            method = "simple_no_replace")
+  init_draws_lst = process_init(init_draws_df,
+                                num_procs = num_procs, model_variables = model_variables)
+  return(init_draws_lst)
+}
+
+#' Performs PSIS resampling on the draws from an approxmation method for inits.
+#' @noRd
+#' @param init A set of draws with `lp__` and `lp_approx__` columns.
+#' @param num_procs Number of inits requested
+#' @param model_variables  A list of all parameters with their types and
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
+#' @param warn_partial Should a warning be thrown if inits are only specified
+#'   for a subset of parameters? Can be controlled by global option
+#'   `cmdstanr_warn_inits`.
+#' @return A character vector of file paths.
+#' @importFrom stats aggregate
+process_init_approx <- function(init, num_procs, model_variables = NULL,
+                                warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                                ...) {
+  validate_fit_init(init, model_variables)
+  # Convert from data.table to data.frame
+  draws_df = init$draws(format = "df")
+  if (is.null(model_variables)) {
+    model_variables = list(parameters = colnames(draws_df)[3:(length(colnames(draws_df)) - 3)])
+  }
+  draws_df$lw = draws_df$lp__ - draws_df$lp_approx__
+  # Calculate unique draws based on 'lw' using base R functions
+  unique_draws = length(unique(draws_df$lw))
+  if (num_procs > unique_draws) {
+    if (inherits(init, " CmdStanPathfinder ")) {
+      algo_name = " Pathfinder "
+      extra_msg = " Try running Pathfinder with psis_resample=FALSE."
+    } else if (inherits(init, "CmdStanVB")) {
+      algo_name = " CmdStanVB "
+      extra_msg = ""
+    } else if (inherits(init, " CmdStanLaplace ")) {
+      algo_name = " CmdStanLaplace "
+      extra_msg = ""
+    } else {
+      algo_name = ""
+      extra_msg = ""
+    }
+    stop(paste0("Not enough distinct draws (", num_procs, ") in", algo_name ,
+      "fit to create inits.", extra_msg))
+  }
+  if (unique_draws < (0.95 * nrow(draws_df))) {
+    temp_df = stats::aggregate(.draw ~ lw, data = draws_df, FUN = min)
+    draws_df = posterior::as_draws_df(merge(temp_df, draws_df, by = 'lw'))
+    draws_df$weight = exp(draws_df$lw - max(draws_df$lw))
+  } else {
+      draws_df$weight = posterior::pareto_smooth(
+        exp(draws_df$lw - max(draws_df$lw)), tail = "right", return_k=FALSE)
+  }
+  init_draws_df = posterior::resample_draws(draws_df, ndraws = num_procs,
+                                            weights = draws_df$weight, method = "simple_no_replace")
+  init_draws_lst = process_init(init_draws_df,
+                                num_procs = num_procs, model_variables = model_variables, warn_partial)
+  return(init_draws_lst)
+}
+
+
+#' Write initial values to files if provided as a `CmdStanPathfinder` class
+#' @noRd
+#' @param init A `CmdStanPathfinder` class
+#' @param num_procs Number of inits requested
+#' @param model_variables  A list of all parameters with their types and
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
+#' @param warn_partial Should a warning be thrown if inits are only specified
+#'   for a subset of parameters? Can be controlled by global option
+#'   `cmdstanr_warn_inits`.
+#' @return A character vector of file paths.
+#' @export
+process_init.CmdStanPathfinder <- function(init, num_procs, model_variables = NULL,
+                                           warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                                           ...) {
+  if (!init$metadata()$calculate_lp) {
+    validate_fit_init(init, model_variables)
+    # Convert from data.table to data.frame
+    draws_df = init$draws(format = "df")
+    if (is.null(model_variables)) {
+      model_variables = list(parameters = colnames(draws_df)[3:(length(colnames(draws_df)) - 3)])
+    }
+    draws_df$weight = rep(1.0, nrow(draws_df))
+    init_draws_df = posterior::resample_draws(draws_df, ndraws = num_procs,
+      weights = draws_df$weight, method = "simple_no_replace")
+    init_draws_lst = process_init(init_draws_df,
+      num_procs = num_procs, model_variables = model_variables, warn_partial)
+    return(init_draws_lst)
+  } else {
+    process_init_approx(init, num_procs, model_variables, warn_partial)
+  }
+}
+
+#' Write initial values to files if provided as a `CmdStanVB` class
+#' @noRd
+#' @param init A `CmdStanVB` class
+#' @param num_procs Number of inits requested
+#' @param model_variables  A list of all parameters with their types and
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
+#' @param warn_partial Should a warning be thrown if inits are only specified
+#'   for a subset of parameters? Can be controlled by global option
+#'   `cmdstanr_warn_inits`.
+#' @return A character vector of file paths.
+#' @export
+process_init.CmdStanVB <- function(init, num_procs, model_variables = NULL,
+                                   warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                                   ...) {
+  process_init_approx(init, num_procs, model_variables, warn_partial)
+}
+
+#' Write initial values to files if provided as a `CmdStanLaplace` class
+#' @noRd
+#' @param init A `CmdStanLaplace` class
+#' @param num_procs Number of inits requested
+#' @param model_variables  A list of all parameters with their types and
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
+#' @param warn_partial Should a warning be thrown if inits are only specified
+#'   for a subset of parameters? Can be controlled by global option
+#'   `cmdstanr_warn_inits`.
+#' @return A character vector of file paths.
+#' @export
+process_init.CmdStanLaplace <- function(init, num_procs, model_variables = NULL,
+                                        warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                                        ...) {
+  process_init_approx(init, num_procs, model_variables, warn_partial)
+}
+
+
+#' Write initial values to files if provided as a `CmdStanMLE` class
+#' @noRd
+#' @param init A `CmdStanMLE` class
+#' @param num_procs Number of inits requested
+#' @param model_variables  A list of all parameters with their types and
+#'   number of dimensions. Typically the output of `model$variables()$parameters`.
+#' @param warn_partial Should a warning be thrown if inits are only specified
+#'   for a subset of parameters? Can be controlled by global option
+#'   `cmdstanr_warn_inits`.
+#' @return A character vector of file paths.
+#' @export
+process_init.CmdStanMLE <- function(init, num_procs, model_variables = NULL,
+                                    warn_partial = getOption("cmdstanr_warn_inits", TRUE),
+                                    ...) {
+  # Convert from data.table to data.frame
+  validate_fit_init(init, model_variables)
+  draws_df = init$draws(format = "df")
+  if (is.null(model_variables)) {
+    model_variables = list(parameters = colnames(draws_df)[2:(length(colnames(draws_df)) - 3)])
+  }
+  init_draws_df = draws_df[rep(1, num_procs),]
+  init_draws_lst_lst = process_init(init_draws_df,
+                                    num_procs = num_procs, model_variables = model_variables, warn_partial)
+  return(init_draws_lst_lst)
 }
 
 #' Validate initial values
diff --git a/R/cmdstanr-package.R b/R/cmdstanr-package.R
index ce6bc525..00455cab 100644
--- a/R/cmdstanr-package.R
+++ b/R/cmdstanr-package.R
@@ -33,6 +33,6 @@
 #' @inherit cmdstan_model examples
 #' @import R6
 #'
-NULL
+"_PACKAGE"
 
 if (getRversion() >= "2.15.1")  utils::globalVariables(c("self", "private", "super"))
diff --git a/R/csv.R b/R/csv.R
index f2e674f8..16fe4807 100644
--- a/R/csv.R
+++ b/R/csv.R
@@ -774,6 +774,7 @@ read_csv_metadata <- function(csv_file) {
           sampling_time <- as.numeric(tmp)
         } else if (grepl("(Total)", tmp, fixed = TRUE)) {
           tmp <- gsub("seconds (Total)", "", tmp, fixed = TRUE)
+          tmp <- trimws(gsub(" Elapsed Time: ", "", tmp, fixed = TRUE))
           total_time <- as.numeric(tmp)
         }
         if (!is.null(csv_file_info$method) &&
diff --git a/R/fit.R b/R/fit.R
index 99feca4c..68215608 100644
--- a/R/fit.R
+++ b/R/fit.R
@@ -518,11 +518,11 @@ unconstrain_variables <- function(variables) {
          " not provided!", call. = FALSE)
   }
 
-  # Remove zero-length parameters from model_variables, otherwise process_init_list
+  # Remove zero-length parameters from model_variables, otherwise process_init
   # warns about missing inputs
   model_variables$parameters <- model_variables$parameters[nonzero_length_params]
 
-  stan_pars <- process_init_list(list(variables), num_procs = 1, model_variables)
+  stan_pars <- process_init(list(variables), num_procs = 1, model_variables)
   private$model_methods_env_$unconstrain_variables(private$model_methods_env_$model_ptr_, stan_pars)
 }
 CmdStanFit$set("public", name = "unconstrain_variables", value = unconstrain_variables)
@@ -594,7 +594,7 @@ unconstrain_draws <- function(files = NULL, draws = NULL,
   # but not in metadata()$variables
   nonzero_length_params <- names(model_variables$parameters) %in% model_par_names
 
-  # Remove zero-length parameters from model_variables, otherwise process_init_list
+  # Remove zero-length parameters from model_variables, otherwise process_init
   # warns about missing inputs
   pars <- names(model_variables$parameters[nonzero_length_params])
 
diff --git a/man-roxygen/model-common-args.R b/man-roxygen/model-common-args.R
index f65eb21c..6b60da19 100644
--- a/man-roxygen/model-common-args.R
+++ b/man-roxygen/model-common-args.R
@@ -44,6 +44,31 @@
 #'  has argument `chain_id` it will be supplied with the chain id (from 1 to
 #'  number of chains) when called to generate the initial values. See
 #'  **Examples**.
+#'  * A [`CmdStanMCMC`], [`CmdStanMLE`], [`CmdStanVB`], [`CmdStanPathfinder`],
+#'  or [`CmdStanLaplace`] fit object.
+#'  If the fit object's parameters are only a subset of the model
+#'  parameters then the other parameters will be drawn by Stan's default
+#'  initialization. The fit object must have at least some parameters that are the
+#'  same name and dimensions as the current Stan model. For the `sample` and
+#'  `pathfinder` method, if the fit object has fewer draws than the requested
+#'  number of chains/paths then the inits will be drawn using sampling with
+#'  replacement. Otherwise sampling without replacement will be used.
+#'  When a [`CmdStanPathfinder`] fit object is used as the init, if
+#'. `psis_resample` was set to `FALSE` and `calculate_lp` was
+#'  set to `TRUE` (default), then resampling without replacement with Pareto
+#'  smoothed weights will be used. If `psis_resample` was set to `TRUE` or
+#'  `calculate_lp` was set to `FALSE` then sampling without replacement with
+#'  uniform weights will be used to select the draws.
+#'  PSIS resampling is used to select the draws for  [`CmdStanVB`],
+#'  and [`CmdStanLaplace`] fit objects.
+#'
+#'  * A type inheriting from `posterior::draws`. If the draws object has less
+#'  samples than the number of requested chains/paths then the inits will be
+#'  drawn using sampling with replacement. Otherwise sampling without
+#'  replacement will be used. If the draws object's parameters are only a subset
+#'  of the model parameters then the other parameters will be drawn by Stan's
+#'  default initialization. The fit object must have at least some parameters
+#'  that are the same name and dimensions as the current Stan model.
 #'
 #' @param save_latent_dynamics (logical) Should auxiliary diagnostic information
 #'   about the latent dynamics be written to temporary diagnostic CSV files?
diff --git a/man/cmdstanr-package.Rd b/man/cmdstanr-package.Rd
index ea43d565..739213ff 100644
--- a/man/cmdstanr-package.Rd
+++ b/man/cmdstanr-package.Rd
@@ -217,6 +217,7 @@ Authors:
 \itemize{
   \item Jonah Gabry \email{jsg2201@columbia.edu}
   \item Rok Češnovar \email{rok.cesnovar@fri.uni-lj.si}
+  \item Steve Bronder
 }
 
 Other contributors:
@@ -228,7 +229,7 @@ Other contributors:
   \item William Michael Landau \email{will.landau@gmail.com} (\href{https://orcid.org/0000-0003-1878-3253}{ORCID}) [contributor]
   \item Jacob Socolar [contributor]
   \item Martin Modrák [contributor]
-  \item Steve Bronder [contributor]
+  \item Ven Popov [contributor]
 }
 
 }
diff --git a/man/model-method-diagnose.Rd b/man/model-method-diagnose.Rd
index 631ce6f6..1deab91d 100644
--- a/man/model-method-diagnose.Rd
+++ b/man/model-method-diagnose.Rd
@@ -61,6 +61,30 @@ take no arguments or a single argument \code{chain_id}. For MCMC, if the functio
 has argument \code{chain_id} it will be supplied with the chain id (from 1 to
 number of chains) when called to generate the initial values. See
 \strong{Examples}.
+\item A \code{\link{CmdStanMCMC}}, \code{\link{CmdStanMLE}}, \code{\link{CmdStanVB}}, \code{\link{CmdStanPathfinder}},
+or \code{\link{CmdStanLaplace}} fit object.
+If the fit object's parameters are only a subset of the model
+parameters then the other parameters will be drawn by Stan's default
+initialization. The fit object must have at least some parameters that are the
+same name and dimensions as the current Stan model. For the \code{sample} and
+\code{pathfinder} method, if the fit object has fewer draws than the requested
+number of chains/paths then the inits will be drawn using sampling with
+replacement. Otherwise sampling without replacement will be used.
+When a \code{\link{CmdStanPathfinder}} fit object is used as the init, if
+. \code{psis_resample} was set to \code{FALSE} and \code{calculate_lp} was
+set to \code{TRUE} (default), then resampling without replacement with Pareto
+smoothed weights will be used. If \code{psis_resample} was set to \code{TRUE} or
+\code{calculate_lp} was set to \code{FALSE} then sampling without replacement with
+uniform weights will be used to select the draws.
+PSIS resampling is used to select the draws for  \code{\link{CmdStanVB}},
+and \code{\link{CmdStanLaplace}} fit objects.
+\item A type inheriting from \code{posterior::draws}. If the draws object has less
+samples than the number of requested chains/paths then the inits will be
+drawn using sampling with replacement. Otherwise sampling without
+replacement will be used. If the draws object's parameters are only a subset
+of the model parameters then the other parameters will be drawn by Stan's
+default initialization. The fit object must have at least some parameters
+that are the same name and dimensions as the current Stan model.
 }}
 
 \item{output_dir}{(string) A path to a directory where CmdStan should write
diff --git a/man/model-method-laplace.Rd b/man/model-method-laplace.Rd
index 2ebb00b3..0bb8bce8 100644
--- a/man/model-method-laplace.Rd
+++ b/man/model-method-laplace.Rd
@@ -74,6 +74,30 @@ take no arguments or a single argument \code{chain_id}. For MCMC, if the functio
 has argument \code{chain_id} it will be supplied with the chain id (from 1 to
 number of chains) when called to generate the initial values. See
 \strong{Examples}.
+\item A \code{\link{CmdStanMCMC}}, \code{\link{CmdStanMLE}}, \code{\link{CmdStanVB}}, \code{\link{CmdStanPathfinder}},
+or \code{\link{CmdStanLaplace}} fit object.
+If the fit object's parameters are only a subset of the model
+parameters then the other parameters will be drawn by Stan's default
+initialization. The fit object must have at least some parameters that are the
+same name and dimensions as the current Stan model. For the \code{sample} and
+\code{pathfinder} method, if the fit object has fewer draws than the requested
+number of chains/paths then the inits will be drawn using sampling with
+replacement. Otherwise sampling without replacement will be used.
+When a \code{\link{CmdStanPathfinder}} fit object is used as the init, if
+. \code{psis_resample} was set to \code{FALSE} and \code{calculate_lp} was
+set to \code{TRUE} (default), then resampling without replacement with Pareto
+smoothed weights will be used. If \code{psis_resample} was set to \code{TRUE} or
+\code{calculate_lp} was set to \code{FALSE} then sampling without replacement with
+uniform weights will be used to select the draws.
+PSIS resampling is used to select the draws for  \code{\link{CmdStanVB}},
+and \code{\link{CmdStanLaplace}} fit objects.
+\item A type inheriting from \code{posterior::draws}. If the draws object has less
+samples than the number of requested chains/paths then the inits will be
+drawn using sampling with replacement. Otherwise sampling without
+replacement will be used. If the draws object's parameters are only a subset
+of the model parameters then the other parameters will be drawn by Stan's
+default initialization. The fit object must have at least some parameters
+that are the same name and dimensions as the current Stan model.
 }}
 
 \item{save_latent_dynamics}{Ignored for this method.}
diff --git a/man/model-method-optimize.Rd b/man/model-method-optimize.Rd
index dbed1150..3718dac3 100644
--- a/man/model-method-optimize.Rd
+++ b/man/model-method-optimize.Rd
@@ -80,6 +80,30 @@ take no arguments or a single argument \code{chain_id}. For MCMC, if the functio
 has argument \code{chain_id} it will be supplied with the chain id (from 1 to
 number of chains) when called to generate the initial values. See
 \strong{Examples}.
+\item A \code{\link{CmdStanMCMC}}, \code{\link{CmdStanMLE}}, \code{\link{CmdStanVB}}, \code{\link{CmdStanPathfinder}},
+or \code{\link{CmdStanLaplace}} fit object.
+If the fit object's parameters are only a subset of the model
+parameters then the other parameters will be drawn by Stan's default
+initialization. The fit object must have at least some parameters that are the
+same name and dimensions as the current Stan model. For the \code{sample} and
+\code{pathfinder} method, if the fit object has fewer draws than the requested
+number of chains/paths then the inits will be drawn using sampling with
+replacement. Otherwise sampling without replacement will be used.
+When a \code{\link{CmdStanPathfinder}} fit object is used as the init, if
+. \code{psis_resample} was set to \code{FALSE} and \code{calculate_lp} was
+set to \code{TRUE} (default), then resampling without replacement with Pareto
+smoothed weights will be used. If \code{psis_resample} was set to \code{TRUE} or
+\code{calculate_lp} was set to \code{FALSE} then sampling without replacement with
+uniform weights will be used to select the draws.
+PSIS resampling is used to select the draws for  \code{\link{CmdStanVB}},
+and \code{\link{CmdStanLaplace}} fit objects.
+\item A type inheriting from \code{posterior::draws}. If the draws object has less
+samples than the number of requested chains/paths then the inits will be
+drawn using sampling with replacement. Otherwise sampling without
+replacement will be used. If the draws object's parameters are only a subset
+of the model parameters then the other parameters will be drawn by Stan's
+default initialization. The fit object must have at least some parameters
+that are the same name and dimensions as the current Stan model.
 }}
 
 \item{save_latent_dynamics}{(logical) Should auxiliary diagnostic information
diff --git a/man/model-method-pathfinder.Rd b/man/model-method-pathfinder.Rd
index d16ad170..1cda0502 100644
--- a/man/model-method-pathfinder.Rd
+++ b/man/model-method-pathfinder.Rd
@@ -85,6 +85,30 @@ take no arguments or a single argument \code{chain_id}. For MCMC, if the functio
 has argument \code{chain_id} it will be supplied with the chain id (from 1 to
 number of chains) when called to generate the initial values. See
 \strong{Examples}.
+\item A \code{\link{CmdStanMCMC}}, \code{\link{CmdStanMLE}}, \code{\link{CmdStanVB}}, \code{\link{CmdStanPathfinder}},
+or \code{\link{CmdStanLaplace}} fit object.
+If the fit object's parameters are only a subset of the model
+parameters then the other parameters will be drawn by Stan's default
+initialization. The fit object must have at least some parameters that are the
+same name and dimensions as the current Stan model. For the \code{sample} and
+\code{pathfinder} method, if the fit object has fewer draws than the requested
+number of chains/paths then the inits will be drawn using sampling with
+replacement. Otherwise sampling without replacement will be used.
+When a \code{\link{CmdStanPathfinder}} fit object is used as the init, if
+. \code{psis_resample} was set to \code{FALSE} and \code{calculate_lp} was
+set to \code{TRUE} (default), then resampling without replacement with Pareto
+smoothed weights will be used. If \code{psis_resample} was set to \code{TRUE} or
+\code{calculate_lp} was set to \code{FALSE} then sampling without replacement with
+uniform weights will be used to select the draws.
+PSIS resampling is used to select the draws for  \code{\link{CmdStanVB}},
+and \code{\link{CmdStanLaplace}} fit objects.
+\item A type inheriting from \code{posterior::draws}. If the draws object has less
+samples than the number of requested chains/paths then the inits will be
+drawn using sampling with replacement. Otherwise sampling without
+replacement will be used. If the draws object's parameters are only a subset
+of the model parameters then the other parameters will be drawn by Stan's
+default initialization. The fit object must have at least some parameters
+that are the same name and dimensions as the current Stan model.
 }}
 
 \item{save_latent_dynamics}{(logical) Should auxiliary diagnostic information
diff --git a/man/model-method-sample.Rd b/man/model-method-sample.Rd
index 51b4ce45..f6d105d1 100644
--- a/man/model-method-sample.Rd
+++ b/man/model-method-sample.Rd
@@ -98,6 +98,30 @@ take no arguments or a single argument \code{chain_id}. For MCMC, if the functio
 has argument \code{chain_id} it will be supplied with the chain id (from 1 to
 number of chains) when called to generate the initial values. See
 \strong{Examples}.
+\item A \code{\link{CmdStanMCMC}}, \code{\link{CmdStanMLE}}, \code{\link{CmdStanVB}}, \code{\link{CmdStanPathfinder}},
+or \code{\link{CmdStanLaplace}} fit object.
+If the fit object's parameters are only a subset of the model
+parameters then the other parameters will be drawn by Stan's default
+initialization. The fit object must have at least some parameters that are the
+same name and dimensions as the current Stan model. For the \code{sample} and
+\code{pathfinder} method, if the fit object has fewer draws than the requested
+number of chains/paths then the inits will be drawn using sampling with
+replacement. Otherwise sampling without replacement will be used.
+When a \code{\link{CmdStanPathfinder}} fit object is used as the init, if
+. \code{psis_resample} was set to \code{FALSE} and \code{calculate_lp} was
+set to \code{TRUE} (default), then resampling without replacement with Pareto
+smoothed weights will be used. If \code{psis_resample} was set to \code{TRUE} or
+\code{calculate_lp} was set to \code{FALSE} then sampling without replacement with
+uniform weights will be used to select the draws.
+PSIS resampling is used to select the draws for  \code{\link{CmdStanVB}},
+and \code{\link{CmdStanLaplace}} fit objects.
+\item A type inheriting from \code{posterior::draws}. If the draws object has less
+samples than the number of requested chains/paths then the inits will be
+drawn using sampling with replacement. Otherwise sampling without
+replacement will be used. If the draws object's parameters are only a subset
+of the model parameters then the other parameters will be drawn by Stan's
+default initialization. The fit object must have at least some parameters
+that are the same name and dimensions as the current Stan model.
 }}
 
 \item{save_latent_dynamics}{(logical) Should auxiliary diagnostic information
diff --git a/man/model-method-sample_mpi.Rd b/man/model-method-sample_mpi.Rd
index 078461bb..34083bd2 100644
--- a/man/model-method-sample_mpi.Rd
+++ b/man/model-method-sample_mpi.Rd
@@ -97,6 +97,30 @@ take no arguments or a single argument \code{chain_id}. For MCMC, if the functio
 has argument \code{chain_id} it will be supplied with the chain id (from 1 to
 number of chains) when called to generate the initial values. See
 \strong{Examples}.
+\item A \code{\link{CmdStanMCMC}}, \code{\link{CmdStanMLE}}, \code{\link{CmdStanVB}}, \code{\link{CmdStanPathfinder}},
+or \code{\link{CmdStanLaplace}} fit object.
+If the fit object's parameters are only a subset of the model
+parameters then the other parameters will be drawn by Stan's default
+initialization. The fit object must have at least some parameters that are the
+same name and dimensions as the current Stan model. For the \code{sample} and
+\code{pathfinder} method, if the fit object has fewer draws than the requested
+number of chains/paths then the inits will be drawn using sampling with
+replacement. Otherwise sampling without replacement will be used.
+When a \code{\link{CmdStanPathfinder}} fit object is used as the init, if
+. \code{psis_resample} was set to \code{FALSE} and \code{calculate_lp} was
+set to \code{TRUE} (default), then resampling without replacement with Pareto
+smoothed weights will be used. If \code{psis_resample} was set to \code{TRUE} or
+\code{calculate_lp} was set to \code{FALSE} then sampling without replacement with
+uniform weights will be used to select the draws.
+PSIS resampling is used to select the draws for  \code{\link{CmdStanVB}},
+and \code{\link{CmdStanLaplace}} fit objects.
+\item A type inheriting from \code{posterior::draws}. If the draws object has less
+samples than the number of requested chains/paths then the inits will be
+drawn using sampling with replacement. Otherwise sampling without
+replacement will be used. If the draws object's parameters are only a subset
+of the model parameters then the other parameters will be drawn by Stan's
+default initialization. The fit object must have at least some parameters
+that are the same name and dimensions as the current Stan model.
 }}
 
 \item{save_latent_dynamics}{(logical) Should auxiliary diagnostic information
diff --git a/man/model-method-variational.Rd b/man/model-method-variational.Rd
index 87decfd9..1d0355ba 100644
--- a/man/model-method-variational.Rd
+++ b/man/model-method-variational.Rd
@@ -81,6 +81,30 @@ take no arguments or a single argument \code{chain_id}. For MCMC, if the functio
 has argument \code{chain_id} it will be supplied with the chain id (from 1 to
 number of chains) when called to generate the initial values. See
 \strong{Examples}.
+\item A \code{\link{CmdStanMCMC}}, \code{\link{CmdStanMLE}}, \code{\link{CmdStanVB}}, \code{\link{CmdStanPathfinder}},
+or \code{\link{CmdStanLaplace}} fit object.
+If the fit object's parameters are only a subset of the model
+parameters then the other parameters will be drawn by Stan's default
+initialization. The fit object must have at least some parameters that are the
+same name and dimensions as the current Stan model. For the \code{sample} and
+\code{pathfinder} method, if the fit object has fewer draws than the requested
+number of chains/paths then the inits will be drawn using sampling with
+replacement. Otherwise sampling without replacement will be used.
+When a \code{\link{CmdStanPathfinder}} fit object is used as the init, if
+. \code{psis_resample} was set to \code{FALSE} and \code{calculate_lp} was
+set to \code{TRUE} (default), then resampling without replacement with Pareto
+smoothed weights will be used. If \code{psis_resample} was set to \code{TRUE} or
+\code{calculate_lp} was set to \code{FALSE} then sampling without replacement with
+uniform weights will be used to select the draws.
+PSIS resampling is used to select the draws for  \code{\link{CmdStanVB}},
+and \code{\link{CmdStanLaplace}} fit objects.
+\item A type inheriting from \code{posterior::draws}. If the draws object has less
+samples than the number of requested chains/paths then the inits will be
+drawn using sampling with replacement. Otherwise sampling without
+replacement will be used. If the draws object's parameters are only a subset
+of the model parameters then the other parameters will be drawn by Stan's
+default initialization. The fit object must have at least some parameters
+that are the same name and dimensions as the current Stan model.
 }}
 
 \item{save_latent_dynamics}{(logical) Should auxiliary diagnostic information
diff --git a/tests/testthat/resources/stan/logistic_simple.stan b/tests/testthat/resources/stan/logistic_simple.stan
new file mode 100644
index 00000000..95dec565
--- /dev/null
+++ b/tests/testthat/resources/stan/logistic_simple.stan
@@ -0,0 +1,11 @@
+data {
+  int<lower=0> N;
+  array[N] int<lower=0, upper=1> y;
+}
+parameters {
+  real alpha;
+}
+model {
+  target += normal_lpdf(alpha | 0, 1);
+  target += bernoulli_logit_lpmf(y | alpha);
+}
diff --git a/tests/testthat/resources/stan/parameter_types.stan b/tests/testthat/resources/stan/parameter_types.stan
new file mode 100644
index 00000000..c54fede6
--- /dev/null
+++ b/tests/testthat/resources/stan/parameter_types.stan
@@ -0,0 +1,15 @@
+parameters {
+  real real_p;
+  vector[2] vector_p;
+  matrix[2, 2] matrix_p;
+  array[2] matrix[2, 2] array_matrix_p;
+  corr_matrix[2] corr_p;
+  array[2, 2] real array_array_real_p;
+  array[2, 2] vector[3] array_array_vector_p;
+  array[2, 2] matrix[3, 3] array_array_matrix_p;
+//  complex complex_p;
+//  complex_matrix[2, 2] complex_matrix_p;
+//  complex_vector[4] complex_vector_p;
+//  tuple(real, vector[3], array[2] matrix[2, 2], complex) tuple_int_vector_arraymatrix_complex_p;
+//  array[2] tuple(real, tuple(vector[2], array[2] tuple(real, complex, matrix[2, 2]))) arraytuple_big_p;
+}
diff --git a/tests/testthat/test-fit-init.R b/tests/testthat/test-fit-init.R
new file mode 100644
index 00000000..509da198
--- /dev/null
+++ b/tests/testthat/test-fit-init.R
@@ -0,0 +1,130 @@
+context("fitted-inits")
+set_cmdstan_path()
+
+
+data_list_schools <- testing_data("schools")
+data_list_logistic <- testing_data("logistic")
+test_inits <- function(mod, fit_init, data_list = NULL) {
+  fit_sample <- mod$sample(data = data_list, chains = 1, init = fit_init,
+    iter_sampling = 100, iter_warmup = 100, refresh = 0, seed = 1234)
+  fit_sample_multi <- mod$sample(data = data_list, chains = 5, init = fit_init,
+    iter_sampling = 100, iter_warmup = 100, refresh = 0, seed = 1234)
+  fit_vb <- mod$variational(data = data_list, refresh = 0, seed = 1234,
+    init = fit_init, algorithm = "fullrank")
+  fit_path <- mod$pathfinder(data = data_list, seed=1234, refresh = 0,
+    num_paths = 4, init = fit_init)
+  fit_laplace <- mod$laplace(data = data_list, seed = 1234, refresh=0,
+    init=fit_init)
+  fit_ml <- mod$optimize(data = data_list, seed = 1234, refresh = 0,
+    init = fit_init, history_size = 400, algorithm = "lbfgs")
+  draws = posterior::as_draws_rvars(fit_init$draws())
+  fit_sample_draws <- mod$sample(data = data_list, chains = 1, init = draws,
+    iter_sampling = 100, iter_warmup = 100, refresh = 0, seed = 1234)
+  return(0)
+}
+
+test_that("Sample method works as init", {
+  mod_params <- testing_model("parameter_types")
+  utils::capture.output(fit_sample_init <- mod_params$sample(chains = 1,
+    iter_warmup = 100, iter_sampling = 100, refresh = 0, seed = 1234))
+  expect_no_error(test_inits(mod_params, fit_sample_init))
+})
+
+test_that("Multi chain Sample method works as init", {
+  mod_params <- testing_model("parameter_types")
+  utils::capture.output(fit_sample_multi_init <- mod_params$sample(chains = 4,
+    iter_warmup = 100, iter_sampling = 100, refresh = 0, seed = 1234))
+  expect_no_error(test_inits(mod_params, fit_sample_multi_init))
+})
+
+test_that("Subsets of parameters are allowed", {
+  mod_logistic_simple <- testing_model("logistic_simple")
+  utils::capture.output(fit_sample_init_simple <- mod_logistic_simple$sample(chains = 1,
+    data = data_list_logistic, iter_warmup = 100, iter_sampling = 100,
+    refresh = 0, seed = 1234))
+  mod_logistic <- testing_model("logistic")
+  expect_no_error(test_inits(mod_logistic, fit_sample_init_simple,
+    data_list_logistic))
+})
+
+test_that("Pathfinder works as init", {
+  mod_logistic <- testing_model("logistic")
+  utils::capture.output(fit_path_init <- mod_logistic$pathfinder(
+    seed=1234, data = data_list_logistic, refresh = 0, num_paths = 1))
+  expect_no_error(test_inits(mod_logistic, fit_path_init, data_list_logistic))
+})
+
+test_that("Multi Pathfinder method works as init", {
+  mod_logistic <- testing_model("logistic")
+  utils::capture.output(fit_path_init <- mod_logistic$pathfinder(seed=1234,
+    data = data_list_logistic, refresh = 0, num_paths = 4))
+  expect_no_error(test_inits(mod_logistic, fit_path_init, data_list_logistic))
+})
+
+test_that("Pathfinder method with psis_resample as false works as init", {
+  mod_logistic <- testing_model("logistic")
+  utils::capture.output(fit_path_init <- mod_logistic$pathfinder(
+    seed=1234, data = data_list_logistic, refresh = 0, num_paths = 1))
+  expect_no_error(test_inits(mod_logistic, fit_path_init, data_list_logistic))
+})
+
+
+test_that("Multi Pathfinder method with psis_resample as false works as init", {
+  mod_logistic <- testing_model("logistic")
+  utils::capture.output(fit_path_init <- mod_logistic$pathfinder(
+    seed=1234, data = data_list_logistic, refresh = 0, num_paths = 4,
+    psis_resample = FALSE))
+  expect_no_error(test_inits(mod_logistic, fit_path_init, data_list_logistic))
+})
+
+
+test_that("Pathfinder method with calculate_lp as false works as init", {
+  mod_logistic <- testing_model("logistic")
+  fit_path_init <- mod_logistic$pathfinder(
+    seed=1234, data = data_list_logistic, refresh = 0, num_paths = 1,
+    psis_resample = FALSE, calculate_lp = FALSE)
+  expect_no_error(test_inits(mod_logistic, fit_path_init, data_list_logistic))
+})
+
+test_that("Multi Pathfinder method with calculate_lp as false works as init", {
+  mod_logistic <- testing_model("logistic")
+  utils::capture.output(fit_path_init <- mod_logistic$pathfinder(
+    seed=1234, data = data_list_logistic, refresh = 0, num_paths = 4,
+    psis_resample = TRUE, calculate_lp = FALSE))
+  expect_no_error(test_inits(mod_logistic, fit_path_init, data_list_logistic))
+})
+
+test_that("Variational method works as init", {
+  mod_logistic <- testing_model("logistic")
+  utils::capture.output(fit_vb_init <- mod_logistic$variational(
+    data = data_list_logistic, seed=1234, refresh = 0))
+  expect_no_error(test_inits(mod_logistic, fit_vb_init, data_list_logistic))
+})
+
+test_that("Optimization method works as init", {
+  mod_logistic <- testing_model("logistic")
+  utils::capture.output(fit_ml_init <- mod_logistic$optimize(
+    data = data_list_logistic, seed=1234, refresh = 0))
+  expect_no_error(test_inits(mod_logistic, fit_ml_init, data_list_logistic))
+})
+
+
+test_that("Draws Object with NA or Inf throws error", {
+  mod_logistic <- testing_model("logistic")
+  utils::capture.output(fit_laplace_init <- mod_logistic$laplace(
+    data = data_list_logistic, seed = 1234, refresh=0))
+  draws_df = fit_laplace_init$draws()
+  draws_df[1, 3] = NA
+  expect_error(mod_logistic$laplace(
+    data = data_list_logistic, seed = 1234, refresh=0, init = draws_df[1, ]), "alpha contains NA or Inf values!")
+  draws_df[1, 4] = NA
+  expect_error(mod_logistic$sample(
+    data = data_list_logistic, seed = 1234, refresh=0, init = draws_df[1:4, ]), "alpha, beta contains NA or Inf values!")
+  draws_df = fit_laplace_init$draws()
+  draws_df[1, 3] = Inf
+  expect_error(mod_logistic$sample(
+    data = data_list_logistic, seed = 1234, refresh=0, init = draws_df[1:4, ]), "alpha contains NA or Inf values!")
+  draws_df[1, 4] = NA
+  expect_error(mod_logistic$sample(
+    data = data_list_logistic, seed = 1234, refresh=0, init = draws_df[1:4, ]), "alpha, beta contains NA or Inf values!")
+})