bebi103.hmc.inferential_calibration_cmdstanpy

bebi103.hmc.inferential_calibration_cmdstanpy(prior_predictive_fun=None, prior_predictive_kwargs={}, prior_sds_fun=None, prior_sds_kwargs={}, posterior_stan_file=None, measured_data=None, additional_posterior_model_data={}, var_names=None, ignored_var_names=(), sample_kwargs=None, diagnostic_check_kwargs=None, cores=1, N=400, samples_dir='inferential_calibration_samples', remove_sample_files=True, df_package='polars', progress_bar=False, output_csv='inferential_calibration_results.csv', overwrite=False)

Perform simulation-based calibration for models using the cmdstanpy interface.

Parameters:

  • prior_predictive_fun (function) – Function with call signature prior_predictive_fun(seed=None, **prior_predictive_kwargs). This function is used to draw prior predictive samples of data sets from the posterior. The function need not necessarily use a provided seed, but must have a keyword argument allowing for seed specification. The function must still run, even if seed is None. The function must return a dictionary of a single prior predictive sample. Each key is the variable name and each value is either a parameter drawn from the prior or a data set drawn from the prior predictive distribution. Both variables and data sets may be entered in the dictionary as Numpy arrays/scalars or as individual elements of Numpy arrays (e.g., with keys ‘my_var[1]’, ‘my_var[2]’, etc.)

  • prior_predictive_kwargs (dict, default {}) – Kwargs to be passed to prior_predictive_fun().

  • prior_sds_fun (function, default None) – Function with call signature prior_sds_fun(**prior_sds_kwargs). This function is used to compute the standard deviations of the marginal prior for each parameter. The function must return a dictionary where each key is the variable name and each value is the standard deviation of the corresponding marginal prior distribution. Note that any parameter must be a scalar. E.g., for a matrix-valued parameter A, the keys must be ‘A[1,1]’, ‘A[1,2]’, and so on. Note the 1-indexing. If None, then prior_predictive_fun() is used to draw samples out of the priors and then standard deviations are estimated from the samples.

  • prior_sds_kwargs (dict, default {}) – Kwargs to be passed to prior_sds_fun(). Ignored if prior_sds_fun is None.

  • posterior_stan_file (str) – Stan file for the posterior that allows sampling.

  • measured_data (dict) – A dictionary where each key is the variable name of measured data. The value of each entry in the measured_data dictionary is the type of the corresponding data, almost always either int or float. Upon generation of a prior predictive data set, a new dictionary is created with keys given by the keys of measured_data and value given by the generated prior predictive data set.

  • additional_posterior_model_data (dict, default {}) – Dictionary with entries specified for the data block of the posterior Stan model. These are all entries in the data block of the posterior model that are not measured data generated by prior_predictive_fun().

  • var_names (None, or tuple or list of strings) – A list of strings containing parameter names to be considered in the inferential calibration analysis. Not all parameters of the model need be considered; only those in var_names have rank statistics calculated. Note that for multidimensional variables, var_names only has the root name. E.g., var_names=[‘x’, ‘y’], not something like var_names=[‘x[0]’, ‘x[1]’, ‘y[0,0]’]. If None, names of variables are inferred from output of prior_predictive_fun().

  • ignored_var_names (tuple or list of strings, default ()) – Variable names to ignore in the inferential calibration procedure. These variables have no rank statistics computed and are omitted from the diagnostic checks. This is typically used for posterior retrodictive variables and log likelihood variables. As with var_names, only the root name is given for multidimensional variables, e.g. ignored_var_names=[‘x_prc’], not ignored_var_names=[‘x_prc[0]’, ‘x_prc[1]’].

  • sample_kwargs (dict, default None) – kwargs to be passed to sm.sample().

  • diagnostic_check_kwargs (dict, default None) – kwargs to pass to check_all_diagnostics(). If verbosity and/or return_diagnostics are given, they are ignored. max_treedepth is inferred from sample_kwargs.

  • cores (int, default 1) – Number of cores to use in the inferential calibration calculation.

  • N (int, 400) – Number of simulations to run.

  • samples_dir (str, default "inferential_calibration_samples") – Path to directory to store .csv and .txt files generated by CmdStan. The directory specified here will NOT be destroyed after a calculation, though the files within it may be, depending on the remove_sample_files kwarg.

  • remove_sample_files (bool, default True) – If True, remove .csv and .txt files generated by CmdStan.

  • df_package (str, either 'polars' (default) or 'pandas') – Which package to use for output data frame

  • progress_bar (bool, default False) – If True, display a progress bar for the calculation using tqdm.

  • output_csv (str, default 'inferential_calibration_results.csv') – Path to CSV file where results will be written. Results are written directly to this file as workers complete (instead of collecting in memory). This reduces memory usage and allows monitoring progress by examining the growing CSV file. The file includes metadata comments at the top with run parameters and start time.

  • overwrite (bool, default False) – If False, raises RuntimeError if output_csv already exists. If True, overwrites existing file with new results.

Returns:

output – Data frame with the output of the inferential calibration analysis. It has one row per scalar parameter per trial and the following columns. If a trial fails, its per-parameter statistics and its HMC diagnostic summaries are null (Polars output) or np.nan (Pandas output), and the error column holds the error message.

  • trial : Unique trial number for the simulation.

  • parameter : The name of the scalar parameter.

  • ground_truth : Value of the parameter used in the simulation. This value was drawn out of the prior distribution.

  • rank_statistic : Value of the rank statistic for the parameter for the trial.

  • mean : mean parameter value based on sampling out of the posterior in the simulation.

  • sd : standard deviation of the parameter value based on sampling out of the posterior in the simulation.

  • shrinkage : The shrinkage for the parameter for the given trial. This is computed as 1 - (sd / sd_prior)**2, where sd_prior is the standard deviation of the parameter’s marginal prior distribution.

  • z_score : The z-score for the parameter for the given trial. This is computed as (mean - ground_truth) / sd, so it is signed.

  • Rhat : Split Rhat for the parameter, computed across all chains.

  • ESS_hat : Total empirical effective sample size for the parameter, summed over all chains.

  • ESS_hat_per_iter : ESS_hat divided by the total number of posterior samples (all chains and draws).

  • ESS_hat_min_per_chain : Minimum per-chain empirical effective sample size for the parameter. This is the quantity the ESS_hat diagnostic warning is based on: a warning is raised (and the corresponding bit set in warning_code) when it falls below min_ess_hat_per_chain. Because ESS_hat is a sum over chains, it can be large even when this per-chain minimum, and hence the warning, indicates poor mixing in a single chain.

  • tail_xi_hat_left : Left tail xi_hat for the parameter, computed over all chains pooled together.

  • tail_xi_hat_right : Right tail xi_hat for the parameter, computed over all chains pooled together.

  • inc_tau_hat : Incremental empirical integrated autocorrelation time for the parameter, computed over all chains pooled together.

  • variance : Empirical variance of the parameter’s posterior samples, computed over all chains pooled together.

  • n_divergences : Total number of divergent transitions in the trial, summed over all chains.

  • n_bad_ebfmi : Number of chains in the trial whose E-BFMI is below efmi_rule_of_thumb.

  • n_max_treedepth : Total number of transitions in the trial that hit the maximum tree depth, summed over all chains.

  • average_accept_proxy : Average proxy acceptance statistic for the trial, taken over all chains.

  • warning_code : Warning code based on diagnostic checks outputted by check_all_diagnostics().

  • prior_predictive_seed : Random number seed passed to prior_predictive_fun() to generate the prior predictive sample for the trial.

  • posterior_seed : Random number seed used by Stan for the posterior MCMC sampling for the trial.

  • L : The number of posterior draws used in computing the rank statistic. The rank statistic should be uniform on the integers [0, L].

  • samples_prefix : Full-path prefix of the sampling files (.csv and .txt) generated by Stan for the posterior MCMC sampling for the trial. If remove_sample_files is True (so the files are deleted), every entry is null (Polars output) or np.nan (Pandas output).

  • error : ‘no error’ if the trial succeeded; otherwise the error message raised during the trial.

Return type:

Polars or Pandas DataFrame

Notes

Each simulation is done by sampling a parameter set out of the prior distribution, using those parameters to generate data from the likelihood, and then performing posterior sampling based on the generated data. A rank statistic for each simulation is computed. This rank statistic should be uniformly distributed over its L possible values. See https://arxiv.org/abs/1804.06788, by Talts, et al., for details.

Prior predictive samples are not stored. The assumption is that prior predictive samples can be drawn very quickly. The seed used to generated the prior predictive sample and other kwargs passed to prior_predictive_fun are stored for each trial, so it is easy to regenerated prior predictive samples as desired.