sbi: simulation-based inference¶
sbi: A Python toolbox for simulation-based inference.
Inference can be run in a single line of code:
+Inference can be run in a single line of code
posterior = infer(simulator, prior, method='SNPE', num_simulations=1000)
and you can choose from a variety of amortized and sequential SBI methods.
+or in a few lines for more flexibility:
+inference = SNPE(prior=prior)
+_ = inference.append_simulations(theta, x).train()
+posterior = inference.build_posterior()
+sbi lets you choose from a variety of amortized and sequential SBI methods:
Amortized methods return a posterior that can be applied to many different observations without retraining, whereas sequential methods focus the inference on one particular observation to be more simulation-efficient. For an overview of implemented methods see below, or checkout or GitHub page.
diff --git a/install/index.html b/install/index.html index 07a4b9903..d779b7f25 100644 --- a/install/index.html +++ b/install/index.html @@ -319,20 +319,6 @@ -- the simulator can simulate batches of parameters and return batches of data.
- the prior does not produce batches and samples and evaluates to Tensor.
- the output shape is a `torch.Size((1,N))` (i.e, has a leading batch dimension 1). If this is not possible, a suitable exception will be raised. Args: simulator: Simulator as provided by the user. prior: Prior as provided by the user. Returns: Tuple (simulator, prior) checked and matching the requirements of sbi. \"\"\" # Check prior, return PyTorch prior. prior , _ , prior_returns_numpy = process_prior ( prior ) # Check simulator, returns PyTorch simulator able to simulate batches. simulator = process_simulator ( simulator , prior , prior_returns_numpy ) # Consistency check after making ready for sbi. check_sbi_inputs ( simulator , prior ) return simulator , prior sbi . inference . base . simulate_for_sbi ( simulator , proposal , num_simulations , num_workers = 1 , simulation_batch_size = 1 , seed = None , show_progress_bar = True ) \u00b6 Returns ( \\(\\theta, x\\) ) pairs obtained from sampling the proposal and simulating. This function performs two steps: Sample parameters \\(\\theta\\) from the proposal . Simulate these parameters to obtain \\(x\\) . Parameters: Name Type Description Default simulator Callable A function that takes parameters \\(\\theta\\) and maps them to simulations, or observations, x , \\(\\text{sim}(\\theta)\\to x\\) . Any regular Python callable (i.e. function or class with __call__ method) can be used. required proposal Any Probability distribution that the parameters \\(\\theta\\) are sampled from. required num_simulations int Number of simulations that are run. required num_workers int Number of parallel workers to use for simulations. 1 simulation_batch_size int Number of parameter sets that the simulator maps to data x at once. If None, we simulate all parameter sets at the same time. If >= 1, the simulator has to process data of shape (simulation_batch_size, parameter_dimension). 1 seed Optional[int] Seed for reproducibility. None show_progress_bar bool Whether to show a progress bar for simulating. This will not affect whether there will be a progressbar while drawing samples from the proposal. True Returns: Sampled parameters \\(\\theta\\) and simulation-outputs \\(x\\) . Source code in sbi/inference/base.py def simulate_for_sbi ( simulator : Callable , proposal : Any , num_simulations : int , num_workers : int = 1 , simulation_batch_size : int = 1 , seed : Optional [ int ] = None , show_progress_bar : bool = True , ) -> Tuple [ Tensor , Tensor ]: r \"\"\"Returns ($\\theta, x$) pairs obtained from sampling the proposal and simulating. This function performs two steps: - Sample parameters $\\theta$ from the `proposal`. - Simulate these parameters to obtain $x$. Args: simulator: A function that takes parameters $\\theta$ and maps them to simulations, or observations, `x`, $\\text{sim}(\\theta)\\to x$. Any regular Python callable (i.e. function or class with `__call__` method) can be used. proposal: Probability distribution that the parameters $\\theta$ are sampled from. num_simulations: Number of simulations that are run. num_workers: Number of parallel workers to use for simulations. simulation_batch_size: Number of parameter sets that the simulator maps to data x at once. If None, we simulate all parameter sets at the same time. If >= 1, the simulator has to process data of shape (simulation_batch_size, parameter_dimension). seed: Seed for reproducibility. show_progress_bar: Whether to show a progress bar for simulating. This will not affect whether there will be a progressbar while drawing samples from the proposal. Returns: Sampled parameters $\\theta$ and simulation-outputs $x$. \"\"\" theta = proposal . sample (( num_simulations ,)) x = simulate_in_batches ( simulator = simulator , theta = theta , sim_batch_size = simulation_batch_size , num_workers = num_workers , seed = seed , show_progress_bars = show_progress_bar , ) return theta , x sbi.inference.snpe.snpe_a.SNPE_A ( PosteriorEstimator ) \u00b6 __init__ ( self , prior = None , density_estimator = 'mdn_snpe_a' , num_components = 10 , device = 'cpu' , logging_level = 'WARNING' , summary_writer = None , show_progress_bars = True ) special \u00b6 SNPE-A [1]. [1] Fast epsilon-free Inference of Simulation Models with Bayesian Conditional Density Estimation , Papamakarios et al., NeurIPS 2016, https://arxiv.org/abs/1605.06376 . This class implements SNPE-A. SNPE-A trains across multiple rounds with a maximum-likelihood-loss. This will make training converge to the proposal posterior instead of the true posterior. To correct for this, SNPE-A applies a post-hoc correction after training. This correction has to be performed analytically. Thus, SNPE-A is limited to Gaussian distributions for all but the last round. In the last round, SNPE-A can use a Mixture of Gaussians. Parameters: Name Type Description Default prior Optional[torch.distributions.distribution.Distribution] A probability distribution that expresses prior knowledge about the parameters, e.g. which ranges are meaningful for them. Any object with .log_prob() and .sample() (for example, a PyTorch distribution) can be used. None density_estimator Union[str, Callable] If it is a string (only \u201cmdn_snpe_a\u201d is valid), use a pre-configured mixture of densities network. Alternatively, a function that builds a custom neural network can be provided. The function will be called with the first batch of simulations (theta, x), which can thus be used for shape inference and potentially for z-scoring. It needs to return a PyTorch nn.Module implementing the density estimator. The density estimator needs to provide the methods .log_prob and .sample() . Note that until the last round only a single (multivariate) Gaussian component is used for training (see Algorithm 1 in [1]). In the last round, this component is replicated num_components times, its parameters are perturbed with a very small noise, and then the last training round is done with the expanded Gaussian mixture as estimator for the proposal posterior. 'mdn_snpe_a' num_components int Number of components of the mixture of Gaussians in the last round. This overrides the num_components value passed to posterior_nn() . 10 device str Training device, e.g., \u201ccpu\u201d, \u201ccuda\u201d or \u201ccuda:{0, 1, \u2026}\u201d. 'cpu' logging_level Union[int, str] Minimum severity of messages to log. One of the strings INFO, WARNING, DEBUG, ERROR and CRITICAL. 'WARNING' summary_writer Optional[Writer] A tensorboard SummaryWriter to control, among others, log file location (default is
- proposal is a `DirectPosterior` with density_estimator `mdn`, as built with `utils.sbi.posterior_nn()`.
- the density estimator is a `mdn`, as built with `utils.sbi.posterior_nn()`.
- `isinstance(prior, MultivariateNormal)` (from `torch.distributions`) or `isinstance(prior, sbi.utils.BoxUniform)` Note that custom implementations of any of these densities (or estimators) will not trigger the non-atomic loss, and the algorithm will fall back onto using the atomic loss. Args: prior: A probability distribution that expresses prior knowledge about the parameters, e.g. which ranges are meaningful for them. density_estimator: If it is a string, use a pre-configured network of the provided type (one of nsf, maf, mdn, made). Alternatively, a function that builds a custom neural network can be provided. The function will be called with the first batch of simulations (theta, x), which can thus be used for shape inference and potentially for z-scoring. It needs to return a PyTorch `nn.Module` implementing the density estimator. The density estimator needs to provide the methods `.log_prob` and `.sample()`. device: Training device, e.g., \"cpu\", \"cuda\" or \"cuda:{0, 1, ...}\". logging_level: Minimum severity of messages to log. One of the strings INFO, WARNING, DEBUG, ERROR and CRITICAL. summary_writer: A tensorboard `SummaryWriter` to control, among others, log file location (default is `
- the simulator can simulate batches of parameters and return batches of data.
- the prior does not produce batches and samples and evaluates to Tensor.
- the output shape is a `torch.Size((1,N))` (i.e, has a leading batch dimension 1). If this is not possible, a suitable exception will be raised. Args: simulator: Simulator as provided by the user. prior: Prior as provided by the user. Returns: Tuple (simulator, prior) checked and matching the requirements of sbi. \"\"\" # Check prior, return PyTorch prior. prior , _ , prior_returns_numpy = process_prior ( prior ) # Check simulator, returns PyTorch simulator able to simulate batches. simulator = process_simulator ( simulator , prior , prior_returns_numpy ) # Consistency check after making ready for sbi. check_sbi_inputs ( simulator , prior ) return simulator , prior","title":"prepare_for_sbi()"},{"location":"reference/#sbi.inference.base.simulate_for_sbi","text":"Returns ( \\(\\theta, x\\) ) pairs obtained from sampling the proposal and simulating. This function performs two steps: Sample parameters \\(\\theta\\) from the proposal . Simulate these parameters to obtain \\(x\\) . Parameters: Name Type Description Default simulator Callable A function that takes parameters \\(\\theta\\) and maps them to simulations, or observations, x , \\(\\text{sim}(\\theta)\\to x\\) . Any regular Python callable (i.e. function or class with __call__ method) can be used. required proposal Any Probability distribution that the parameters \\(\\theta\\) are sampled from. required num_simulations int Number of simulations that are run. required num_workers int Number of parallel workers to use for simulations. 1 simulation_batch_size int Number of parameter sets that the simulator maps to data x at once. If None, we simulate all parameter sets at the same time. If >= 1, the simulator has to process data of shape (simulation_batch_size, parameter_dimension). 1 seed Optional[int] Seed for reproducibility. None show_progress_bar bool Whether to show a progress bar for simulating. This will not affect whether there will be a progressbar while drawing samples from the proposal. True Returns: Sampled parameters \\(\\theta\\) and simulation-outputs \\(x\\) . Source code in sbi/inference/base.py def simulate_for_sbi ( simulator : Callable , proposal : Any , num_simulations : int , num_workers : int = 1 , simulation_batch_size : int = 1 , seed : Optional [ int ] = None , show_progress_bar : bool = True , ) -> Tuple [ Tensor , Tensor ]: r \"\"\"Returns ($\\theta, x$) pairs obtained from sampling the proposal and simulating. This function performs two steps: - Sample parameters $\\theta$ from the `proposal`. - Simulate these parameters to obtain $x$. Args: simulator: A function that takes parameters $\\theta$ and maps them to simulations, or observations, `x`, $\\text{sim}(\\theta)\\to x$. Any regular Python callable (i.e. function or class with `__call__` method) can be used. proposal: Probability distribution that the parameters $\\theta$ are sampled from. num_simulations: Number of simulations that are run. num_workers: Number of parallel workers to use for simulations. simulation_batch_size: Number of parameter sets that the simulator maps to data x at once. If None, we simulate all parameter sets at the same time. If >= 1, the simulator has to process data of shape (simulation_batch_size, parameter_dimension). seed: Seed for reproducibility. show_progress_bar: Whether to show a progress bar for simulating. This will not affect whether there will be a progressbar while drawing samples from the proposal. Returns: Sampled parameters $\\theta$ and simulation-outputs $x$. \"\"\" theta = proposal . sample (( num_simulations ,)) x = simulate_in_batches ( simulator = simulator , theta = theta , sim_batch_size = simulation_batch_size , num_workers = num_workers , seed = seed , show_progress_bars = show_progress_bar , ) return theta , x","title":"simulate_for_sbi()"},{"location":"reference/#sbi.inference.snpe.snpe_a.SNPE_A","text":"","title":"SNPE_A"},{"location":"reference/#sbi.inference.snpe.snpe_c.SNPE_C","text":"","title":"SNPE_C"},{"location":"reference/#sbi.inference.snle.snle_a.SNLE_A","text":"","title":"SNLE_A"},{"location":"reference/#sbi.inference.snre.snre_a.SNRE_A","text":"","title":"SNRE_A"},{"location":"reference/#sbi.inference.snre.snre_b.SNRE_B","text":"","title":"SNRE_B"},{"location":"reference/#sbi.inference.snre.snre_c.SNRE_C","text":"","title":"SNRE_C"},{"location":"reference/#sbi.inference.snre.bnre.BNRE","text":"","title":"BNRE"},{"location":"reference/#sbi.inference.abc.mcabc.MCABC","text":"","title":"MCABC"},{"location":"reference/#sbi.inference.abc.smcabc.SMCABC","text":"","title":"SMCABC"},{"location":"reference/#posteriors","text":"","title":"Posteriors"},{"location":"reference/#sbi.inference.posteriors.direct_posterior.DirectPosterior","text":"Posterior \\(p(\\theta|x_o)\\) with log_prob() and sample() methods, only applicable to SNPE. SNPE trains a neural network to directly approximate the posterior distribution. However, for bounded priors, the neural network can have leakage: it puts non-zero mass in regions where the prior is zero. The DirectPosterior class wraps the trained network to deal with these cases. Specifically, this class offers the following functionality: - correct the calculation of the log probability such that it compensates for the leakage. - reject samples that lie outside of the prior bounds. This class can not be used in combination with SNLE or SNRE.","title":"DirectPosterior"},{"location":"reference/#sbi.inference.posteriors.importance_posterior.ImportanceSamplingPosterior","text":"Provides importance sampling to sample from the posterior. SNLE or SNRE train neural networks to approximate the likelihood(-ratios). ImportanceSamplingPosterior allows to estimate the posterior log-probability by estimating the normlalization constant with importance sampling. It also allows to perform importance sampling (with .sample() ) and to draw approximate samples with sampling-importance-resampling (SIR) (with .sir_sample() )","title":"ImportanceSamplingPosterior"},{"location":"reference/#sbi.inference.posteriors.mcmc_posterior.MCMCPosterior","text":"Provides MCMC to sample from the posterior. SNLE or SNRE train neural networks to approximate the likelihood(-ratios). MCMCPosterior allows to sample from the posterior with MCMC.","title":"MCMCPosterior"},{"location":"reference/#sbi.inference.posteriors.rejection_posterior.RejectionPosterior","text":"Provides rejection sampling to sample from the posterior. SNLE or SNRE train neural networks to approximate the likelihood(-ratios). RejectionPosterior allows to sample from the posterior with rejection sampling.","title":"RejectionPosterior"},{"location":"reference/#sbi.inference.posteriors.vi_posterior.VIPosterior","text":"Provides VI (Variational Inference) to sample from the posterior. SNLE or SNRE train neural networks to approximate the likelihood(-ratios). VIPosterior allows to learn a tractable variational posterior \\(q(\\theta)\\) which approximates the true posterior \\(p(\\theta|x_o)\\) . After this second training stage, we can produce approximate posterior samples, by just sampling from q with no additional cost. For additional information see [1] and [2]. References: [1] Variational methods for simulation-based inference, Manuel Gl\u00f6ckler, Michael Deistler, Jakob Macke, 2022, https://openreview.net/forum?id=kZ0UYdhqkNY [2] Sequential Neural Posterior and Likelihood Approximation, Samuel Wiqvist, Jes Frellsen, Umberto Picchini, 2021, https://arxiv.org/abs/2102.06522","title":"VIPosterior"},{"location":"reference/#models","text":"","title":"Models"},{"location":"reference/#sbi.utils.get_nn_models.posterior_nn","text":"Returns a function that builds a density estimator for learning the posterior. This function will usually be used for SNPE. The returned function is to be passed to the inference class when using the flexible interface. Parameters: Name Type Description Default model str The type of density estimator that will be created. One of [ mdn , made , maf , maf_rqs , nsf ]. required z_score_theta Optional[str] Whether to z-score parameters \\(\\theta\\) before passing them into the network, can take one of the following: - none , or None: do not z-score. - independent : z-score each dimension independently. - structured : treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. 'independent' z_score_x Optional[str] Whether to z-score simulation outputs \\(x\\) before passing them into the network, same options as z_score_theta. 'independent' hidden_features int Number of hidden features. 50 num_transforms int Number of transforms when a flow is used. Only relevant if density estimator is a normalizing flow (i.e. currently either a maf or a nsf ). Ignored if density estimator is a mdn or made . 5 num_bins int Number of bins used for the splines in nsf . Ignored if density estimator not nsf . 10 embedding_net Module Optional embedding network for simulation outputs \\(x\\) . This embedding net allows to learn features from potentially high-dimensional simulation outputs. Identity() num_components int Number of mixture components for a mixture of Gaussians. Ignored if density estimator is not an mdn. 10 kwargs additional custom arguments passed to downstream build functions. {} Source code in sbi/utils/get_nn_models.py def posterior_nn ( model : str , z_score_theta : Optional [ str ] = \"independent\" , z_score_x : Optional [ str ] = \"independent\" , hidden_features : int = 50 , num_transforms : int = 5 , num_bins : int = 10 , embedding_net : nn . Module = nn . Identity (), num_components : int = 10 , ** kwargs , ) -> Callable : r \"\"\" Returns a function that builds a density estimator for learning the posterior. This function will usually be used for SNPE. The returned function is to be passed to the inference class when using the flexible interface. Args: model: The type of density estimator that will be created. One of [`mdn`, `made`, `maf`, `maf_rqs`, `nsf`]. z_score_theta: Whether to z-score parameters $\\theta$ before passing them into the network, can take one of the following: - `none`, or None: do not z-score. - `independent`: z-score each dimension independently. - `structured`: treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. z_score_x: Whether to z-score simulation outputs $x$ before passing them into the network, same options as z_score_theta. hidden_features: Number of hidden features. num_transforms: Number of transforms when a flow is used. Only relevant if density estimator is a normalizing flow (i.e. currently either a `maf` or a `nsf`). Ignored if density estimator is a `mdn` or `made`. num_bins: Number of bins used for the splines in `nsf`. Ignored if density estimator not `nsf`. embedding_net: Optional embedding network for simulation outputs $x$. This embedding net allows to learn features from potentially high-dimensional simulation outputs. num_components: Number of mixture components for a mixture of Gaussians. Ignored if density estimator is not an mdn. kwargs: additional custom arguments passed to downstream build functions. \"\"\" kwargs = dict ( zip ( ( \"z_score_x\" , \"z_score_y\" , \"hidden_features\" , \"num_transforms\" , \"num_bins\" , \"embedding_net\" , \"num_components\" , ), ( z_score_theta , z_score_x , hidden_features , num_transforms , num_bins , embedding_net , num_components , ), ), ** kwargs , ) def build_fn_snpe_a ( batch_theta , batch_x , num_components ): \"\"\"Build function for SNPE-A Extract the number of components from the kwargs, such that they are exposed as a kwargs, offering the possibility to later override this kwarg with `functools.partial`. This is necessary in order to make sure that the MDN in SNPE-A only has one component when running the Algorithm 1 part. \"\"\" return build_mdn ( batch_x = batch_theta , batch_y = batch_x , num_components = num_components , ** kwargs , ) def build_fn ( batch_theta , batch_x ): if model == \"mdn\" : return build_mdn ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) elif model == \"made\" : return build_made ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) elif model == \"maf\" : return build_maf ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) elif model == \"maf_rqs\" : return build_maf_rqs ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) elif model == \"nsf\" : return build_nsf ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) else : raise NotImplementedError if model == \"mdn_snpe_a\" : if num_components != 10 : raise ValueError ( \"You set `num_components`. For SNPE-A, this has to be done at \" \"instantiation of the inference object, i.e. \" \"`inference = SNPE_A(..., num_components=20)`\" ) kwargs . pop ( \"num_components\" ) return build_fn_snpe_a if model == \"mdn_snpe_a\" else build_fn","title":"posterior_nn()"},{"location":"reference/#sbi.utils.get_nn_models.likelihood_nn","text":"Returns a function that builds a density estimator for learning the likelihood. This function will usually be used for SNLE. The returned function is to be passed to the inference class when using the flexible interface. Parameters: Name Type Description Default model str The type of density estimator that will be created. One of [ mdn , made , maf , maf_rqs , nsf ]. required z_score_theta Optional[str] Whether to z-score parameters \\(\\theta\\) before passing them into the network, can take one of the following: - none , or None: do not z-score. - independent : z-score each dimension independently. - structured : treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. 'independent' z_score_x Optional[str] Whether to z-score simulation outputs \\(x\\) before passing them into the network, same options as z_score_theta. 'independent' hidden_features int Number of hidden features. 50 num_transforms int Number of transforms when a flow is used. Only relevant if density estimator is a normalizing flow (i.e. currently either a maf or a nsf ). Ignored if density estimator is a mdn or made . 5 num_bins int Number of bins used for the splines in nsf . Ignored if density estimator not nsf . 10 embedding_net Module Optional embedding network for parameters \\(\\theta\\) . Identity() num_components int Number of mixture components for a mixture of Gaussians. Ignored if density estimator is not an mdn. 10 kwargs additional custom arguments passed to downstream build functions. {} Source code in sbi/utils/get_nn_models.py def likelihood_nn ( model : str , z_score_theta : Optional [ str ] = \"independent\" , z_score_x : Optional [ str ] = \"independent\" , hidden_features : int = 50 , num_transforms : int = 5 , num_bins : int = 10 , embedding_net : nn . Module = nn . Identity (), num_components : int = 10 , ** kwargs , ) -> Callable : r \"\"\" Returns a function that builds a density estimator for learning the likelihood. This function will usually be used for SNLE. The returned function is to be passed to the inference class when using the flexible interface. Args: model: The type of density estimator that will be created. One of [`mdn`, `made`, `maf`, `maf_rqs`, `nsf`]. z_score_theta: Whether to z-score parameters $\\theta$ before passing them into the network, can take one of the following: - `none`, or None: do not z-score. - `independent`: z-score each dimension independently. - `structured`: treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. z_score_x: Whether to z-score simulation outputs $x$ before passing them into the network, same options as z_score_theta. hidden_features: Number of hidden features. num_transforms: Number of transforms when a flow is used. Only relevant if density estimator is a normalizing flow (i.e. currently either a `maf` or a `nsf`). Ignored if density estimator is a `mdn` or `made`. num_bins: Number of bins used for the splines in `nsf`. Ignored if density estimator not `nsf`. embedding_net: Optional embedding network for parameters $\\theta$. num_components: Number of mixture components for a mixture of Gaussians. Ignored if density estimator is not an mdn. kwargs: additional custom arguments passed to downstream build functions. \"\"\" kwargs = dict ( zip ( ( \"z_score_x\" , \"z_score_y\" , \"hidden_features\" , \"num_transforms\" , \"num_bins\" , \"embedding_net\" , \"num_components\" , ), ( z_score_x , z_score_theta , hidden_features , num_transforms , num_bins , embedding_net , num_components , ), ), ** kwargs , ) def build_fn ( batch_theta , batch_x ): if model == \"mdn\" : return build_mdn ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"made\" : return build_made ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"maf\" : return build_maf ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"maf_rqs\" : return build_maf_rqs ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"nsf\" : return build_nsf ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"mnle\" : return build_mnle ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) else : raise NotImplementedError return build_fn","title":"likelihood_nn()"},{"location":"reference/#sbi.utils.get_nn_models.classifier_nn","text":"Returns a function that builds a classifier for learning density ratios. This function will usually be used for SNRE. The returned function is to be passed to the inference class when using the flexible interface. Note that in the view of the SNRE classifier we build below, x=theta and y=x. Parameters: Name Type Description Default model str The type of classifier that will be created. One of [ linear , mlp , resnet ]. required z_score_theta Optional[str] Whether to z-score parameters \\(\\theta\\) before passing them into the network, can take one of the following: - none , or None: do not z-score. - independent : z-score each dimension independently. - structured : treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. 'independent' z_score_x Optional[str] Whether to z-score simulation outputs \\(x\\) before passing them into the network, same options as z_score_theta. 'independent' hidden_features int Number of hidden features. 50 embedding_net_theta Module Optional embedding network for parameters \\(\\theta\\) . Identity() embedding_net_x Module Optional embedding network for simulation outputs \\(x\\) . This embedding net allows to learn features from potentially high-dimensional simulation outputs. Identity() kwargs additional custom arguments passed to downstream build functions. {} Source code in sbi/utils/get_nn_models.py def classifier_nn ( model : str , z_score_theta : Optional [ str ] = \"independent\" , z_score_x : Optional [ str ] = \"independent\" , hidden_features : int = 50 , embedding_net_theta : nn . Module = nn . Identity (), embedding_net_x : nn . Module = nn . Identity (), ** kwargs , ) -> Callable : r \"\"\" Returns a function that builds a classifier for learning density ratios. This function will usually be used for SNRE. The returned function is to be passed to the inference class when using the flexible interface. Note that in the view of the SNRE classifier we build below, x=theta and y=x. Args: model: The type of classifier that will be created. One of [`linear`, `mlp`, `resnet`]. z_score_theta: Whether to z-score parameters $\\theta$ before passing them into the network, can take one of the following: - `none`, or None: do not z-score. - `independent`: z-score each dimension independently. - `structured`: treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. z_score_x: Whether to z-score simulation outputs $x$ before passing them into the network, same options as z_score_theta. hidden_features: Number of hidden features. embedding_net_theta: Optional embedding network for parameters $\\theta$. embedding_net_x: Optional embedding network for simulation outputs $x$. This embedding net allows to learn features from potentially high-dimensional simulation outputs. kwargs: additional custom arguments passed to downstream build functions. \"\"\" kwargs = dict ( zip ( ( \"z_score_x\" , \"z_score_y\" , \"hidden_features\" , \"embedding_net_x\" , \"embedding_net_y\" , ), ( z_score_theta , z_score_x , hidden_features , embedding_net_theta , embedding_net_x , ), ), ** kwargs , ) def build_fn ( batch_theta , batch_x ): if model == \"linear\" : return build_linear_classifier ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) if model == \"mlp\" : return build_mlp_classifier ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) if model == \"resnet\" : return build_resnet_classifier ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) else : raise NotImplementedError return build_fn","title":"classifier_nn()"},{"location":"reference/#potentials","text":"","title":"Potentials"},{"location":"reference/#sbi.inference.potentials.posterior_based_potential.posterior_estimator_based_potential","text":"Returns the potential for posterior-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. The potential is the same as the log-probability of the posterior_estimator , but it is set to \\(-\\inf\\) outside of the prior bounds. Parameters: Name Type Description Default posterior_estimator Module The neural network modelling the posterior. required prior Distribution The prior distribution. required x_o Optional[torch.Tensor] The observed data at which to evaluate the posterior. required enable_transform bool Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for theta_transform . True Returns: Type Description Tuple[Callable, torch Transform] The potential function and a transformation that maps to unconstrained space. Source code in sbi/inference/potentials/posterior_based_potential.py def posterior_estimator_based_potential ( posterior_estimator : nn . Module , prior : Distribution , x_o : Optional [ Tensor ], enable_transform : bool = True , ) -> Tuple [ Callable , TorchTransform ]: r \"\"\"Returns the potential for posterior-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. The potential is the same as the log-probability of the `posterior_estimator`, but it is set to $-\\inf$ outside of the prior bounds. Args: posterior_estimator: The neural network modelling the posterior. prior: The prior distribution. x_o: The observed data at which to evaluate the posterior. enable_transform: Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for `theta_transform`. Returns: The potential function and a transformation that maps to unconstrained space. \"\"\" device = str ( next ( posterior_estimator . parameters ()) . device ) potential_fn = PosteriorBasedPotential ( posterior_estimator , prior , x_o , device = device ) theta_transform = mcmc_transform ( prior , device = device , enable_transform = enable_transform ) return potential_fn , theta_transform","title":"posterior_estimator_based_potential()"},{"location":"reference/#sbi.inference.potentials.likelihood_based_potential.likelihood_estimator_based_potential","text":"Returns potential \\(\\log(p(x_o|\\theta)p(\\theta))\\) for likelihood-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. Parameters: Name Type Description Default likelihood_estimator Module The neural network modelling the likelihood. required prior Distribution The prior distribution. required x_o Optional[torch.Tensor] The observed data at which to evaluate the likelihood. required enable_transform bool Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for theta_transform . True Returns: Type Description Tuple[Callable, torch Transform] The potential function \\(p(x_o|\\theta)p(\\theta)\\) and a transformation that maps to unconstrained space. Source code in sbi/inference/potentials/likelihood_based_potential.py def likelihood_estimator_based_potential ( likelihood_estimator : nn . Module , prior : Distribution , x_o : Optional [ Tensor ], enable_transform : bool = True , ) -> Tuple [ Callable , TorchTransform ]: r \"\"\"Returns potential $\\log(p(x_o|\\theta)p(\\theta))$ for likelihood-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. Args: likelihood_estimator: The neural network modelling the likelihood. prior: The prior distribution. x_o: The observed data at which to evaluate the likelihood. enable_transform: Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for `theta_transform`. Returns: The potential function $p(x_o|\\theta)p(\\theta)$ and a transformation that maps to unconstrained space. \"\"\" device = str ( next ( likelihood_estimator . parameters ()) . device ) potential_fn = LikelihoodBasedPotential ( likelihood_estimator , prior , x_o , device = device ) theta_transform = mcmc_transform ( prior , device = device , enable_transform = enable_transform ) return potential_fn , theta_transform","title":"likelihood_estimator_based_potential()"},{"location":"reference/#sbi.inference.potentials.ratio_based_potential.ratio_estimator_based_potential","text":"Returns the potential for ratio-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. Parameters: Name Type Description Default ratio_estimator Module The neural network modelling likelihood-to-evidence ratio. required prior Distribution The prior distribution. required x_o Optional[torch.Tensor] The observed data at which to evaluate the likelihood-to-evidence ratio. required enable_transform bool Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for theta_transform . True Returns: Type Description Tuple[Callable, torch Transform] The potential function and a transformation that maps to unconstrained space. Source code in sbi/inference/potentials/ratio_based_potential.py def ratio_estimator_based_potential ( ratio_estimator : nn . Module , prior : Distribution , x_o : Optional [ Tensor ], enable_transform : bool = True , ) -> Tuple [ Callable , TorchTransform ]: r \"\"\"Returns the potential for ratio-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. Args: ratio_estimator: The neural network modelling likelihood-to-evidence ratio. prior: The prior distribution. x_o: The observed data at which to evaluate the likelihood-to-evidence ratio. enable_transform: Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for `theta_transform`. Returns: The potential function and a transformation that maps to unconstrained space. \"\"\" device = str ( next ( ratio_estimator . parameters ()) . device ) potential_fn = RatioBasedPotential ( ratio_estimator , prior , x_o , device = device ) theta_transform = mcmc_transform ( prior , device = device , enable_transform = enable_transform ) return potential_fn , theta_transform","title":"ratio_estimator_based_potential()"},{"location":"reference/#analysis","text":"","title":"Analysis"},{"location":"reference/#sbi.analysis.plot.pairplot","text":"Plot samples in a 2D grid showing marginals and pairwise marginals. Each of the diagonal plots can be interpreted as a 1D-marginal of the distribution that the samples were drawn from. Each upper-diagonal plot can be interpreted as a 2D-marginal of the distribution. Parameters: Name Type Description Default samples Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] Samples used to build the histogram. required points Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] List of additional points to scatter. None limits Union[List, torch.Tensor] Array containing the plot xlim for each parameter dimension. If None, just use the min and max of the passed samples None subset Optional[List[int]] List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1 st and 3 rd dimension but will discard the 0 th and 2 nd (and, if they exist, the 4 th , 5 th and so on). None offdiag Union[str, List[str]] Plotting style for upper diagonal, {hist, scatter, contour, cond, None}. 'hist' upper Optional[str] deprecated, use offdiag instead. None diag Union[str, List[str]] Plotting style for diagonal, {hist, cond, None}. 'hist' figsize Tuple Size of the entire figure. (10, 10) labels Optional[List[str]] List of strings specifying the names of the parameters. None ticks Union[List, torch.Tensor] Position of the ticks. [] fig matplotlib figure to plot on. None axes matplotlib axes corresponding to fig. None **kwargs Additional arguments to adjust the plot, e.g., samples_colors , points_colors and many more, see the source code in _get_default_opts() in sbi.analysis.plot for details. {} Returns: figure and axis of posterior distribution plot Source code in sbi/analysis/plot.py def pairplot ( samples : Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ], points : Optional [ Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ] ] = None , limits : Optional [ Union [ List , torch . Tensor ]] = None , subset : Optional [ List [ int ]] = None , offdiag : Optional [ Union [ List [ str ], str ]] = \"hist\" , diag : Optional [ Union [ List [ str ], str ]] = \"hist\" , figsize : Tuple = ( 10 , 10 ), labels : Optional [ List [ str ]] = None , ticks : Union [ List , torch . Tensor ] = [], upper : Optional [ str ] = None , fig = None , axes = None , ** kwargs , ): \"\"\" Plot samples in a 2D grid showing marginals and pairwise marginals. Each of the diagonal plots can be interpreted as a 1D-marginal of the distribution that the samples were drawn from. Each upper-diagonal plot can be interpreted as a 2D-marginal of the distribution. Args: samples: Samples used to build the histogram. points: List of additional points to scatter. limits: Array containing the plot xlim for each parameter dimension. If None, just use the min and max of the passed samples subset: List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1st and 3rd dimension but will discard the 0th and 2nd (and, if they exist, the 4th, 5th and so on). offdiag: Plotting style for upper diagonal, {hist, scatter, contour, cond, None}. upper: deprecated, use offdiag instead. diag: Plotting style for diagonal, {hist, cond, None}. figsize: Size of the entire figure. labels: List of strings specifying the names of the parameters. ticks: Position of the ticks. fig: matplotlib figure to plot on. axes: matplotlib axes corresponding to fig. **kwargs: Additional arguments to adjust the plot, e.g., `samples_colors`, `points_colors` and many more, see the source code in `_get_default_opts()` in `sbi.analysis.plot` for details. Returns: figure and axis of posterior distribution plot \"\"\" # TODO: add color map support # TODO: automatically determine good bin sizes for histograms # TODO: add legend (if legend is True) opts = _get_default_opts () # update the defaults dictionary by the current values of the variables (passed by # the user) opts = _update ( opts , locals ()) opts = _update ( opts , kwargs ) samples , dim , limits = prepare_for_plot ( samples , limits ) # checks. if opts [ \"legend\" ]: assert len ( opts [ \"samples_labels\" ]) >= len ( samples ), \"Provide at least as many labels as samples.\" if opts [ \"upper\" ] is not None : warn ( \"upper is deprecated, use offdiag instead.\" ) opts [ \"offdiag\" ] = opts [ \"upper\" ] # Prepare diag/upper/lower if type ( opts [ \"diag\" ]) is not list : opts [ \"diag\" ] = [ opts [ \"diag\" ] for _ in range ( len ( samples ))] if type ( opts [ \"offdiag\" ]) is not list : opts [ \"offdiag\" ] = [ opts [ \"offdiag\" ] for _ in range ( len ( samples ))] # if type(opts['lower']) is not list: # opts['lower'] = [opts['lower'] for _ in range(len(samples))] opts [ \"lower\" ] = None diag_func = get_diag_func ( samples , limits , opts , ** kwargs ) def offdiag_func ( row , col , limits , ** kwargs ): if len ( samples ) > 0 : for n , v in enumerate ( samples ): if opts [ \"offdiag\" ][ n ] == \"hist\" or opts [ \"offdiag\" ][ n ] == \"hist2d\" : hist , xedges , yedges = np . histogram2d ( v [:, col ], v [:, row ], range = [ [ limits [ col ][ 0 ], limits [ col ][ 1 ]], [ limits [ row ][ 0 ], limits [ row ][ 1 ]], ], ** opts [ \"hist_offdiag\" ], ) plt . imshow ( hist . T , origin = \"lower\" , extent = ( xedges [ 0 ], xedges [ - 1 ], yedges [ 0 ], yedges [ - 1 ], ), aspect = \"auto\" , ) elif opts [ \"offdiag\" ][ n ] in [ \"kde\" , \"kde2d\" , \"contour\" , \"contourf\" , ]: density = gaussian_kde ( v [:, [ col , row ]] . T , bw_method = opts [ \"kde_offdiag\" ][ \"bw_method\" ], ) X , Y = np . meshgrid ( np . linspace ( limits [ col ][ 0 ], limits [ col ][ 1 ], opts [ \"kde_offdiag\" ][ \"bins\" ], ), np . linspace ( limits [ row ][ 0 ], limits [ row ][ 1 ], opts [ \"kde_offdiag\" ][ \"bins\" ], ), ) positions = np . vstack ([ X . ravel (), Y . ravel ()]) Z = np . reshape ( density ( positions ) . T , X . shape ) if opts [ \"offdiag\" ][ n ] == \"kde\" or opts [ \"offdiag\" ][ n ] == \"kde2d\" : plt . imshow ( Z , extent = ( limits [ col ][ 0 ], limits [ col ][ 1 ], limits [ row ][ 0 ], limits [ row ][ 1 ], ), origin = \"lower\" , aspect = \"auto\" , ) elif opts [ \"offdiag\" ][ n ] == \"contour\" : if opts [ \"contour_offdiag\" ][ \"percentile\" ]: Z = probs2contours ( Z , opts [ \"contour_offdiag\" ][ \"levels\" ]) else : Z = ( Z - Z . min ()) / ( Z . max () - Z . min ()) plt . contour ( X , Y , Z , origin = \"lower\" , extent = [ limits [ col ][ 0 ], limits [ col ][ 1 ], limits [ row ][ 0 ], limits [ row ][ 1 ], ], colors = opts [ \"samples_colors\" ][ n ], levels = opts [ \"contour_offdiag\" ][ \"levels\" ], ) else : pass elif opts [ \"offdiag\" ][ n ] == \"scatter\" : plt . scatter ( v [:, col ], v [:, row ], color = opts [ \"samples_colors\" ][ n ], ** opts [ \"scatter_offdiag\" ], ) elif opts [ \"offdiag\" ][ n ] == \"plot\" : plt . plot ( v [:, col ], v [:, row ], color = opts [ \"samples_colors\" ][ n ], ** opts [ \"plot_offdiag\" ], ) else : pass return _arrange_plots ( diag_func , offdiag_func , dim , limits , points , opts , fig = fig , axes = axes )","title":"pairplot()"},{"location":"reference/#sbi.analysis.plot.marginal_plot","text":"Plot samples in a row showing 1D marginals of selected dimensions. Each of the plots can be interpreted as a 1D-marginal of the distribution that the samples were drawn from. Parameters: Name Type Description Default samples Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] Samples used to build the histogram. required points Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] List of additional points to scatter. None limits Union[List, torch.Tensor] Array containing the plot xlim for each parameter dimension. If None, just use the min and max of the passed samples None subset Optional[List[int]] List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1 st and 3 rd dimension but will discard the 0 th and 2 nd (and, if they exist, the 4 th , 5 th and so on). None diag Optional[str] Plotting style for 1D marginals, {hist, kde cond, None}. 'hist' figsize Tuple Size of the entire figure. (10, 10) labels Optional[List[str]] List of strings specifying the names of the parameters. None ticks Union[List, torch.Tensor] Position of the ticks. [] points_colors Colors of the points . required fig matplotlib figure to plot on. None axes matplotlib axes corresponding to fig. None **kwargs Additional arguments to adjust the plot, e.g., samples_colors , points_colors and many more, see the source code in _get_default_opts() in sbi.analysis.plot for details. {} Returns: figure and axis of posterior distribution plot Source code in sbi/analysis/plot.py def marginal_plot ( samples : Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ], points : Optional [ Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ] ] = None , limits : Optional [ Union [ List , torch . Tensor ]] = None , subset : Optional [ List [ int ]] = None , diag : Optional [ str ] = \"hist\" , figsize : Tuple = ( 10 , 10 ), labels : Optional [ List [ str ]] = None , ticks : Union [ List , torch . Tensor ] = [], fig = None , axes = None , ** kwargs , ): \"\"\" Plot samples in a row showing 1D marginals of selected dimensions. Each of the plots can be interpreted as a 1D-marginal of the distribution that the samples were drawn from. Args: samples: Samples used to build the histogram. points: List of additional points to scatter. limits: Array containing the plot xlim for each parameter dimension. If None, just use the min and max of the passed samples subset: List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1st and 3rd dimension but will discard the 0th and 2nd (and, if they exist, the 4th, 5th and so on). diag: Plotting style for 1D marginals, {hist, kde cond, None}. figsize: Size of the entire figure. labels: List of strings specifying the names of the parameters. ticks: Position of the ticks. points_colors: Colors of the `points`. fig: matplotlib figure to plot on. axes: matplotlib axes corresponding to fig. **kwargs: Additional arguments to adjust the plot, e.g., `samples_colors`, `points_colors` and many more, see the source code in `_get_default_opts()` in `sbi.analysis.plot` for details. Returns: figure and axis of posterior distribution plot \"\"\" opts = _get_default_opts () # update the defaults dictionary by the current values of the variables (passed by # the user) opts = _update ( opts , locals ()) opts = _update ( opts , kwargs ) samples , dim , limits = prepare_for_plot ( samples , limits ) # Prepare diag/upper/lower if type ( opts [ \"diag\" ]) is not list : opts [ \"diag\" ] = [ opts [ \"diag\" ] for _ in range ( len ( samples ))] diag_func = get_diag_func ( samples , limits , opts , ** kwargs ) return _arrange_plots ( diag_func , None , dim , limits , points , opts , fig = fig , axes = axes )","title":"marginal_plot()"},{"location":"reference/#sbi.analysis.plot.conditional_pairplot","text":"Plot conditional distribution given all other parameters. The conditionals can be interpreted as slices through the density at a location given by condition . For example: Say we have a 3D density with parameters \\(\\theta_0\\) , \\(\\theta_1\\) , \\(\\theta_2\\) and a condition \\(c\\) passed by the user in the condition argument. For the plot of \\(\\theta_0\\) on the diagonal, this will plot the conditional \\(p(\\theta_0 | \\theta_1=c[1], \\theta_2=c[2])\\) . For the upper diagonal of \\(\\theta_1\\) and \\(\\theta_2\\) , it will plot \\(p(\\theta_1, \\theta_2 | \\theta_0=c[0])\\) . All other diagonals and upper-diagonals are built in the corresponding way. Parameters: Name Type Description Default density Any Probability density with a log_prob() method. required condition Tensor Condition that all but the one/two regarded parameters are fixed to. The condition should be of shape (1, dim_theta), i.e. it could e.g. be a sample from the posterior distribution. required limits Union[List, torch.Tensor] Limits in between which each parameter will be evaluated. required points Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] Additional points to scatter. None subset Optional[List[int]] List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1 st and 3 rd dimension but will discard the 0 th and 2 nd (and, if they exist, the 4 th , 5 th and so on) None resolution int Resolution of the grid at which we evaluate the pdf . 50 figsize Tuple Size of the entire figure. (10, 10) labels Optional[List[str]] List of strings specifying the names of the parameters. None ticks Union[List, torch.Tensor] Position of the ticks. [] points_colors Colors of the points . required fig matplotlib figure to plot on. None axes matplotlib axes corresponding to fig. None **kwargs Additional arguments to adjust the plot, e.g., samples_colors , points_colors and many more, see the source code in _get_default_opts() in sbi.analysis.plot for details. {} Returns: figure and axis of posterior distribution plot Source code in sbi/analysis/plot.py def conditional_pairplot ( density : Any , condition : torch . Tensor , limits : Union [ List , torch . Tensor ], points : Optional [ Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ] ] = None , subset : Optional [ List [ int ]] = None , resolution : int = 50 , figsize : Tuple = ( 10 , 10 ), labels : Optional [ List [ str ]] = None , ticks : Union [ List , torch . Tensor ] = [], fig = None , axes = None , ** kwargs , ): r \"\"\" Plot conditional distribution given all other parameters. The conditionals can be interpreted as slices through the `density` at a location given by `condition`. For example: Say we have a 3D density with parameters $\\theta_0$, $\\theta_1$, $\\theta_2$ and a condition $c$ passed by the user in the `condition` argument. For the plot of $\\theta_0$ on the diagonal, this will plot the conditional $p(\\theta_0 | \\theta_1=c[1], \\theta_2=c[2])$. For the upper diagonal of $\\theta_1$ and $\\theta_2$, it will plot $p(\\theta_1, \\theta_2 | \\theta_0=c[0])$. All other diagonals and upper-diagonals are built in the corresponding way. Args: density: Probability density with a `log_prob()` method. condition: Condition that all but the one/two regarded parameters are fixed to. The condition should be of shape (1, dim_theta), i.e. it could e.g. be a sample from the posterior distribution. limits: Limits in between which each parameter will be evaluated. points: Additional points to scatter. subset: List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1st and 3rd dimension but will discard the 0th and 2nd (and, if they exist, the 4th, 5th and so on) resolution: Resolution of the grid at which we evaluate the `pdf`. figsize: Size of the entire figure. labels: List of strings specifying the names of the parameters. ticks: Position of the ticks. points_colors: Colors of the `points`. fig: matplotlib figure to plot on. axes: matplotlib axes corresponding to fig. **kwargs: Additional arguments to adjust the plot, e.g., `samples_colors`, `points_colors` and many more, see the source code in `_get_default_opts()` in `sbi.analysis.plot` for details. Returns: figure and axis of posterior distribution plot \"\"\" device = density . _device if hasattr ( density , \"_device\" ) else \"cpu\" # Setting these is required because _pairplot_scaffold will check if opts['diag'] is # `None`. This would break if opts has no key 'diag'. Same for 'upper'. diag = \"cond\" offdiag = \"cond\" opts = _get_default_opts () # update the defaults dictionary by the current values of the variables (passed by # the user) opts = _update ( opts , locals ()) opts = _update ( opts , kwargs ) opts [ \"lower\" ] = None dim , limits , eps_margins = prepare_for_conditional_plot ( condition , opts ) diag_func = get_conditional_diag_func ( opts , limits , eps_margins , resolution ) def offdiag_func ( row , col , ** kwargs ): p_image = ( eval_conditional_density ( opts [ \"density\" ], opts [ \"condition\" ] . to ( device ), limits . to ( device ), row , col , resolution = resolution , eps_margins1 = eps_margins [ row ], eps_margins2 = eps_margins [ col ], ) . to ( \"cpu\" ) . numpy () ) plt . imshow ( p_image . T , origin = \"lower\" , extent = ( limits [ col , 0 ] . item (), limits [ col , 1 ] . item (), limits [ row , 0 ] . item (), limits [ row , 1 ] . item (), ), aspect = \"auto\" , ) return _arrange_plots ( diag_func , offdiag_func , dim , limits , points , opts , fig = fig , axes = axes )","title":"conditional_pairplot()"},{"location":"reference/#sbi.analysis.conditional_density.conditional_corrcoeff","text":"Returns the conditional correlation matrix of a distribution. To compute the conditional distribution, we condition all but two parameters to values from condition , and then compute the Pearson correlation coefficient \\(\\rho\\) between the remaining two parameters under the distribution density . We do so for any pair of parameters specified in subset , thus creating a matrix containing conditional correlations between any pair of parameters. If condition is a batch of conditions, this function computes the conditional correlation matrix for each one of them and returns the mean. Parameters: Name Type Description Default density Any Probability density function with .log_prob() function. required limits Tensor Limits within which to evaluate the density . required condition Tensor Values to condition the density on. If a batch of conditions is passed, we compute the conditional correlation matrix for each of them and return the average conditional correlation matrix. required subset Optional[List[int]] Evaluate the conditional distribution only on a subset of dimensions. If None this function uses all dimensions. None resolution int Number of grid points on which the conditional distribution is evaluated. A higher value increases the accuracy of the estimated correlation but also increases the computational cost. 50 Returns: Average conditional correlation matrix of shape either (num_dim, num_dim) or (len(subset), len(subset)) if subset was specified. Source code in sbi/analysis/conditional_density.py def conditional_corrcoeff ( density : Any , limits : Tensor , condition : Tensor , subset : Optional [ List [ int ]] = None , resolution : int = 50 , ) -> Tensor : r \"\"\"Returns the conditional correlation matrix of a distribution. To compute the conditional distribution, we condition all but two parameters to values from `condition`, and then compute the Pearson correlation coefficient $\\rho$ between the remaining two parameters under the distribution `density`. We do so for any pair of parameters specified in `subset`, thus creating a matrix containing conditional correlations between any pair of parameters. If `condition` is a batch of conditions, this function computes the conditional correlation matrix for each one of them and returns the mean. Args: density: Probability density function with `.log_prob()` function. limits: Limits within which to evaluate the `density`. condition: Values to condition the `density` on. If a batch of conditions is passed, we compute the conditional correlation matrix for each of them and return the average conditional correlation matrix. subset: Evaluate the conditional distribution only on a subset of dimensions. If `None` this function uses all dimensions. resolution: Number of grid points on which the conditional distribution is evaluated. A higher value increases the accuracy of the estimated correlation but also increases the computational cost. Returns: Average conditional correlation matrix of shape either `(num_dim, num_dim)` or `(len(subset), len(subset))` if `subset` was specified. \"\"\" device = density . _device if hasattr ( density , \"_device\" ) else \"cpu\" subset_ = subset if subset is not None else range ( condition . shape [ 1 ]) correlation_matrices = [] for cond in condition : correlation_matrices . append ( torch . stack ( [ compute_corrcoeff ( eval_conditional_density ( density , cond . to ( device ), limits . to ( device ), dim1 = dim1 , dim2 = dim2 , resolution = resolution , ), limits [[ dim1 , dim2 ]] . to ( device ), ) for dim1 in subset_ for dim2 in subset_ if dim1 < dim2 ] ) ) average_correlations = torch . mean ( torch . stack ( correlation_matrices ), dim = 0 ) # `average_correlations` is still a vector containing the upper triangular entries. # Below, assemble them into a matrix: av_correlation_matrix = torch . zeros (( len ( subset_ ), len ( subset_ )), device = device ) triu_indices = torch . triu_indices ( row = len ( subset_ ), col = len ( subset_ ), offset = 1 , device = device ) av_correlation_matrix [ triu_indices [ 0 ], triu_indices [ 1 ]] = average_correlations # Make the matrix symmetric by copying upper diagonal to lower diagonal. av_correlation_matrix = torch . triu ( av_correlation_matrix ) + torch . tril ( av_correlation_matrix . T ) av_correlation_matrix . fill_diagonal_ ( 1.0 ) return av_correlation_matrix","title":"conditional_corrcoeff()"},{"location":"examples/00_HH_simulator/","text":"Inference on Hodgkin-Huxley model: tutorial \u00b6 In this tutorial, we use sbi to do inference on a Hodgkin-Huxley model from neuroscience (Hodgkin and Huxley, 1952). We will learn two parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ) based on a current-clamp recording, that we generate synthetically (in practice, this would be an experimental observation). Note, you find the original version of this notebook at https://github.com/sbi-dev/sbi/blob/main/examples/00_HH_simulator.ipynb in the sbi repository. First we are going to import basic packages. import numpy as np import torch # visualization import matplotlib as mpl import matplotlib.pyplot as plt # sbi from sbi import utils as utils from sbi import analysis as analysis from sbi.inference.base import infer # remove top and right axis from plots mpl . rcParams [ \"axes.spines.right\" ] = False mpl . rcParams [ \"axes.spines.top\" ] = False Different required components \u00b6 Before running inference, let us define the different required components: observed data prior over model parameters simulator 1. Observed data \u00b6 Let us assume we current-clamped a neuron and recorded the following voltage trace: In fact, this voltage trace was not measured experimentally but synthetically generated by simulating a Hodgkin-Huxley model with particular parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ). We will come back to this point later in the tutorial. 2. Simulator \u00b6 We would like to infer the posterior over the two parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ) of a Hodgkin-Huxley model, given the observed electrophysiological recording above. The model has channel kinetics as in Pospischil et al. 2008 , and is defined by the following set of differential equations (parameters of interest highlighted in orange): \\[ \\scriptsize \\begin{align} C_m\\frac{dV}{dt}&=g_1\\left(E_1-V\\right)+ \\color{orange}{\\bar{g}_{Na}}m^3h\\left(E_{Na}-V\\right)+ \\color{orange}{\\bar{g}_{K}}n^4\\left(E_K-V\\right)+ \\bar{g}_Mp\\left(E_K-V\\right)+ I_{inj}+ \\sigma\\eta\\left(t\\right)\\\\ \\frac{dq}{dt}&=\\frac{q_\\infty\\left(V\\right)-q}{\\tau_q\\left(V\\right)},\\;q\\in\\{m,h,n,p\\} \\end{align} \\] Above, \\(V\\) represents the membrane potential, \\(C_m\\) is the membrane capacitance, \\(g_{\\text{l}}\\) is the leak conductance, \\(E_{\\text{l}}\\) is the membrane reversal potential, \\(\\bar{g}_c\\) is the density of channels of type \\(c\\) ( \\(\\text{Na}^+\\) , \\(\\text{K}^+\\) , M), \\(E_c\\) is the reversal potential of \\(c\\) , ( \\(m\\) , \\(h\\) , \\(n\\) , \\(p\\) ) are the respective channel gating kinetic variables, and \\(\\sigma \\eta(t)\\) is the intrinsic neural noise. The right hand side of the voltage dynamics is composed of a leak current, a voltage-dependent \\(\\text{Na}^+\\) current, a delayed-rectifier \\(\\text{K}^+\\) current, a slow voltage-dependent \\(\\text{K}^+\\) current responsible for spike-frequency adaptation, and an injected current \\(I_{\\text{inj}}\\) . Channel gating variables \\(q\\) have dynamics fully characterized by the neuron membrane potential \\(V\\) , given the respective steady-state \\(q_{\\infty}(V)\\) and time constant \\(\\tau_{q}(V)\\) (details in Pospischil et al. 2008). The input current \\(I_{\\text{inj}}\\) is defined as from HH_helper_functions import syn_current I , t_on , t_off , dt , t , A_soma = syn_current () The Hodgkin-Huxley simulator is given by: from HH_helper_functions import HHsimulator Putting the input current and the simulator together: def run_HH_model ( params ): params = np . asarray ( params ) # input current, time step I , t_on , t_off , dt , t , A_soma = syn_current () t = np . arange ( 0 , len ( I ), 1 ) * dt # initial voltage V0 = - 70 states = HHsimulator ( V0 , params . reshape ( 1 , - 1 ), dt , t , I ) return dict ( data = states . reshape ( - 1 ), time = t , dt = dt , I = I . reshape ( - 1 )) To get an idea of the output of the Hodgkin-Huxley model, let us generate some voltage traces for different parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), given the input current \\(I_{\\text{inj}}\\) : # three sets of (g_Na, g_K) params = np . array ([[ 50.0 , 1.0 ], [ 4.0 , 1.5 ], [ 20.0 , 15.0 ]]) num_samples = len ( params [:, 0 ]) sim_samples = np . zeros (( num_samples , len ( I ))) for i in range ( num_samples ): sim_samples [ i , :] = run_HH_model ( params = params [ i , :])[ \"data\" ] # colors for traces col_min = 2 num_colors = num_samples + col_min cm1 = mpl . cm . Blues col1 = [ cm1 ( 1.0 * i / num_colors ) for i in range ( col_min , num_colors )] fig = plt . figure ( figsize = ( 7 , 5 )) gs = mpl . gridspec . GridSpec ( 2 , 1 , height_ratios = [ 4 , 1 ]) ax = plt . subplot ( gs [ 0 ]) for i in range ( num_samples ): plt . plot ( t , sim_samples [ i , :], color = col1 [ i ], lw = 2 ) plt . ylabel ( \"voltage (mV)\" ) ax . set_xticks ([]) ax . set_yticks ([ - 80 , - 20 , 40 ]) ax = plt . subplot ( gs [ 1 ]) plt . plot ( t , I * A_soma * 1e3 , \"k\" , lw = 2 ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"input (nA)\" ) ax . set_xticks ([ 0 , max ( t ) / 2 , max ( t )]) ax . set_yticks ([ 0 , 1.1 * np . max ( I * A_soma * 1e3 )]) ax . yaxis . set_major_formatter ( mpl . ticker . FormatStrFormatter ( \" %.2f \" )) plt . show () As can be seen, the voltage traces can be quite diverse for different parameter values. Often, we are not interested in matching the exact trace, but only in matching certain features thereof. In this example of the Hodgkin-Huxley model, the summary features are the number of spikes, the mean resting potential, the standard deviation of the resting potential, and the first four voltage moments: mean, standard deviation, skewness and kurtosis. Using the function calculate_summary_statistics() imported below, we obtain these statistics from the output of the Hodgkin Huxley simulator. from HH_helper_functions import calculate_summary_statistics Lastly, we define a function that performs all of the above steps at once. The function simulation_wrapper takes in conductance values, runs the Hodgkin Huxley model and then returns the summary statistics. def simulation_wrapper ( params ): \"\"\" Returns summary statistics from conductance values in `params`. Summarizes the output of the HH simulator and converts it to `torch.Tensor`. \"\"\" obs = run_HH_model ( params ) summstats = torch . as_tensor ( calculate_summary_statistics ( obs )) return summstats sbi takes any function as simulator. Thus, sbi also has the flexibility to use simulators that utilize external packages, e.g., Brian ( http://briansimulator.org/ ), nest ( https://www.nest-simulator.org/ ), or NEURON ( https://neuron.yale.edu/neuron/ ). External simulators do not even need to be Python-based as long as they store simulation outputs in a format that can be read from Python. All that is necessary is to wrap your external simulator of choice into a Python callable that takes a parameter set and outputs a set of summary statistics we want to fit the parameters to. 3. Prior over model parameters \u00b6 Now that we have the simulator, we need to define a function with the prior over the model parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), which in this case is chosen to be a Uniform distribution: prior_min = [ 0.5 , 1e-4 ] prior_max = [ 80.0 , 15.0 ] prior = utils . torchutils . BoxUniform ( low = torch . as_tensor ( prior_min ), high = torch . as_tensor ( prior_max ) ) Inference \u00b6 Now that we have all the required components, we can run inference with SNPE to identify parameters whose activity matches this trace. posterior = infer ( simulation_wrapper , prior , method = \"SNPE\" , num_simulations = 300 , num_workers = 4 ) HBox(children=(FloatProgress(value=0.0, description='Running 300 simulations in 300 batches.', max=300.0, styl\u2026 Neural network successfully converged after 233 epochs. Note sbi can parallelize your simulator. If you experience problems with parallelization, try setting num_workers=1 and please give us an error report as a GitHub issue . Coming back to the observed data \u00b6 As mentioned at the beginning of the tutorial, the observed data are generated by the Hodgkin-Huxley model with a set of known parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ). To illustrate how to compute the summary statistics of the observed data, let us regenerate the observed data: # true parameters and respective labels true_params = np . array ([ 50.0 , 5.0 ]) labels_params = [ r \"$g_ {Na} $\" , r \"$g_ {K} $\" ] observation_trace = run_HH_model ( true_params ) observation_summary_statistics = calculate_summary_statistics ( observation_trace ) As we already shown above, the observed voltage traces look as follows: fig = plt . figure ( figsize = ( 7 , 5 )) gs = mpl . gridspec . GridSpec ( 2 , 1 , height_ratios = [ 4 , 1 ]) ax = plt . subplot ( gs [ 0 ]) plt . plot ( observation_trace [ \"time\" ], observation_trace [ \"data\" ]) plt . ylabel ( \"voltage (mV)\" ) plt . title ( \"observed data\" ) plt . setp ( ax , xticks = [], yticks = [ - 80 , - 20 , 40 ]) ax = plt . subplot ( gs [ 1 ]) plt . plot ( observation_trace [ \"time\" ], I * A_soma * 1e3 , \"k\" , lw = 2 ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"input (nA)\" ) ax . set_xticks ([ 0 , max ( observation_trace [ \"time\" ]) / 2 , max ( observation_trace [ \"time\" ])]) ax . set_yticks ([ 0 , 1.1 * np . max ( I * A_soma * 1e3 )]) ax . yaxis . set_major_formatter ( mpl . ticker . FormatStrFormatter ( \" %.2f \" )) Analysis of the posterior given the observed data \u00b6 After running the inference algorithm, let us inspect the inferred posterior distribution over the parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), given the observed trace. To do so, we first draw samples (i.e. consistent parameter sets) from the posterior: samples = posterior . sample (( 10000 ,), x = observation_summary_statistics ) HBox(children=(FloatProgress(value=0.0, description='Drawing 10000 posterior samples', max=10000.0, style=Prog\u2026 fig , axes = analysis . pairplot ( samples , limits = [[ 0.5 , 80 ], [ 1e-4 , 15.0 ]], ticks = [[ 0.5 , 80 ], [ 1e-4 , 15.0 ]], figsize = ( 5 , 5 ), points = true_params , points_offdiag = { \"markersize\" : 6 }, points_colors = \"r\" , ); As can be seen, the inferred posterior contains the ground-truth parameters (red) in a high-probability region. Now, let us sample parameters from the posterior distribution, simulate the Hodgkin-Huxley model for this parameter set and compare the simulations with the observed data: # Draw a sample from the posterior and convert to numpy for plotting. posterior_sample = posterior . sample (( 1 ,), x = observation_summary_statistics ) . numpy () HBox(children=(FloatProgress(value=0.0, description='Drawing 1 posterior samples', max=1.0, style=ProgressStyl\u2026 fig = plt . figure ( figsize = ( 7 , 5 )) # plot observation t = observation_trace [ \"time\" ] y_obs = observation_trace [ \"data\" ] plt . plot ( t , y_obs , lw = 2 , label = \"observation\" ) # simulate and plot samples x = run_HH_model ( posterior_sample ) plt . plot ( t , x [ \"data\" ], \"--\" , lw = 2 , label = \"posterior sample\" ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"voltage (mV)\" ) ax = plt . gca () handles , labels = ax . get_legend_handles_labels () ax . legend ( handles [:: - 1 ], labels [:: - 1 ], bbox_to_anchor = ( 1.3 , 1 ), loc = \"upper right\" ) ax . set_xticks ([ 0 , 60 , 120 ]) ax . set_yticks ([ - 80 , - 20 , 40 ]); As can be seen, the sample from the inferred posterior leads to simulations that closely resemble the observed data, confirming that SNPE did a good job at capturing the observed data in this simple case. References \u00b6 A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117(4):500\u2013544, 1952. M. Pospischil, M. Toledo-Rodriguez, C. Monier, Z. Piwkowska, T. Bal, Y. Fr\u00e9gnac, H. Markram, and A. Destexhe. Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biological Cybernetics, 99(4-5), 2008.","title":"Hodgkin-Huxley example"},{"location":"examples/00_HH_simulator/#inference-on-hodgkin-huxley-model-tutorial","text":"In this tutorial, we use sbi to do inference on a Hodgkin-Huxley model from neuroscience (Hodgkin and Huxley, 1952). We will learn two parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ) based on a current-clamp recording, that we generate synthetically (in practice, this would be an experimental observation). Note, you find the original version of this notebook at https://github.com/sbi-dev/sbi/blob/main/examples/00_HH_simulator.ipynb in the sbi repository. First we are going to import basic packages. import numpy as np import torch # visualization import matplotlib as mpl import matplotlib.pyplot as plt # sbi from sbi import utils as utils from sbi import analysis as analysis from sbi.inference.base import infer # remove top and right axis from plots mpl . rcParams [ \"axes.spines.right\" ] = False mpl . rcParams [ \"axes.spines.top\" ] = False","title":"Inference on Hodgkin-Huxley model: tutorial"},{"location":"examples/00_HH_simulator/#different-required-components","text":"Before running inference, let us define the different required components: observed data prior over model parameters simulator","title":"Different required components"},{"location":"examples/00_HH_simulator/#1-observed-data","text":"Let us assume we current-clamped a neuron and recorded the following voltage trace: In fact, this voltage trace was not measured experimentally but synthetically generated by simulating a Hodgkin-Huxley model with particular parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ). We will come back to this point later in the tutorial.","title":"1. Observed data"},{"location":"examples/00_HH_simulator/#2-simulator","text":"We would like to infer the posterior over the two parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ) of a Hodgkin-Huxley model, given the observed electrophysiological recording above. The model has channel kinetics as in Pospischil et al. 2008 , and is defined by the following set of differential equations (parameters of interest highlighted in orange): \\[ \\scriptsize \\begin{align} C_m\\frac{dV}{dt}&=g_1\\left(E_1-V\\right)+ \\color{orange}{\\bar{g}_{Na}}m^3h\\left(E_{Na}-V\\right)+ \\color{orange}{\\bar{g}_{K}}n^4\\left(E_K-V\\right)+ \\bar{g}_Mp\\left(E_K-V\\right)+ I_{inj}+ \\sigma\\eta\\left(t\\right)\\\\ \\frac{dq}{dt}&=\\frac{q_\\infty\\left(V\\right)-q}{\\tau_q\\left(V\\right)},\\;q\\in\\{m,h,n,p\\} \\end{align} \\] Above, \\(V\\) represents the membrane potential, \\(C_m\\) is the membrane capacitance, \\(g_{\\text{l}}\\) is the leak conductance, \\(E_{\\text{l}}\\) is the membrane reversal potential, \\(\\bar{g}_c\\) is the density of channels of type \\(c\\) ( \\(\\text{Na}^+\\) , \\(\\text{K}^+\\) , M), \\(E_c\\) is the reversal potential of \\(c\\) , ( \\(m\\) , \\(h\\) , \\(n\\) , \\(p\\) ) are the respective channel gating kinetic variables, and \\(\\sigma \\eta(t)\\) is the intrinsic neural noise. The right hand side of the voltage dynamics is composed of a leak current, a voltage-dependent \\(\\text{Na}^+\\) current, a delayed-rectifier \\(\\text{K}^+\\) current, a slow voltage-dependent \\(\\text{K}^+\\) current responsible for spike-frequency adaptation, and an injected current \\(I_{\\text{inj}}\\) . Channel gating variables \\(q\\) have dynamics fully characterized by the neuron membrane potential \\(V\\) , given the respective steady-state \\(q_{\\infty}(V)\\) and time constant \\(\\tau_{q}(V)\\) (details in Pospischil et al. 2008). The input current \\(I_{\\text{inj}}\\) is defined as from HH_helper_functions import syn_current I , t_on , t_off , dt , t , A_soma = syn_current () The Hodgkin-Huxley simulator is given by: from HH_helper_functions import HHsimulator Putting the input current and the simulator together: def run_HH_model ( params ): params = np . asarray ( params ) # input current, time step I , t_on , t_off , dt , t , A_soma = syn_current () t = np . arange ( 0 , len ( I ), 1 ) * dt # initial voltage V0 = - 70 states = HHsimulator ( V0 , params . reshape ( 1 , - 1 ), dt , t , I ) return dict ( data = states . reshape ( - 1 ), time = t , dt = dt , I = I . reshape ( - 1 )) To get an idea of the output of the Hodgkin-Huxley model, let us generate some voltage traces for different parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), given the input current \\(I_{\\text{inj}}\\) : # three sets of (g_Na, g_K) params = np . array ([[ 50.0 , 1.0 ], [ 4.0 , 1.5 ], [ 20.0 , 15.0 ]]) num_samples = len ( params [:, 0 ]) sim_samples = np . zeros (( num_samples , len ( I ))) for i in range ( num_samples ): sim_samples [ i , :] = run_HH_model ( params = params [ i , :])[ \"data\" ] # colors for traces col_min = 2 num_colors = num_samples + col_min cm1 = mpl . cm . Blues col1 = [ cm1 ( 1.0 * i / num_colors ) for i in range ( col_min , num_colors )] fig = plt . figure ( figsize = ( 7 , 5 )) gs = mpl . gridspec . GridSpec ( 2 , 1 , height_ratios = [ 4 , 1 ]) ax = plt . subplot ( gs [ 0 ]) for i in range ( num_samples ): plt . plot ( t , sim_samples [ i , :], color = col1 [ i ], lw = 2 ) plt . ylabel ( \"voltage (mV)\" ) ax . set_xticks ([]) ax . set_yticks ([ - 80 , - 20 , 40 ]) ax = plt . subplot ( gs [ 1 ]) plt . plot ( t , I * A_soma * 1e3 , \"k\" , lw = 2 ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"input (nA)\" ) ax . set_xticks ([ 0 , max ( t ) / 2 , max ( t )]) ax . set_yticks ([ 0 , 1.1 * np . max ( I * A_soma * 1e3 )]) ax . yaxis . set_major_formatter ( mpl . ticker . FormatStrFormatter ( \" %.2f \" )) plt . show () As can be seen, the voltage traces can be quite diverse for different parameter values. Often, we are not interested in matching the exact trace, but only in matching certain features thereof. In this example of the Hodgkin-Huxley model, the summary features are the number of spikes, the mean resting potential, the standard deviation of the resting potential, and the first four voltage moments: mean, standard deviation, skewness and kurtosis. Using the function calculate_summary_statistics() imported below, we obtain these statistics from the output of the Hodgkin Huxley simulator. from HH_helper_functions import calculate_summary_statistics Lastly, we define a function that performs all of the above steps at once. The function simulation_wrapper takes in conductance values, runs the Hodgkin Huxley model and then returns the summary statistics. def simulation_wrapper ( params ): \"\"\" Returns summary statistics from conductance values in `params`. Summarizes the output of the HH simulator and converts it to `torch.Tensor`. \"\"\" obs = run_HH_model ( params ) summstats = torch . as_tensor ( calculate_summary_statistics ( obs )) return summstats sbi takes any function as simulator. Thus, sbi also has the flexibility to use simulators that utilize external packages, e.g., Brian ( http://briansimulator.org/ ), nest ( https://www.nest-simulator.org/ ), or NEURON ( https://neuron.yale.edu/neuron/ ). External simulators do not even need to be Python-based as long as they store simulation outputs in a format that can be read from Python. All that is necessary is to wrap your external simulator of choice into a Python callable that takes a parameter set and outputs a set of summary statistics we want to fit the parameters to.","title":"2. Simulator"},{"location":"examples/00_HH_simulator/#3-prior-over-model-parameters","text":"Now that we have the simulator, we need to define a function with the prior over the model parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), which in this case is chosen to be a Uniform distribution: prior_min = [ 0.5 , 1e-4 ] prior_max = [ 80.0 , 15.0 ] prior = utils . torchutils . BoxUniform ( low = torch . as_tensor ( prior_min ), high = torch . as_tensor ( prior_max ) )","title":"3. Prior over model parameters"},{"location":"examples/00_HH_simulator/#inference","text":"Now that we have all the required components, we can run inference with SNPE to identify parameters whose activity matches this trace. posterior = infer ( simulation_wrapper , prior , method = \"SNPE\" , num_simulations = 300 , num_workers = 4 ) HBox(children=(FloatProgress(value=0.0, description='Running 300 simulations in 300 batches.', max=300.0, styl\u2026 Neural network successfully converged after 233 epochs. Note sbi can parallelize your simulator. If you experience problems with parallelization, try setting num_workers=1 and please give us an error report as a GitHub issue .","title":"Inference"},{"location":"examples/00_HH_simulator/#coming-back-to-the-observed-data","text":"As mentioned at the beginning of the tutorial, the observed data are generated by the Hodgkin-Huxley model with a set of known parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ). To illustrate how to compute the summary statistics of the observed data, let us regenerate the observed data: # true parameters and respective labels true_params = np . array ([ 50.0 , 5.0 ]) labels_params = [ r \"$g_ {Na} $\" , r \"$g_ {K} $\" ] observation_trace = run_HH_model ( true_params ) observation_summary_statistics = calculate_summary_statistics ( observation_trace ) As we already shown above, the observed voltage traces look as follows: fig = plt . figure ( figsize = ( 7 , 5 )) gs = mpl . gridspec . GridSpec ( 2 , 1 , height_ratios = [ 4 , 1 ]) ax = plt . subplot ( gs [ 0 ]) plt . plot ( observation_trace [ \"time\" ], observation_trace [ \"data\" ]) plt . ylabel ( \"voltage (mV)\" ) plt . title ( \"observed data\" ) plt . setp ( ax , xticks = [], yticks = [ - 80 , - 20 , 40 ]) ax = plt . subplot ( gs [ 1 ]) plt . plot ( observation_trace [ \"time\" ], I * A_soma * 1e3 , \"k\" , lw = 2 ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"input (nA)\" ) ax . set_xticks ([ 0 , max ( observation_trace [ \"time\" ]) / 2 , max ( observation_trace [ \"time\" ])]) ax . set_yticks ([ 0 , 1.1 * np . max ( I * A_soma * 1e3 )]) ax . yaxis . set_major_formatter ( mpl . ticker . FormatStrFormatter ( \" %.2f \" ))","title":"Coming back to the observed data"},{"location":"examples/00_HH_simulator/#analysis-of-the-posterior-given-the-observed-data","text":"After running the inference algorithm, let us inspect the inferred posterior distribution over the parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), given the observed trace. To do so, we first draw samples (i.e. consistent parameter sets) from the posterior: samples = posterior . sample (( 10000 ,), x = observation_summary_statistics ) HBox(children=(FloatProgress(value=0.0, description='Drawing 10000 posterior samples', max=10000.0, style=Prog\u2026 fig , axes = analysis . pairplot ( samples , limits = [[ 0.5 , 80 ], [ 1e-4 , 15.0 ]], ticks = [[ 0.5 , 80 ], [ 1e-4 , 15.0 ]], figsize = ( 5 , 5 ), points = true_params , points_offdiag = { \"markersize\" : 6 }, points_colors = \"r\" , ); As can be seen, the inferred posterior contains the ground-truth parameters (red) in a high-probability region. Now, let us sample parameters from the posterior distribution, simulate the Hodgkin-Huxley model for this parameter set and compare the simulations with the observed data: # Draw a sample from the posterior and convert to numpy for plotting. posterior_sample = posterior . sample (( 1 ,), x = observation_summary_statistics ) . numpy () HBox(children=(FloatProgress(value=0.0, description='Drawing 1 posterior samples', max=1.0, style=ProgressStyl\u2026 fig = plt . figure ( figsize = ( 7 , 5 )) # plot observation t = observation_trace [ \"time\" ] y_obs = observation_trace [ \"data\" ] plt . plot ( t , y_obs , lw = 2 , label = \"observation\" ) # simulate and plot samples x = run_HH_model ( posterior_sample ) plt . plot ( t , x [ \"data\" ], \"--\" , lw = 2 , label = \"posterior sample\" ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"voltage (mV)\" ) ax = plt . gca () handles , labels = ax . get_legend_handles_labels () ax . legend ( handles [:: - 1 ], labels [:: - 1 ], bbox_to_anchor = ( 1.3 , 1 ), loc = \"upper right\" ) ax . set_xticks ([ 0 , 60 , 120 ]) ax . set_yticks ([ - 80 , - 20 , 40 ]); As can be seen, the sample from the inferred posterior leads to simulations that closely resemble the observed data, confirming that SNPE did a good job at capturing the observed data in this simple case.","title":"Analysis of the posterior given the observed data"},{"location":"examples/00_HH_simulator/#references","text":"A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117(4):500\u2013544, 1952. M. Pospischil, M. Toledo-Rodriguez, C. Monier, Z. Piwkowska, T. Bal, Y. Fr\u00e9gnac, H. Markram, and A. Destexhe. Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biological Cybernetics, 99(4-5), 2008.","title":"References"},{"location":"examples/01_decision_making_model/","text":"SBI for decision-making models \u00b6 In a previous tutorial , we showed how to use SBI with trial-based iid data. Such scenarios can arise, for example, in models of perceptual decision making. In addition to trial-based iid data points, these models often come with mixed data types and varying experimental conditions. Here, we show how sbi can be used to perform inference in such models with the MNLE method. Trial-based SBI with mixed data types \u00b6 In some cases, models with trial-based data additionally return data with mixed data types, e.g., continous and discrete data. For example, most computational models of decision-making have continuous reaction times and discrete choices as output. This can induce a problem when performing trial-based SBI that relies on learning a neural likelihood: It is challenging for most density estimators to handle both, continuous and discrete data at the same time. However, there is a recent SBI method for solving this problem, it\u2019s called Mixed Neural Likelihood Estimation (MNLE). It works just like NLE, but with mixed data types. The trick is that it learns two separate density estimators, one for the discrete part of the data, and one for the continuous part, and combines the two to obtain the final neural likelihood. Crucially, the continuous density estimator is trained conditioned on the output of the discrete one, such that statistical dependencies between the discrete and continuous data (e.g., between choices and reaction times) are modeled as well. The interested reader is referred to the original paper available here . MNLE was recently added to sbi (see this PR and also issue ) and follows the same API as SNLE . In this tutorial we will show how to apply MNLE to mixed data, and how to deal with varying experimental conditions. Toy problem for MNLE \u00b6 To illustrate MNLE we set up a toy simulator that outputs mixed data and for which we know the likelihood such we can obtain reference posterior samples via MCMC. Simulator : To simulate mixed data we do the following Sample reaction time from inverse Gamma Sample choices from Binomial Return reaction time \\(rt \\in (0, \\infty)\\) and choice index \\(c \\in \\{0, 1\\}\\) \\[ c \\sim \\text{Binomial}(\\rho) \\\\ rt \\sim \\text{InverseGamma}(\\alpha=2, \\beta) \\\\ \\] Prior : The priors of the two parameters \\(\\rho\\) and \\(\\beta\\) are independent. We define a Beta prior over the probabilty parameter of the Binomial used in the simulator and a Gamma prior over the shape-parameter of the inverse Gamma used in the simulator: \\[ p(\\beta, \\rho) = p(\\beta) \\; p(\\rho) ; \\\\ p(\\beta) = \\text{Gamma}(1, 0.5) \\\\ p(\\text{probs}) = \\text{Beta}(2, 2) \\] Because the InverseGamma and the Binomial likelihoods are well-defined we can perform MCMC on this problem and obtain reference-posterior samples. import matplotlib.pyplot as plt import torch from torch import Tensor from sbi.inference import MNLE from pyro.distributions import InverseGamma from torch.distributions import Beta , Binomial , Categorical , Gamma from sbi.utils import MultipleIndependent from sbi.utils.metrics import c2st from sbi.analysis import pairplot from sbi.inference import MCMCPosterior from sbi.utils.torchutils import atleast_2d from sbi.inference.potentials.likelihood_based_potential import ( MixedLikelihoodBasedPotential , ) from sbi.utils.conditional_density_utils import ConditionedPotential from sbi.utils import mcmc_transform from sbi.inference.potentials.base_potential import BasePotential # Toy simulator for mixed data def mixed_simulator ( theta : Tensor , concentration_scaling : float = 1.0 ): \"\"\"Returns a sample from a mixed distribution given parameters theta. Args: theta: batch of parameters, shape (batch_size, 2) concentration_scaling: scaling factor for the concentration parameter of the InverseGamma distribution, mimics an experimental condition. \"\"\" beta , ps = theta [:, : 1 ], theta [:, 1 :] choices = Binomial ( probs = ps ) . sample () rts = InverseGamma ( concentration = concentration_scaling * torch . ones_like ( beta ), rate = beta ) . sample () return torch . cat (( rts , choices ), dim = 1 ) # The potential function defines the ground truth likelihood and allows us to obtain reference posterior samples via MCMC. class PotentialFunctionProvider ( BasePotential ): allow_iid_x = True # type: ignore def __init__ ( self , prior , x_o , concentration_scaling = 1.0 , device = \"cpu\" ): super () . __init__ ( prior , x_o , device ) self . concentration_scaling = concentration_scaling def __call__ ( self , theta , track_gradients : bool = True ): theta = atleast_2d ( theta ) with torch . set_grad_enabled ( track_gradients ): iid_ll = self . iid_likelihood ( theta ) return iid_ll + self . prior . log_prob ( theta ) def iid_likelihood ( self , theta ): lp_choices = torch . stack ( [ Binomial ( probs = th . reshape ( 1 , - 1 )) . log_prob ( self . x_o [:, 1 :]) for th in theta [:, 1 :] ], dim = 1 , ) lp_rts = torch . stack ( [ InverseGamma ( concentration = self . concentration_scaling * torch . ones_like ( beta_i ), rate = beta_i , ) . log_prob ( self . x_o [:, : 1 ]) for beta_i in theta [:, : 1 ] ], dim = 1 , ) joint_likelihood = ( lp_choices + lp_rts ) . squeeze () assert joint_likelihood . shape == torch . Size ([ self . x_o . shape [ 0 ], theta . shape [ 0 ]]) return joint_likelihood . sum ( 0 ) # Define independent prior. prior = MultipleIndependent ( [ Gamma ( torch . tensor ([ 1.0 ]), torch . tensor ([ 0.5 ])), Beta ( torch . tensor ([ 2.0 ]), torch . tensor ([ 2.0 ])), ], validate_args = False , ) Obtain reference-posterior samples via analytical likelihood and MCMC \u00b6 torch . manual_seed ( 42 ) num_trials = 10 num_samples = 1000 theta_o = prior . sample (( 1 ,)) x_o = mixed_simulator ( theta_o . repeat ( num_trials , 1 )) mcmc_kwargs = dict ( num_chains = 20 , warmup_steps = 50 , method = \"slice_np_vectorized\" , init_strategy = \"proposal\" , ) true_posterior = MCMCPosterior ( potential_fn = PotentialFunctionProvider ( prior , x_o ), proposal = prior , theta_transform = mcmc_transform ( prior , enable_transform = True ), ** mcmc_kwargs , ) true_samples = true_posterior . sample (( num_samples ,)) /Users/janbolts/qode/sbi/sbi/utils/sbiutils.py:342: UserWarning: An x with a batch size of 10 was passed. It will be interpreted as a batch of independent and identically distributed data X={x_1, ..., x_n}, i.e., data generated based on the same underlying (unknown) parameter. The resulting posterior will be with respect to entire batch, i.e,. p(theta | X). warnings.warn( Running vectorized MCMC with 20 chains: 0%| | 0/20000 [00:00 1 , you might experience an error that a certain object from your simulator could not be pickled (an example can be found here ). This can be fixed by forcing sbi to pickle with dill instead of the default cloudpickle . To do so, adjust your code as follows: Install dill : pip install dill At the very beginning of your python script, set the pickler to dill : from joblib.externals.loky import set_loky_pickler set_loky_pickler ( \"dill\" ) Move all imports required by your simulator into the simulator: # Imports specified outside of the simulator will break dill: import torch def my_simulator ( parameters ): return torch . ones ( 1 , 10 ) # Therefore, move the imports into the simulator: def my_simulator ( parameters ): import torch return torch . ones ( 1 , 10 ) Alternative: parallelize yourself \u00b6 You can also write your own code to parallelize simulations with whatever multiprocessing framework you prefer. You can then simulate your data outside of sbi and pass the simulated data as shown in the flexible interface : Some more background \u00b6 sbi uses joblib to parallelize simulations, which in turn uses pickle or cloudpickle to serialize the simulator. Almost all simulators will be picklable with cloudpickle , but we have experienced issues e.g. with neuron simulators, see here .","title":"When using multiple workers, I get a pickling error. Can I still use multiprocessing?"},{"location":"faq/question_03/#when-using-multiple-workers-i-get-a-pickling-error-can-i-still-use-multiprocessing","text":"Yes, but you will have to make a few adjustments to your code. Some background: when using num_workers > 1 , you might experience an error that a certain object from your simulator could not be pickled (an example can be found here ). This can be fixed by forcing sbi to pickle with dill instead of the default cloudpickle . To do so, adjust your code as follows: Install dill : pip install dill At the very beginning of your python script, set the pickler to dill : from joblib.externals.loky import set_loky_pickler set_loky_pickler ( \"dill\" ) Move all imports required by your simulator into the simulator: # Imports specified outside of the simulator will break dill: import torch def my_simulator ( parameters ): return torch . ones ( 1 , 10 ) # Therefore, move the imports into the simulator: def my_simulator ( parameters ): import torch return torch . ones ( 1 , 10 )","title":"When using multiple workers, I get a pickling error. Can I still use multiprocessing?"},{"location":"faq/question_03/#alternative-parallelize-yourself","text":"You can also write your own code to parallelize simulations with whatever multiprocessing framework you prefer. You can then simulate your data outside of sbi and pass the simulated data as shown in the flexible interface :","title":"Alternative: parallelize yourself"},{"location":"faq/question_03/#some-more-background","text":"sbi uses joblib to parallelize simulations, which in turn uses pickle or cloudpickle to serialize the simulator. Almost all simulators will be picklable with cloudpickle , but we have experienced issues e.g. with neuron simulators, see here .","title":"Some more background"},{"location":"faq/question_04/","text":"Can I use the GPU for training the density estimator? \u00b6 TLDR; Yes, by passing device=\"cuda\" and by passing a prior that lives on the device name your passed. But no speed-ups for default density estimators. Yes. When creating the inference object in the flexible interface, you can pass the device as an argument, e.g., inference = SNPE ( prior , device = \"cuda\" , density_estimator = \"maf\" ) The device is set to \"cpu\" by default, and it can be set to anything, as long as it maps to an existing PyTorch CUDA device. sbi will take care of copying the net and the training data to and from the device . Note that the prior must be on the training device already, e.g., when passing device=\"cuda:0\" , make sure to pass a prior object that was created on that device, e.g., prior = torch.distributions.MultivariateNormal(loc=torch.zeros(2, device=\"cuda:0\"), covariance_matrix=torch.eye(2, device=\"cuda:0\")) . Performance \u00b6 Whether or not you reduce your training time when training on a GPU depends on the problem at hand. We provide a couple of default density estimators for SNPE , SNLE and SNRE , e.g., a mixture density network ( density_estimator=\"mdn\" ) or a Masked Autoregressive Flow ( density_estimator=\"maf\" ). For those default density estimators we do not expect a speed up. This is because the underlying neural networks are quite shallow and not tall, e.g., they do not have many parameters or matrix operations that profit a lot from being executed on the GPU. A speed up through training on the GPU will most likely become visible when you are using convolutional modules in your neural networks. E.g., when passing an embedding net for image processing like in this example: https://github.com/sbi-dev/sbi/blob/main/tutorials/05_embedding_net.ipynb .","title":"Can I use the GPU for training the density estimator?"},{"location":"faq/question_04/#can-i-use-the-gpu-for-training-the-density-estimator","text":"TLDR; Yes, by passing device=\"cuda\" and by passing a prior that lives on the device name your passed. But no speed-ups for default density estimators. Yes. When creating the inference object in the flexible interface, you can pass the device as an argument, e.g., inference = SNPE ( prior , device = \"cuda\" , density_estimator = \"maf\" ) The device is set to \"cpu\" by default, and it can be set to anything, as long as it maps to an existing PyTorch CUDA device. sbi will take care of copying the net and the training data to and from the device . Note that the prior must be on the training device already, e.g., when passing device=\"cuda:0\" , make sure to pass a prior object that was created on that device, e.g., prior = torch.distributions.MultivariateNormal(loc=torch.zeros(2, device=\"cuda:0\"), covariance_matrix=torch.eye(2, device=\"cuda:0\")) .","title":"Can I use the GPU for training the density estimator?"},{"location":"faq/question_04/#performance","text":"Whether or not you reduce your training time when training on a GPU depends on the problem at hand. We provide a couple of default density estimators for SNPE , SNLE and SNRE , e.g., a mixture density network ( density_estimator=\"mdn\" ) or a Masked Autoregressive Flow ( density_estimator=\"maf\" ). For those default density estimators we do not expect a speed up. This is because the underlying neural networks are quite shallow and not tall, e.g., they do not have many parameters or matrix operations that profit a lot from being executed on the GPU. A speed up through training on the GPU will most likely become visible when you are using convolutional modules in your neural networks. E.g., when passing an embedding net for image processing like in this example: https://github.com/sbi-dev/sbi/blob/main/tutorials/05_embedding_net.ipynb .","title":"Performance"},{"location":"faq/question_05/","text":"How should I save and load objects in sbi ? \u00b6 NeuralPosterior objects are picklable. import pickle # ... run inference posterior = inference . build_posterior () with open ( \"/path/to/my_posterior.pkl\" , \"wb\" ) as handle : pickle . dump ( posterior , handle ) Note: posterior objects that were saved under sbi v0.17.2 or older can not be loaded under sbi v0.18.0 or newer. Note: if you try to load a posterior that was saved under sbi v0.14.x or earlier under sbi v0.15.x until sbi v0.17.x , you have to add: import sys from sbi.utils import user_input_checks_utils sys . modules [ \"sbi.user_input.user_input_checks_utils\" ] = user_input_checks_utils to your script before loading the posterior. As of sbi v0.18.0 , NeuralInference objects are also picklable. import pickle # ... run inference posterior = inference . build_posterior () with open ( \"/path/to/my_inference.pkl\" , \"wb\" ) as handle : pickle . dump ( inference , handle ) However, saving and loading the inference object will slightly modify the object (in order to make it serializable). These modifications lead to the following two changes in behaviour: 1) Retraining from scratch is not supported, i.e. .train(..., retrain_from_scratch=True) does not work. 2) When the loaded object calls the .train() method, it generates a new tensorboard summary writer (instead of appending to the current one).","title":"How should I save and load objects in sbi?"},{"location":"faq/question_05/#how-should-i-save-and-load-objects-in-sbi","text":"NeuralPosterior objects are picklable. import pickle # ... run inference posterior = inference . build_posterior () with open ( \"/path/to/my_posterior.pkl\" , \"wb\" ) as handle : pickle . dump ( posterior , handle ) Note: posterior objects that were saved under sbi v0.17.2 or older can not be loaded under sbi v0.18.0 or newer. Note: if you try to load a posterior that was saved under sbi v0.14.x or earlier under sbi v0.15.x until sbi v0.17.x , you have to add: import sys from sbi.utils import user_input_checks_utils sys . modules [ \"sbi.user_input.user_input_checks_utils\" ] = user_input_checks_utils to your script before loading the posterior. As of sbi v0.18.0 , NeuralInference objects are also picklable. import pickle # ... run inference posterior = inference . build_posterior () with open ( \"/path/to/my_inference.pkl\" , \"wb\" ) as handle : pickle . dump ( inference , handle ) However, saving and loading the inference object will slightly modify the object (in order to make it serializable). These modifications lead to the following two changes in behaviour: 1) Retraining from scratch is not supported, i.e. .train(..., retrain_from_scratch=True) does not work. 2) When the loaded object calls the .train() method, it generates a new tensorboard summary writer (instead of appending to the current one).","title":"How should I save and load objects in sbi?"},{"location":"faq/question_06/","text":"Can I stop neural network training and resume it later? \u00b6 Many clusters have a time limit and sbi might exceed this limit. You can circumvent this problem by using the flexible interface . After simulations are finished, sbi trains a neural network. If this process takes too long, you can stop training and resume it later. The syntax is: inference = SNPE ( prior = prior ) inference = inference . append_simulations ( theta , x ) inference . train ( max_num_epochs = 300 ) # Pick `max_num_epochs` such that it does not exceed the runtime. with open ( \"path/to/my/inference.pkl\" , \"wb\" ) as handle : pickle . dump ( inference , handle ) # To resume training: with open ( \"path/to/my/inference.pkl\" , \"rb\" ) as handle : inference_from_disk = pickle . load ( handle ) inference_from_disk . train ( resume_training = True , max_num_epochs = 600 ) # Run epochs 301 until 600 (or stop early). posterior = inference_from_disk . build_posterior ()","title":"Can I stop neural network training and resume it later?"},{"location":"faq/question_06/#can-i-stop-neural-network-training-and-resume-it-later","text":"Many clusters have a time limit and sbi might exceed this limit. You can circumvent this problem by using the flexible interface . After simulations are finished, sbi trains a neural network. If this process takes too long, you can stop training and resume it later. The syntax is: inference = SNPE ( prior = prior ) inference = inference . append_simulations ( theta , x ) inference . train ( max_num_epochs = 300 ) # Pick `max_num_epochs` such that it does not exceed the runtime. with open ( \"path/to/my/inference.pkl\" , \"wb\" ) as handle : pickle . dump ( inference , handle ) # To resume training: with open ( \"path/to/my/inference.pkl\" , \"rb\" ) as handle : inference_from_disk = pickle . load ( handle ) inference_from_disk . train ( resume_training = True , max_num_epochs = 600 ) # Run epochs 301 until 600 (or stop early). posterior = inference_from_disk . build_posterior ()","title":"Can I stop neural network training and resume it later?"},{"location":"faq/question_07/","text":"Can I use a custom prior with sbi? \u00b6 sbi works with torch distributions only so we recommend to use those whenever possible. For example, when you are used to using scipy.stats distributions as priors then we recommend using the corresponding torch.distributions , most common distributions are implemented there. In case you want to use a custom prior that is not in the set of common distributions that\u2019s possible as well: You need to write a prior class that mimicks the behaviour of a torch.distributions.Distribution class. Then sbi will wrap this class to make it a fully functional torch Distribution . Essentially, the class needs two methods: .sample(sample_shape) , where sample_shape is a shape tuple, e.g., (n,) , and returns a batch of n samples, e.g., of shape (n, 2)` for a two dimenional prior. .log_prob(value) method that returns the \u201clog probs\u201d of parameters under the prior, e.g., for a batches of n parameters with shape (n, ndims) it should return a log probs array of shape (n,) . For sbi > 0.17.2 this could look like the following: class CustomUniformPrior : \"\"\"User defined numpy uniform prior. Custom prior with user-defined valid .sample and .log_prob methods. \"\"\" def __init__ ( self , lower : Tensor , upper : Tensor , return_numpy : bool = False ): self . lower = lower self . upper = upper self . dist = BoxUniform ( lower , upper ) self . return_numpy = return_numpy def sample ( self , sample_shape = torch . Size ([])): samples = self . dist . sample ( sample_shape ) return samples . numpy () if self . return_numpy else samples def log_prob ( self , values ): if self . return_numpy : values = torch . as_tensor ( values ) log_probs = self . dist . log_prob ( values ) return log_probs . numpy () if self . return_numpy else log_probs Once you have such a class you can wrap into a Distribution using the process_prior function sbi provides: from sbi.utils import process_prior custom_prior = CustomUniformPrior ( torch . zeros ( 2 ), torch . ones ( 2 )) prior , * _ = process_prior ( custom_prior ) # Keeping only the first return. # use this wrapped prior in sbi... In sbi it is sometimes required to check the support of the prior, e.g., when the prior support is bounded and one wants to reject samples from the posterior density estimator that lie outside the prior support. In torch Distributions this is handled automatically, however, when using a custom prior it is not. Thus, if your prior has bounded support (like the one above) it makes sense to pass the bounds to the wrapper function such that sbi can pass them to torch Distributions : from sbi.utils import process_prior custom_prior = CustomUniformPrior ( torch . zeros ( 2 ), torch . ones ( 2 )) prior = process_prior ( custom_prior , custom_prior_wrapper_kwargs = dict ( lower_bound = torch . zeros ( 2 ), upper_bound = torch . ones ( 2 ))) # use this wrapped prior in sbi... Note that in custom_prior_wrapper_kwargs you can pass additinal arguments for the wrapper, e.g., validate_args or arg_constraints see the Distribution documentation for more details. If you are running sbi < 0.17.2 and use SNLE the code above will produce a NotImplementedError (see #581 ). In this case you need to update to a newer version of sbi or use SNPE instead.","title":"Can I use a custom prior with sbi?"},{"location":"faq/question_07/#can-i-use-a-custom-prior-with-sbi","text":"sbi works with torch distributions only so we recommend to use those whenever possible. For example, when you are used to using scipy.stats distributions as priors then we recommend using the corresponding torch.distributions , most common distributions are implemented there. In case you want to use a custom prior that is not in the set of common distributions that\u2019s possible as well: You need to write a prior class that mimicks the behaviour of a torch.distributions.Distribution class. Then sbi will wrap this class to make it a fully functional torch Distribution . Essentially, the class needs two methods: .sample(sample_shape) , where sample_shape is a shape tuple, e.g., (n,) , and returns a batch of n samples, e.g., of shape (n, 2)` for a two dimenional prior. .log_prob(value) method that returns the \u201clog probs\u201d of parameters under the prior, e.g., for a batches of n parameters with shape (n, ndims) it should return a log probs array of shape (n,) . For sbi > 0.17.2 this could look like the following: class CustomUniformPrior : \"\"\"User defined numpy uniform prior. Custom prior with user-defined valid .sample and .log_prob methods. \"\"\" def __init__ ( self , lower : Tensor , upper : Tensor , return_numpy : bool = False ): self . lower = lower self . upper = upper self . dist = BoxUniform ( lower , upper ) self . return_numpy = return_numpy def sample ( self , sample_shape = torch . Size ([])): samples = self . dist . sample ( sample_shape ) return samples . numpy () if self . return_numpy else samples def log_prob ( self , values ): if self . return_numpy : values = torch . as_tensor ( values ) log_probs = self . dist . log_prob ( values ) return log_probs . numpy () if self . return_numpy else log_probs Once you have such a class you can wrap into a Distribution using the process_prior function sbi provides: from sbi.utils import process_prior custom_prior = CustomUniformPrior ( torch . zeros ( 2 ), torch . ones ( 2 )) prior , * _ = process_prior ( custom_prior ) # Keeping only the first return. # use this wrapped prior in sbi... In sbi it is sometimes required to check the support of the prior, e.g., when the prior support is bounded and one wants to reject samples from the posterior density estimator that lie outside the prior support. In torch Distributions this is handled automatically, however, when using a custom prior it is not. Thus, if your prior has bounded support (like the one above) it makes sense to pass the bounds to the wrapper function such that sbi can pass them to torch Distributions : from sbi.utils import process_prior custom_prior = CustomUniformPrior ( torch . zeros ( 2 ), torch . ones ( 2 )) prior = process_prior ( custom_prior , custom_prior_wrapper_kwargs = dict ( lower_bound = torch . zeros ( 2 ), upper_bound = torch . ones ( 2 ))) # use this wrapped prior in sbi... Note that in custom_prior_wrapper_kwargs you can pass additinal arguments for the wrapper, e.g., validate_args or arg_constraints see the Distribution documentation for more details. If you are running sbi < 0.17.2 and use SNLE the code above will produce a NotImplementedError (see #581 ). In this case you need to update to a newer version of sbi or use SNPE instead.","title":"Can I use a custom prior with sbi?"},{"location":"tutorial/00_getting_started/","text":"Getting started with sbi \u00b6 Note, you can find the original version of this notebook at https://github.com/sbi-dev/sbi/blob/main/tutorials/00_getting_started.ipynb in the sbi repository. import torch from sbi import utils as utils from sbi import analysis as analysis from sbi.inference.base import infer Running the inference procedure \u00b6 sbi provides a simple interface to run state-of-the-art algorithms for simulation-based inference. For inference, you need to provide two ingredients: 1) a prior distribution that allows to sample parameter sets. 2) a simulator that takes parameter sets and produces simulation outputs. For example, we can have a 3-dimensional parameter space with a uniform prior between [-1,1] and a simple simulator that for the sake of example adds 1.0 and some Gaussian noise to the parameter set: num_dim = 3 prior = utils . BoxUniform ( low =- 2 * torch . ones ( num_dim ), high = 2 * torch . ones ( num_dim )) def simulator ( parameter_set ): return 1.0 + parameter_set + torch . randn ( parameter_set . shape ) * 0.1 sbi can then run inference: posterior = infer ( simulator , prior , method = \"SNPE\" , num_simulations = 1000 ) Running 1000 simulations.: 0%| | 0/1000 [00:001 , the posterior is no longer amortized: it will give good results when sampled around x=observation , but possibly bad results for other x . Once we have obtained the posterior, we can .sample() , .log_prob() , or .pairplot() in the same way as for the simple interface. posterior_samples = posterior . sample (( 10000 ,), x = x_o ) # plot posterior samples _ = analysis . pairplot ( posterior_samples , limits = [[ - 2 , 2 ], [ - 2 , 2 ], [ - 2 , 2 ]], figsize = ( 5 , 5 ) ) Drawing 10000 posterior samples: 0%| | 0/10000 [00:001 , the posterior is no longer amortized: it will give good results when sampled around x=observation , but possibly bad results for other x . Once we have obtained the posterior, we can .sample() , .log_prob() , or .pairplot() in the same way as for the simple interface. posterior_samples = posterior . sample (( 10000 ,), x = x_o ) # plot posterior samples _ = analysis . pairplot ( posterior_samples , limits = [[ - 2 , 2 ], [ - 2 , 2 ], [ - 2 , 2 ]], figsize = ( 5 , 5 ) ) Drawing 10000 posterior samples: 0%| | 0/10000 [00:00 The simulator model \u00b6 The simulator model that we consider has two parameters: \\(r\\) and \\(\\theta\\) . On each run, it generates 100 two-dimensional points centered around \\((r \\cos(\\theta), r \\sin(\\theta))\\) and perturbed by a Gaussian noise with variance 0.01. Instead of simply outputting the \\((x,y)\\) coordinates of each data point, the model generates a grayscale image of the scattered points with dimensions 32 by 32. This image is further perturbed by an uniform noise with values betweeen 0 and 0.2. The code below defines such model. def simulator_model ( parameter , return_points = False ): \"\"\"Simulator model with two-dimensional input parameter and 1024-dimensional output This simulator serves as a basic example for using a neural net for learning summary features. It has only two input parameters but generates high-dimensional output vectors. The data is generated as follows: (-) Input: parameter = [r, theta] (1) Generate 100 two-dimensional points centered around (r cos(theta),r sin(theta)) and perturbed by a Gaussian noise with variance 0.01 (2) Create a grayscale image I of the scattered points with dimensions 32 by 32 (3) Perturb I with an uniform noise with values betweeen 0 and 0.2 (-) Output: I Parameters ---------- parameter : array-like, shape (2) The two input parameters of the model, ordered as [r, theta] return_points : bool (default: False) Whether the simulator should return the coordinates of the simulated data points as well Returns ------- I: torch tensor, shape (1, 1024) Output flattened image (optional) points: array-like, shape (100, 2) Coordinates of the 2D simulated data points \"\"\" r = parameter [ 0 ] theta = parameter [ 1 ] sigma_points = 0.10 npoints = 100 points = [] for _ in range ( npoints ): x = r * torch . cos ( theta ) + sigma_points * torch . randn ( 1 ) y = r * torch . sin ( theta ) + sigma_points * torch . randn ( 1 ) points . append ([ x , y ]) points = torch . as_tensor ( points ) nx = 32 ny = 32 sigma_image = 0.20 I = torch . zeros ( nx , ny ) for point in points : pi = int (( point [ 0 ] - ( - 1 )) / (( + 1 ) - ( - 1 )) * nx ) pj = int (( point [ 1 ] - ( - 1 )) / (( + 1 ) - ( - 1 )) * ny ) if ( pi < nx ) and ( pj < ny ): I [ pi , pj ] = 1 I = I + sigma_image * torch . rand ( nx , ny ) I = I . T I = I . reshape ( 1 , - 1 ) if return_points : return I , points else : return I The figure below shows an example of the output of the simulator when \\(r = 0.70\\) and \\(\\theta = \\pi/4\\) # simulate samples true_parameter = torch . tensor ([ 0.70 , torch . pi / 4 ]) x_observed , x_points = simulator_model ( true_parameter , return_points = True ) # plot the observation fig , ax = plt . subplots ( facecolor = \"white\" , figsize = ( 11.15 , 5.61 ), ncols = 2 , constrained_layout = True ) circle = plt . Circle (( 0 , 0 ), 1.0 , color = \"k\" , ls = \"--\" , lw = 0.8 , fill = False ) ax [ 0 ] . add_artist ( circle ) ax [ 0 ] . scatter ( x_points [:, 0 ], x_points [:, 1 ], s = 20 ) ax [ 0 ] . set_xlabel ( \"x\" ) ax [ 0 ] . set_ylabel ( \"y\" ) ax [ 0 ] . set_xlim ( - 1 , + 1 ) ax [ 0 ] . set_xticks ([ - 1 , 0.0 , + 1.0 ]) ax [ 0 ] . set_ylim ( - 1 , + 1 ) ax [ 0 ] . set_yticks ([ - 1 , 0.0 , + 1.0 ]) ax [ 0 ] . set_title ( r \"original simulated points with $r = 0.70$ and $\\theta = \\pi/4$\" ) ax [ 1 ] . imshow ( x_observed . view ( 32 , 32 ), origin = \"lower\" , cmap = \"gray\" ) ax [ 1 ] . set_xticks ([]) ax [ 1 ] . set_yticks ([]) ax [ 1 ] . set_title ( \"noisy observed data (gray image with 32 x 32 pixels)\" ) Text(0.5, 1.0, 'noisy observed data (gray image with 32 x 32 pixels)') Defining an embedding_net \u00b6 An inference procedure applied to the output data from this simulator model determines the posterior distribution of \\(r\\) and \\(\\theta\\) given an observation of \\(x\\) , which lives in a 1024 dimensional space (32 x 32 = 1024). To avoid working directly on these high-dimensional vectors, one can use a convolutional neural network (CNN) that takes the 32x32 images as input and encodes them into 8-dimensional feature vectors. This CNN is trained along with the neural density estimator of the inference procedure and serves as an automatic summary statistics extractor. We define and instantiate the CNN as follows: class SummaryNet ( nn . Module ): def __init__ ( self ): super () . __init__ () # 2D convolutional layer self . conv1 = nn . Conv2d ( in_channels = 1 , out_channels = 6 , kernel_size = 5 , padding = 2 ) # Maxpool layer that reduces 32x32 image to 4x4 self . pool = nn . MaxPool2d ( kernel_size = 8 , stride = 8 ) # Fully connected layer taking as input the 6 flattened output arrays from the maxpooling layer self . fc = nn . Linear ( in_features = 6 * 4 * 4 , out_features = 8 ) def forward ( self , x ): x = x . view ( - 1 , 1 , 32 , 32 ) x = self . pool ( F . relu ( self . conv1 ( x ))) x = x . view ( - 1 , 6 * 4 * 4 ) x = F . relu ( self . fc ( x )) return x embedding_net = SummaryNet () The inference procedure \u00b6 With the embedding_net defined and instantiated, we can follow the usual workflow of an inference procedure in sbi . The embedding_net object appears as an input argument when instantiating the neural density estimator with utils.posterior_nn . # set prior distribution for the parameters prior = utils . BoxUniform ( low = torch . tensor ([ 0.0 , 0.0 ]), high = torch . tensor ([ 1.0 , 2 * torch . pi ]) ) # make a SBI-wrapper on the simulator object for compatibility simulator_wrapper , prior = prepare_for_sbi ( simulator_model , prior ) # instantiate the neural density estimator neural_posterior = utils . posterior_nn ( model = \"maf\" , embedding_net = embedding_net , hidden_features = 10 , num_transforms = 2 ) # setup the inference procedure with the SNPE-C procedure inference = SNPE ( prior = prior , density_estimator = neural_posterior ) # run the inference procedure on one round and 10000 simulated data points theta , x = simulate_for_sbi ( simulator_wrapper , prior , num_simulations = 10000 ) Running 10000 simulations.: 0%| | 0/10000 [00:00","title":"Learning summary statistics with a neural net"},{"location":"tutorial/05_embedding_net/#the-simulator-model","text":"The simulator model that we consider has two parameters: \\(r\\) and \\(\\theta\\) . On each run, it generates 100 two-dimensional points centered around \\((r \\cos(\\theta), r \\sin(\\theta))\\) and perturbed by a Gaussian noise with variance 0.01. Instead of simply outputting the \\((x,y)\\) coordinates of each data point, the model generates a grayscale image of the scattered points with dimensions 32 by 32. This image is further perturbed by an uniform noise with values betweeen 0 and 0.2. The code below defines such model. def simulator_model ( parameter , return_points = False ): \"\"\"Simulator model with two-dimensional input parameter and 1024-dimensional output This simulator serves as a basic example for using a neural net for learning summary features. It has only two input parameters but generates high-dimensional output vectors. The data is generated as follows: (-) Input: parameter = [r, theta] (1) Generate 100 two-dimensional points centered around (r cos(theta),r sin(theta)) and perturbed by a Gaussian noise with variance 0.01 (2) Create a grayscale image I of the scattered points with dimensions 32 by 32 (3) Perturb I with an uniform noise with values betweeen 0 and 0.2 (-) Output: I Parameters ---------- parameter : array-like, shape (2) The two input parameters of the model, ordered as [r, theta] return_points : bool (default: False) Whether the simulator should return the coordinates of the simulated data points as well Returns ------- I: torch tensor, shape (1, 1024) Output flattened image (optional) points: array-like, shape (100, 2) Coordinates of the 2D simulated data points \"\"\" r = parameter [ 0 ] theta = parameter [ 1 ] sigma_points = 0.10 npoints = 100 points = [] for _ in range ( npoints ): x = r * torch . cos ( theta ) + sigma_points * torch . randn ( 1 ) y = r * torch . sin ( theta ) + sigma_points * torch . randn ( 1 ) points . append ([ x , y ]) points = torch . as_tensor ( points ) nx = 32 ny = 32 sigma_image = 0.20 I = torch . zeros ( nx , ny ) for point in points : pi = int (( point [ 0 ] - ( - 1 )) / (( + 1 ) - ( - 1 )) * nx ) pj = int (( point [ 1 ] - ( - 1 )) / (( + 1 ) - ( - 1 )) * ny ) if ( pi < nx ) and ( pj < ny ): I [ pi , pj ] = 1 I = I + sigma_image * torch . rand ( nx , ny ) I = I . T I = I . reshape ( 1 , - 1 ) if return_points : return I , points else : return I The figure below shows an example of the output of the simulator when \\(r = 0.70\\) and \\(\\theta = \\pi/4\\) # simulate samples true_parameter = torch . tensor ([ 0.70 , torch . pi / 4 ]) x_observed , x_points = simulator_model ( true_parameter , return_points = True ) # plot the observation fig , ax = plt . subplots ( facecolor = \"white\" , figsize = ( 11.15 , 5.61 ), ncols = 2 , constrained_layout = True ) circle = plt . Circle (( 0 , 0 ), 1.0 , color = \"k\" , ls = \"--\" , lw = 0.8 , fill = False ) ax [ 0 ] . add_artist ( circle ) ax [ 0 ] . scatter ( x_points [:, 0 ], x_points [:, 1 ], s = 20 ) ax [ 0 ] . set_xlabel ( \"x\" ) ax [ 0 ] . set_ylabel ( \"y\" ) ax [ 0 ] . set_xlim ( - 1 , + 1 ) ax [ 0 ] . set_xticks ([ - 1 , 0.0 , + 1.0 ]) ax [ 0 ] . set_ylim ( - 1 , + 1 ) ax [ 0 ] . set_yticks ([ - 1 , 0.0 , + 1.0 ]) ax [ 0 ] . set_title ( r \"original simulated points with $r = 0.70$ and $\\theta = \\pi/4$\" ) ax [ 1 ] . imshow ( x_observed . view ( 32 , 32 ), origin = \"lower\" , cmap = \"gray\" ) ax [ 1 ] . set_xticks ([]) ax [ 1 ] . set_yticks ([]) ax [ 1 ] . set_title ( \"noisy observed data (gray image with 32 x 32 pixels)\" ) Text(0.5, 1.0, 'noisy observed data (gray image with 32 x 32 pixels)')","title":"The simulator model"},{"location":"tutorial/05_embedding_net/#defining-an-embedding_net","text":"An inference procedure applied to the output data from this simulator model determines the posterior distribution of \\(r\\) and \\(\\theta\\) given an observation of \\(x\\) , which lives in a 1024 dimensional space (32 x 32 = 1024). To avoid working directly on these high-dimensional vectors, one can use a convolutional neural network (CNN) that takes the 32x32 images as input and encodes them into 8-dimensional feature vectors. This CNN is trained along with the neural density estimator of the inference procedure and serves as an automatic summary statistics extractor. We define and instantiate the CNN as follows: class SummaryNet ( nn . Module ): def __init__ ( self ): super () . __init__ () # 2D convolutional layer self . conv1 = nn . Conv2d ( in_channels = 1 , out_channels = 6 , kernel_size = 5 , padding = 2 ) # Maxpool layer that reduces 32x32 image to 4x4 self . pool = nn . MaxPool2d ( kernel_size = 8 , stride = 8 ) # Fully connected layer taking as input the 6 flattened output arrays from the maxpooling layer self . fc = nn . Linear ( in_features = 6 * 4 * 4 , out_features = 8 ) def forward ( self , x ): x = x . view ( - 1 , 1 , 32 , 32 ) x = self . pool ( F . relu ( self . conv1 ( x ))) x = x . view ( - 1 , 6 * 4 * 4 ) x = F . relu ( self . fc ( x )) return x embedding_net = SummaryNet ()","title":"Defining an embedding_net"},{"location":"tutorial/05_embedding_net/#the-inference-procedure","text":"With the embedding_net defined and instantiated, we can follow the usual workflow of an inference procedure in sbi . The embedding_net object appears as an input argument when instantiating the neural density estimator with utils.posterior_nn . # set prior distribution for the parameters prior = utils . BoxUniform ( low = torch . tensor ([ 0.0 , 0.0 ]), high = torch . tensor ([ 1.0 , 2 * torch . pi ]) ) # make a SBI-wrapper on the simulator object for compatibility simulator_wrapper , prior = prepare_for_sbi ( simulator_model , prior ) # instantiate the neural density estimator neural_posterior = utils . posterior_nn ( model = \"maf\" , embedding_net = embedding_net , hidden_features = 10 , num_transforms = 2 ) # setup the inference procedure with the SNPE-C procedure inference = SNPE ( prior = prior , density_estimator = neural_posterior ) # run the inference procedure on one round and 10000 simulated data points theta , x = simulate_for_sbi ( simulator_wrapper , prior , num_simulations = 10000 ) Running 10000 simulations.: 0%| | 0/10000 [00:00] 1.3 Summary statistics \u00b6 We will compare two methods for defining summary statistics. One method uses three summary statistics which are function evaluations at three points in time. The other method uses a single summary statistic: the mean squared error between the observed and the simulated trace. In the second case, one then tries to obtain the posterior \\(p(\\theta | 0)\\) , i.e. the error being zero. These two methods are implemented below: \\(\\textbf{get_3_values()}\\) returns 3 function evaluations at \\(x=-0.5, x=0\\) and \\(x=0.75\\) . \\(\\textbf{get_MSE()}\\) returns the mean squared error between true and a quadratic function corresponding to a prior distributions sample. def get_3_values ( theta , seed = None ): \"\"\" Return 3 'y' values corresponding to x=-0.5,0,0.75 as summary statistic vector \"\"\" return np . array ( [ eval ( theta , - 0.5 , seed = seed ), eval ( theta , 0 , seed = seed ), eval ( theta , 0.75 , seed = seed ), ] ) . T def get_MSE ( theta , theta_o , seed = None ): \"\"\" Return the mean-squared error (MSE) i.e. Euclidean distance from the observation function \"\"\" _ , y = create_t_x ( theta_o , seed = seed ) # truth _ , y_ = create_t_x ( theta , seed = seed ) # simulations return np . mean ( np . square ( y_ - y ), axis = 0 , keepdims = True ) . T # MSE Let\u2019s try a couple of samples from our prior and see their summary statistics. Notice that these indeed change in small amounts every time you rerun it due to the noise, except if you set the seed. 1.4 Simulating data \u00b6 Let us see various plots of prior samples and their summary statistics versus the truth, i.e. our artificial observation. t , x_truth = create_t_x ( theta_o ) plt . plot ( t , x_truth , \"k\" , zorder = 1 , label = \"truth\" ) n_samples = 100 theta = prior . sample (( n_samples ,)) t , x = create_t_x ( theta . numpy ()) plt . plot ( t , x , \"grey\" , zorder = 0 ) plt . legend ()
- the simulator can simulate batches of parameters and return batches of data.
- the prior does not produce batches and samples and evaluates to Tensor.
- the output shape is a `torch.Size((1,N))` (i.e, has a leading batch dimension 1). If this is not possible, a suitable exception will be raised. Args: simulator: Simulator as provided by the user. prior: Prior as provided by the user. Returns: Tuple (simulator, prior) checked and matching the requirements of sbi. \"\"\" # Check prior, return PyTorch prior. prior , _ , prior_returns_numpy = process_prior ( prior ) # Check simulator, returns PyTorch simulator able to simulate batches. simulator = process_simulator ( simulator , prior , prior_returns_numpy ) # Consistency check after making ready for sbi. check_sbi_inputs ( simulator , prior ) return simulator , prior sbi . inference . base . simulate_for_sbi ( simulator , proposal , num_simulations , num_workers = 1 , simulation_batch_size = 1 , seed = None , show_progress_bar = True ) \u00b6 Returns ( \\(\\theta, x\\) ) pairs obtained from sampling the proposal and simulating. This function performs two steps: Sample parameters \\(\\theta\\) from the proposal . Simulate these parameters to obtain \\(x\\) . Parameters: Name Type Description Default simulator Callable A function that takes parameters \\(\\theta\\) and maps them to simulations, or observations, x , \\(\\text{sim}(\\theta)\\to x\\) . Any regular Python callable (i.e. function or class with __call__ method) can be used. required proposal Any Probability distribution that the parameters \\(\\theta\\) are sampled from. required num_simulations int Number of simulations that are run. required num_workers int Number of parallel workers to use for simulations. 1 simulation_batch_size int Number of parameter sets that the simulator maps to data x at once. If None, we simulate all parameter sets at the same time. If >= 1, the simulator has to process data of shape (simulation_batch_size, parameter_dimension). 1 seed Optional[int] Seed for reproducibility. None show_progress_bar bool Whether to show a progress bar for simulating. This will not affect whether there will be a progressbar while drawing samples from the proposal. True Returns: Sampled parameters \\(\\theta\\) and simulation-outputs \\(x\\) . Source code in sbi/inference/base.py def simulate_for_sbi ( simulator : Callable , proposal : Any , num_simulations : int , num_workers : int = 1 , simulation_batch_size : int = 1 , seed : Optional [ int ] = None , show_progress_bar : bool = True , ) -> Tuple [ Tensor , Tensor ]: r \"\"\"Returns ($\\theta, x$) pairs obtained from sampling the proposal and simulating. This function performs two steps: - Sample parameters $\\theta$ from the `proposal`. - Simulate these parameters to obtain $x$. Args: simulator: A function that takes parameters $\\theta$ and maps them to simulations, or observations, `x`, $\\text{sim}(\\theta)\\to x$. Any regular Python callable (i.e. function or class with `__call__` method) can be used. proposal: Probability distribution that the parameters $\\theta$ are sampled from. num_simulations: Number of simulations that are run. num_workers: Number of parallel workers to use for simulations. simulation_batch_size: Number of parameter sets that the simulator maps to data x at once. If None, we simulate all parameter sets at the same time. If >= 1, the simulator has to process data of shape (simulation_batch_size, parameter_dimension). seed: Seed for reproducibility. show_progress_bar: Whether to show a progress bar for simulating. This will not affect whether there will be a progressbar while drawing samples from the proposal. Returns: Sampled parameters $\\theta$ and simulation-outputs $x$. \"\"\" theta = proposal . sample (( num_simulations ,)) x = simulate_in_batches ( simulator = simulator , theta = theta , sim_batch_size = simulation_batch_size , num_workers = num_workers , seed = seed , show_progress_bars = show_progress_bar , ) return theta , x sbi.inference.snpe.snpe_a.SNPE_A ( PosteriorEstimator ) \u00b6 __init__ ( self , prior = None , density_estimator = 'mdn_snpe_a' , num_components = 10 , device = 'cpu' , logging_level = 'WARNING' , summary_writer = None , show_progress_bars = True ) special \u00b6 SNPE-A [1]. [1] Fast epsilon-free Inference of Simulation Models with Bayesian Conditional Density Estimation , Papamakarios et al., NeurIPS 2016, https://arxiv.org/abs/1605.06376 . This class implements SNPE-A. SNPE-A trains across multiple rounds with a maximum-likelihood-loss. This will make training converge to the proposal posterior instead of the true posterior. To correct for this, SNPE-A applies a post-hoc correction after training. This correction has to be performed analytically. Thus, SNPE-A is limited to Gaussian distributions for all but the last round. In the last round, SNPE-A can use a Mixture of Gaussians. Parameters: Name Type Description Default prior Optional[torch.distributions.distribution.Distribution] A probability distribution that expresses prior knowledge about the parameters, e.g. which ranges are meaningful for them. Any object with .log_prob() and .sample() (for example, a PyTorch distribution) can be used. None density_estimator Union[str, Callable] If it is a string (only \u201cmdn_snpe_a\u201d is valid), use a pre-configured mixture of densities network. Alternatively, a function that builds a custom neural network can be provided. The function will be called with the first batch of simulations (theta, x), which can thus be used for shape inference and potentially for z-scoring. It needs to return a PyTorch nn.Module implementing the density estimator. The density estimator needs to provide the methods .log_prob and .sample() . Note that until the last round only a single (multivariate) Gaussian component is used for training (see Algorithm 1 in [1]). In the last round, this component is replicated num_components times, its parameters are perturbed with a very small noise, and then the last training round is done with the expanded Gaussian mixture as estimator for the proposal posterior. 'mdn_snpe_a' num_components int Number of components of the mixture of Gaussians in the last round. This overrides the num_components value passed to posterior_nn() . 10 device str Training device, e.g., \u201ccpu\u201d, \u201ccuda\u201d or \u201ccuda:{0, 1, \u2026}\u201d. 'cpu' logging_level Union[int, str] Minimum severity of messages to log. One of the strings INFO, WARNING, DEBUG, ERROR and CRITICAL. 'WARNING' summary_writer Optional[Writer] A tensorboard SummaryWriter to control, among others, log file location (default is
- proposal is a `DirectPosterior` with density_estimator `mdn`, as built with `utils.sbi.posterior_nn()`.
- the density estimator is a `mdn`, as built with `utils.sbi.posterior_nn()`.
- `isinstance(prior, MultivariateNormal)` (from `torch.distributions`) or `isinstance(prior, sbi.utils.BoxUniform)` Note that custom implementations of any of these densities (or estimators) will not trigger the non-atomic loss, and the algorithm will fall back onto using the atomic loss. Args: prior: A probability distribution that expresses prior knowledge about the parameters, e.g. which ranges are meaningful for them. density_estimator: If it is a string, use a pre-configured network of the provided type (one of nsf, maf, mdn, made). Alternatively, a function that builds a custom neural network can be provided. The function will be called with the first batch of simulations (theta, x), which can thus be used for shape inference and potentially for z-scoring. It needs to return a PyTorch `nn.Module` implementing the density estimator. The density estimator needs to provide the methods `.log_prob` and `.sample()`. device: Training device, e.g., \"cpu\", \"cuda\" or \"cuda:{0, 1, ...}\". logging_level: Minimum severity of messages to log. One of the strings INFO, WARNING, DEBUG, ERROR and CRITICAL. summary_writer: A tensorboard `SummaryWriter` to control, among others, log file location (default is `
- the simulator can simulate batches of parameters and return batches of data.
- the prior does not produce batches and samples and evaluates to Tensor.
- the output shape is a `torch.Size((1,N))` (i.e, has a leading batch dimension 1). If this is not possible, a suitable exception will be raised. Args: simulator: Simulator as provided by the user. prior: Prior as provided by the user. Returns: Tuple (simulator, prior) checked and matching the requirements of sbi. \"\"\" # Check prior, return PyTorch prior. prior , _ , prior_returns_numpy = process_prior ( prior ) # Check simulator, returns PyTorch simulator able to simulate batches. simulator = process_simulator ( simulator , prior , prior_returns_numpy ) # Consistency check after making ready for sbi. check_sbi_inputs ( simulator , prior ) return simulator , prior","title":"prepare_for_sbi()"},{"location":"reference/#sbi.inference.base.simulate_for_sbi","text":"Returns ( \\(\\theta, x\\) ) pairs obtained from sampling the proposal and simulating. This function performs two steps: Sample parameters \\(\\theta\\) from the proposal . Simulate these parameters to obtain \\(x\\) . Parameters: Name Type Description Default simulator Callable A function that takes parameters \\(\\theta\\) and maps them to simulations, or observations, x , \\(\\text{sim}(\\theta)\\to x\\) . Any regular Python callable (i.e. function or class with __call__ method) can be used. required proposal Any Probability distribution that the parameters \\(\\theta\\) are sampled from. required num_simulations int Number of simulations that are run. required num_workers int Number of parallel workers to use for simulations. 1 simulation_batch_size int Number of parameter sets that the simulator maps to data x at once. If None, we simulate all parameter sets at the same time. If >= 1, the simulator has to process data of shape (simulation_batch_size, parameter_dimension). 1 seed Optional[int] Seed for reproducibility. None show_progress_bar bool Whether to show a progress bar for simulating. This will not affect whether there will be a progressbar while drawing samples from the proposal. True Returns: Sampled parameters \\(\\theta\\) and simulation-outputs \\(x\\) . Source code in sbi/inference/base.py def simulate_for_sbi ( simulator : Callable , proposal : Any , num_simulations : int , num_workers : int = 1 , simulation_batch_size : int = 1 , seed : Optional [ int ] = None , show_progress_bar : bool = True , ) -> Tuple [ Tensor , Tensor ]: r \"\"\"Returns ($\\theta, x$) pairs obtained from sampling the proposal and simulating. This function performs two steps: - Sample parameters $\\theta$ from the `proposal`. - Simulate these parameters to obtain $x$. Args: simulator: A function that takes parameters $\\theta$ and maps them to simulations, or observations, `x`, $\\text{sim}(\\theta)\\to x$. Any regular Python callable (i.e. function or class with `__call__` method) can be used. proposal: Probability distribution that the parameters $\\theta$ are sampled from. num_simulations: Number of simulations that are run. num_workers: Number of parallel workers to use for simulations. simulation_batch_size: Number of parameter sets that the simulator maps to data x at once. If None, we simulate all parameter sets at the same time. If >= 1, the simulator has to process data of shape (simulation_batch_size, parameter_dimension). seed: Seed for reproducibility. show_progress_bar: Whether to show a progress bar for simulating. This will not affect whether there will be a progressbar while drawing samples from the proposal. Returns: Sampled parameters $\\theta$ and simulation-outputs $x$. \"\"\" theta = proposal . sample (( num_simulations ,)) x = simulate_in_batches ( simulator = simulator , theta = theta , sim_batch_size = simulation_batch_size , num_workers = num_workers , seed = seed , show_progress_bars = show_progress_bar , ) return theta , x","title":"simulate_for_sbi()"},{"location":"reference/#sbi.inference.snpe.snpe_a.SNPE_A","text":"","title":"SNPE_A"},{"location":"reference/#sbi.inference.snpe.snpe_c.SNPE_C","text":"","title":"SNPE_C"},{"location":"reference/#sbi.inference.snle.snle_a.SNLE_A","text":"","title":"SNLE_A"},{"location":"reference/#sbi.inference.snre.snre_a.SNRE_A","text":"","title":"SNRE_A"},{"location":"reference/#sbi.inference.snre.snre_b.SNRE_B","text":"","title":"SNRE_B"},{"location":"reference/#sbi.inference.snre.snre_c.SNRE_C","text":"","title":"SNRE_C"},{"location":"reference/#sbi.inference.snre.bnre.BNRE","text":"","title":"BNRE"},{"location":"reference/#sbi.inference.abc.mcabc.MCABC","text":"","title":"MCABC"},{"location":"reference/#sbi.inference.abc.smcabc.SMCABC","text":"","title":"SMCABC"},{"location":"reference/#posteriors","text":"","title":"Posteriors"},{"location":"reference/#sbi.inference.posteriors.direct_posterior.DirectPosterior","text":"Posterior \\(p(\\theta|x_o)\\) with log_prob() and sample() methods, only applicable to SNPE. SNPE trains a neural network to directly approximate the posterior distribution. However, for bounded priors, the neural network can have leakage: it puts non-zero mass in regions where the prior is zero. The DirectPosterior class wraps the trained network to deal with these cases. Specifically, this class offers the following functionality: - correct the calculation of the log probability such that it compensates for the leakage. - reject samples that lie outside of the prior bounds. This class can not be used in combination with SNLE or SNRE.","title":"DirectPosterior"},{"location":"reference/#sbi.inference.posteriors.importance_posterior.ImportanceSamplingPosterior","text":"Provides importance sampling to sample from the posterior. SNLE or SNRE train neural networks to approximate the likelihood(-ratios). ImportanceSamplingPosterior allows to estimate the posterior log-probability by estimating the normlalization constant with importance sampling. It also allows to perform importance sampling (with .sample() ) and to draw approximate samples with sampling-importance-resampling (SIR) (with .sir_sample() )","title":"ImportanceSamplingPosterior"},{"location":"reference/#sbi.inference.posteriors.mcmc_posterior.MCMCPosterior","text":"Provides MCMC to sample from the posterior. SNLE or SNRE train neural networks to approximate the likelihood(-ratios). MCMCPosterior allows to sample from the posterior with MCMC.","title":"MCMCPosterior"},{"location":"reference/#sbi.inference.posteriors.rejection_posterior.RejectionPosterior","text":"Provides rejection sampling to sample from the posterior. SNLE or SNRE train neural networks to approximate the likelihood(-ratios). RejectionPosterior allows to sample from the posterior with rejection sampling.","title":"RejectionPosterior"},{"location":"reference/#sbi.inference.posteriors.vi_posterior.VIPosterior","text":"Provides VI (Variational Inference) to sample from the posterior. SNLE or SNRE train neural networks to approximate the likelihood(-ratios). VIPosterior allows to learn a tractable variational posterior \\(q(\\theta)\\) which approximates the true posterior \\(p(\\theta|x_o)\\) . After this second training stage, we can produce approximate posterior samples, by just sampling from q with no additional cost. For additional information see [1] and [2]. References: [1] Variational methods for simulation-based inference, Manuel Gl\u00f6ckler, Michael Deistler, Jakob Macke, 2022, https://openreview.net/forum?id=kZ0UYdhqkNY [2] Sequential Neural Posterior and Likelihood Approximation, Samuel Wiqvist, Jes Frellsen, Umberto Picchini, 2021, https://arxiv.org/abs/2102.06522","title":"VIPosterior"},{"location":"reference/#models","text":"","title":"Models"},{"location":"reference/#sbi.utils.get_nn_models.posterior_nn","text":"Returns a function that builds a density estimator for learning the posterior. This function will usually be used for SNPE. The returned function is to be passed to the inference class when using the flexible interface. Parameters: Name Type Description Default model str The type of density estimator that will be created. One of [ mdn , made , maf , maf_rqs , nsf ]. required z_score_theta Optional[str] Whether to z-score parameters \\(\\theta\\) before passing them into the network, can take one of the following: - none , or None: do not z-score. - independent : z-score each dimension independently. - structured : treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. 'independent' z_score_x Optional[str] Whether to z-score simulation outputs \\(x\\) before passing them into the network, same options as z_score_theta. 'independent' hidden_features int Number of hidden features. 50 num_transforms int Number of transforms when a flow is used. Only relevant if density estimator is a normalizing flow (i.e. currently either a maf or a nsf ). Ignored if density estimator is a mdn or made . 5 num_bins int Number of bins used for the splines in nsf . Ignored if density estimator not nsf . 10 embedding_net Module Optional embedding network for simulation outputs \\(x\\) . This embedding net allows to learn features from potentially high-dimensional simulation outputs. Identity() num_components int Number of mixture components for a mixture of Gaussians. Ignored if density estimator is not an mdn. 10 kwargs additional custom arguments passed to downstream build functions. {} Source code in sbi/utils/get_nn_models.py def posterior_nn ( model : str , z_score_theta : Optional [ str ] = \"independent\" , z_score_x : Optional [ str ] = \"independent\" , hidden_features : int = 50 , num_transforms : int = 5 , num_bins : int = 10 , embedding_net : nn . Module = nn . Identity (), num_components : int = 10 , ** kwargs , ) -> Callable : r \"\"\" Returns a function that builds a density estimator for learning the posterior. This function will usually be used for SNPE. The returned function is to be passed to the inference class when using the flexible interface. Args: model: The type of density estimator that will be created. One of [`mdn`, `made`, `maf`, `maf_rqs`, `nsf`]. z_score_theta: Whether to z-score parameters $\\theta$ before passing them into the network, can take one of the following: - `none`, or None: do not z-score. - `independent`: z-score each dimension independently. - `structured`: treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. z_score_x: Whether to z-score simulation outputs $x$ before passing them into the network, same options as z_score_theta. hidden_features: Number of hidden features. num_transforms: Number of transforms when a flow is used. Only relevant if density estimator is a normalizing flow (i.e. currently either a `maf` or a `nsf`). Ignored if density estimator is a `mdn` or `made`. num_bins: Number of bins used for the splines in `nsf`. Ignored if density estimator not `nsf`. embedding_net: Optional embedding network for simulation outputs $x$. This embedding net allows to learn features from potentially high-dimensional simulation outputs. num_components: Number of mixture components for a mixture of Gaussians. Ignored if density estimator is not an mdn. kwargs: additional custom arguments passed to downstream build functions. \"\"\" kwargs = dict ( zip ( ( \"z_score_x\" , \"z_score_y\" , \"hidden_features\" , \"num_transforms\" , \"num_bins\" , \"embedding_net\" , \"num_components\" , ), ( z_score_theta , z_score_x , hidden_features , num_transforms , num_bins , embedding_net , num_components , ), ), ** kwargs , ) def build_fn_snpe_a ( batch_theta , batch_x , num_components ): \"\"\"Build function for SNPE-A Extract the number of components from the kwargs, such that they are exposed as a kwargs, offering the possibility to later override this kwarg with `functools.partial`. This is necessary in order to make sure that the MDN in SNPE-A only has one component when running the Algorithm 1 part. \"\"\" return build_mdn ( batch_x = batch_theta , batch_y = batch_x , num_components = num_components , ** kwargs , ) def build_fn ( batch_theta , batch_x ): if model == \"mdn\" : return build_mdn ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) elif model == \"made\" : return build_made ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) elif model == \"maf\" : return build_maf ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) elif model == \"maf_rqs\" : return build_maf_rqs ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) elif model == \"nsf\" : return build_nsf ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) else : raise NotImplementedError if model == \"mdn_snpe_a\" : if num_components != 10 : raise ValueError ( \"You set `num_components`. For SNPE-A, this has to be done at \" \"instantiation of the inference object, i.e. \" \"`inference = SNPE_A(..., num_components=20)`\" ) kwargs . pop ( \"num_components\" ) return build_fn_snpe_a if model == \"mdn_snpe_a\" else build_fn","title":"posterior_nn()"},{"location":"reference/#sbi.utils.get_nn_models.likelihood_nn","text":"Returns a function that builds a density estimator for learning the likelihood. This function will usually be used for SNLE. The returned function is to be passed to the inference class when using the flexible interface. Parameters: Name Type Description Default model str The type of density estimator that will be created. One of [ mdn , made , maf , maf_rqs , nsf ]. required z_score_theta Optional[str] Whether to z-score parameters \\(\\theta\\) before passing them into the network, can take one of the following: - none , or None: do not z-score. - independent : z-score each dimension independently. - structured : treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. 'independent' z_score_x Optional[str] Whether to z-score simulation outputs \\(x\\) before passing them into the network, same options as z_score_theta. 'independent' hidden_features int Number of hidden features. 50 num_transforms int Number of transforms when a flow is used. Only relevant if density estimator is a normalizing flow (i.e. currently either a maf or a nsf ). Ignored if density estimator is a mdn or made . 5 num_bins int Number of bins used for the splines in nsf . Ignored if density estimator not nsf . 10 embedding_net Module Optional embedding network for parameters \\(\\theta\\) . Identity() num_components int Number of mixture components for a mixture of Gaussians. Ignored if density estimator is not an mdn. 10 kwargs additional custom arguments passed to downstream build functions. {} Source code in sbi/utils/get_nn_models.py def likelihood_nn ( model : str , z_score_theta : Optional [ str ] = \"independent\" , z_score_x : Optional [ str ] = \"independent\" , hidden_features : int = 50 , num_transforms : int = 5 , num_bins : int = 10 , embedding_net : nn . Module = nn . Identity (), num_components : int = 10 , ** kwargs , ) -> Callable : r \"\"\" Returns a function that builds a density estimator for learning the likelihood. This function will usually be used for SNLE. The returned function is to be passed to the inference class when using the flexible interface. Args: model: The type of density estimator that will be created. One of [`mdn`, `made`, `maf`, `maf_rqs`, `nsf`]. z_score_theta: Whether to z-score parameters $\\theta$ before passing them into the network, can take one of the following: - `none`, or None: do not z-score. - `independent`: z-score each dimension independently. - `structured`: treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. z_score_x: Whether to z-score simulation outputs $x$ before passing them into the network, same options as z_score_theta. hidden_features: Number of hidden features. num_transforms: Number of transforms when a flow is used. Only relevant if density estimator is a normalizing flow (i.e. currently either a `maf` or a `nsf`). Ignored if density estimator is a `mdn` or `made`. num_bins: Number of bins used for the splines in `nsf`. Ignored if density estimator not `nsf`. embedding_net: Optional embedding network for parameters $\\theta$. num_components: Number of mixture components for a mixture of Gaussians. Ignored if density estimator is not an mdn. kwargs: additional custom arguments passed to downstream build functions. \"\"\" kwargs = dict ( zip ( ( \"z_score_x\" , \"z_score_y\" , \"hidden_features\" , \"num_transforms\" , \"num_bins\" , \"embedding_net\" , \"num_components\" , ), ( z_score_x , z_score_theta , hidden_features , num_transforms , num_bins , embedding_net , num_components , ), ), ** kwargs , ) def build_fn ( batch_theta , batch_x ): if model == \"mdn\" : return build_mdn ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"made\" : return build_made ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"maf\" : return build_maf ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"maf_rqs\" : return build_maf_rqs ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"nsf\" : return build_nsf ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) elif model == \"mnle\" : return build_mnle ( batch_x = batch_x , batch_y = batch_theta , ** kwargs ) else : raise NotImplementedError return build_fn","title":"likelihood_nn()"},{"location":"reference/#sbi.utils.get_nn_models.classifier_nn","text":"Returns a function that builds a classifier for learning density ratios. This function will usually be used for SNRE. The returned function is to be passed to the inference class when using the flexible interface. Note that in the view of the SNRE classifier we build below, x=theta and y=x. Parameters: Name Type Description Default model str The type of classifier that will be created. One of [ linear , mlp , resnet ]. required z_score_theta Optional[str] Whether to z-score parameters \\(\\theta\\) before passing them into the network, can take one of the following: - none , or None: do not z-score. - independent : z-score each dimension independently. - structured : treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. 'independent' z_score_x Optional[str] Whether to z-score simulation outputs \\(x\\) before passing them into the network, same options as z_score_theta. 'independent' hidden_features int Number of hidden features. 50 embedding_net_theta Module Optional embedding network for parameters \\(\\theta\\) . Identity() embedding_net_x Module Optional embedding network for simulation outputs \\(x\\) . This embedding net allows to learn features from potentially high-dimensional simulation outputs. Identity() kwargs additional custom arguments passed to downstream build functions. {} Source code in sbi/utils/get_nn_models.py def classifier_nn ( model : str , z_score_theta : Optional [ str ] = \"independent\" , z_score_x : Optional [ str ] = \"independent\" , hidden_features : int = 50 , embedding_net_theta : nn . Module = nn . Identity (), embedding_net_x : nn . Module = nn . Identity (), ** kwargs , ) -> Callable : r \"\"\" Returns a function that builds a classifier for learning density ratios. This function will usually be used for SNRE. The returned function is to be passed to the inference class when using the flexible interface. Note that in the view of the SNRE classifier we build below, x=theta and y=x. Args: model: The type of classifier that will be created. One of [`linear`, `mlp`, `resnet`]. z_score_theta: Whether to z-score parameters $\\theta$ before passing them into the network, can take one of the following: - `none`, or None: do not z-score. - `independent`: z-score each dimension independently. - `structured`: treat dimensions as related, therefore compute mean and std over the entire batch, instead of per-dimension. Should be used when each sample is, for example, a time series or an image. z_score_x: Whether to z-score simulation outputs $x$ before passing them into the network, same options as z_score_theta. hidden_features: Number of hidden features. embedding_net_theta: Optional embedding network for parameters $\\theta$. embedding_net_x: Optional embedding network for simulation outputs $x$. This embedding net allows to learn features from potentially high-dimensional simulation outputs. kwargs: additional custom arguments passed to downstream build functions. \"\"\" kwargs = dict ( zip ( ( \"z_score_x\" , \"z_score_y\" , \"hidden_features\" , \"embedding_net_x\" , \"embedding_net_y\" , ), ( z_score_theta , z_score_x , hidden_features , embedding_net_theta , embedding_net_x , ), ), ** kwargs , ) def build_fn ( batch_theta , batch_x ): if model == \"linear\" : return build_linear_classifier ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) if model == \"mlp\" : return build_mlp_classifier ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) if model == \"resnet\" : return build_resnet_classifier ( batch_x = batch_theta , batch_y = batch_x , ** kwargs ) else : raise NotImplementedError return build_fn","title":"classifier_nn()"},{"location":"reference/#potentials","text":"","title":"Potentials"},{"location":"reference/#sbi.inference.potentials.posterior_based_potential.posterior_estimator_based_potential","text":"Returns the potential for posterior-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. The potential is the same as the log-probability of the posterior_estimator , but it is set to \\(-\\inf\\) outside of the prior bounds. Parameters: Name Type Description Default posterior_estimator Module The neural network modelling the posterior. required prior Distribution The prior distribution. required x_o Optional[torch.Tensor] The observed data at which to evaluate the posterior. required enable_transform bool Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for theta_transform . True Returns: Type Description Tuple[Callable, torch Transform] The potential function and a transformation that maps to unconstrained space. Source code in sbi/inference/potentials/posterior_based_potential.py def posterior_estimator_based_potential ( posterior_estimator : nn . Module , prior : Distribution , x_o : Optional [ Tensor ], enable_transform : bool = True , ) -> Tuple [ Callable , TorchTransform ]: r \"\"\"Returns the potential for posterior-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. The potential is the same as the log-probability of the `posterior_estimator`, but it is set to $-\\inf$ outside of the prior bounds. Args: posterior_estimator: The neural network modelling the posterior. prior: The prior distribution. x_o: The observed data at which to evaluate the posterior. enable_transform: Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for `theta_transform`. Returns: The potential function and a transformation that maps to unconstrained space. \"\"\" device = str ( next ( posterior_estimator . parameters ()) . device ) potential_fn = PosteriorBasedPotential ( posterior_estimator , prior , x_o , device = device ) theta_transform = mcmc_transform ( prior , device = device , enable_transform = enable_transform ) return potential_fn , theta_transform","title":"posterior_estimator_based_potential()"},{"location":"reference/#sbi.inference.potentials.likelihood_based_potential.likelihood_estimator_based_potential","text":"Returns potential \\(\\log(p(x_o|\\theta)p(\\theta))\\) for likelihood-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. Parameters: Name Type Description Default likelihood_estimator Module The neural network modelling the likelihood. required prior Distribution The prior distribution. required x_o Optional[torch.Tensor] The observed data at which to evaluate the likelihood. required enable_transform bool Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for theta_transform . True Returns: Type Description Tuple[Callable, torch Transform] The potential function \\(p(x_o|\\theta)p(\\theta)\\) and a transformation that maps to unconstrained space. Source code in sbi/inference/potentials/likelihood_based_potential.py def likelihood_estimator_based_potential ( likelihood_estimator : nn . Module , prior : Distribution , x_o : Optional [ Tensor ], enable_transform : bool = True , ) -> Tuple [ Callable , TorchTransform ]: r \"\"\"Returns potential $\\log(p(x_o|\\theta)p(\\theta))$ for likelihood-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. Args: likelihood_estimator: The neural network modelling the likelihood. prior: The prior distribution. x_o: The observed data at which to evaluate the likelihood. enable_transform: Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for `theta_transform`. Returns: The potential function $p(x_o|\\theta)p(\\theta)$ and a transformation that maps to unconstrained space. \"\"\" device = str ( next ( likelihood_estimator . parameters ()) . device ) potential_fn = LikelihoodBasedPotential ( likelihood_estimator , prior , x_o , device = device ) theta_transform = mcmc_transform ( prior , device = device , enable_transform = enable_transform ) return potential_fn , theta_transform","title":"likelihood_estimator_based_potential()"},{"location":"reference/#sbi.inference.potentials.ratio_based_potential.ratio_estimator_based_potential","text":"Returns the potential for ratio-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. Parameters: Name Type Description Default ratio_estimator Module The neural network modelling likelihood-to-evidence ratio. required prior Distribution The prior distribution. required x_o Optional[torch.Tensor] The observed data at which to evaluate the likelihood-to-evidence ratio. required enable_transform bool Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for theta_transform . True Returns: Type Description Tuple[Callable, torch Transform] The potential function and a transformation that maps to unconstrained space. Source code in sbi/inference/potentials/ratio_based_potential.py def ratio_estimator_based_potential ( ratio_estimator : nn . Module , prior : Distribution , x_o : Optional [ Tensor ], enable_transform : bool = True , ) -> Tuple [ Callable , TorchTransform ]: r \"\"\"Returns the potential for ratio-based methods. It also returns a transformation that can be used to transform the potential into unconstrained space. Args: ratio_estimator: The neural network modelling likelihood-to-evidence ratio. prior: The prior distribution. x_o: The observed data at which to evaluate the likelihood-to-evidence ratio. enable_transform: Whether to transform parameters to unconstrained space. When False, an identity transform will be returned for `theta_transform`. Returns: The potential function and a transformation that maps to unconstrained space. \"\"\" device = str ( next ( ratio_estimator . parameters ()) . device ) potential_fn = RatioBasedPotential ( ratio_estimator , prior , x_o , device = device ) theta_transform = mcmc_transform ( prior , device = device , enable_transform = enable_transform ) return potential_fn , theta_transform","title":"ratio_estimator_based_potential()"},{"location":"reference/#analysis","text":"","title":"Analysis"},{"location":"reference/#sbi.analysis.plot.pairplot","text":"Plot samples in a 2D grid showing marginals and pairwise marginals. Each of the diagonal plots can be interpreted as a 1D-marginal of the distribution that the samples were drawn from. Each upper-diagonal plot can be interpreted as a 2D-marginal of the distribution. Parameters: Name Type Description Default samples Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] Samples used to build the histogram. required points Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] List of additional points to scatter. None limits Union[List, torch.Tensor] Array containing the plot xlim for each parameter dimension. If None, just use the min and max of the passed samples None subset Optional[List[int]] List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1 st and 3 rd dimension but will discard the 0 th and 2 nd (and, if they exist, the 4 th , 5 th and so on). None offdiag Union[str, List[str]] Plotting style for upper diagonal, {hist, scatter, contour, cond, None}. 'hist' upper Optional[str] deprecated, use offdiag instead. None diag Union[str, List[str]] Plotting style for diagonal, {hist, cond, None}. 'hist' figsize Tuple Size of the entire figure. (10, 10) labels Optional[List[str]] List of strings specifying the names of the parameters. None ticks Union[List, torch.Tensor] Position of the ticks. [] fig matplotlib figure to plot on. None axes matplotlib axes corresponding to fig. None **kwargs Additional arguments to adjust the plot, e.g., samples_colors , points_colors and many more, see the source code in _get_default_opts() in sbi.analysis.plot for details. {} Returns: figure and axis of posterior distribution plot Source code in sbi/analysis/plot.py def pairplot ( samples : Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ], points : Optional [ Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ] ] = None , limits : Optional [ Union [ List , torch . Tensor ]] = None , subset : Optional [ List [ int ]] = None , offdiag : Optional [ Union [ List [ str ], str ]] = \"hist\" , diag : Optional [ Union [ List [ str ], str ]] = \"hist\" , figsize : Tuple = ( 10 , 10 ), labels : Optional [ List [ str ]] = None , ticks : Union [ List , torch . Tensor ] = [], upper : Optional [ str ] = None , fig = None , axes = None , ** kwargs , ): \"\"\" Plot samples in a 2D grid showing marginals and pairwise marginals. Each of the diagonal plots can be interpreted as a 1D-marginal of the distribution that the samples were drawn from. Each upper-diagonal plot can be interpreted as a 2D-marginal of the distribution. Args: samples: Samples used to build the histogram. points: List of additional points to scatter. limits: Array containing the plot xlim for each parameter dimension. If None, just use the min and max of the passed samples subset: List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1st and 3rd dimension but will discard the 0th and 2nd (and, if they exist, the 4th, 5th and so on). offdiag: Plotting style for upper diagonal, {hist, scatter, contour, cond, None}. upper: deprecated, use offdiag instead. diag: Plotting style for diagonal, {hist, cond, None}. figsize: Size of the entire figure. labels: List of strings specifying the names of the parameters. ticks: Position of the ticks. fig: matplotlib figure to plot on. axes: matplotlib axes corresponding to fig. **kwargs: Additional arguments to adjust the plot, e.g., `samples_colors`, `points_colors` and many more, see the source code in `_get_default_opts()` in `sbi.analysis.plot` for details. Returns: figure and axis of posterior distribution plot \"\"\" # TODO: add color map support # TODO: automatically determine good bin sizes for histograms # TODO: add legend (if legend is True) opts = _get_default_opts () # update the defaults dictionary by the current values of the variables (passed by # the user) opts = _update ( opts , locals ()) opts = _update ( opts , kwargs ) samples , dim , limits = prepare_for_plot ( samples , limits ) # checks. if opts [ \"legend\" ]: assert len ( opts [ \"samples_labels\" ]) >= len ( samples ), \"Provide at least as many labels as samples.\" if opts [ \"upper\" ] is not None : warn ( \"upper is deprecated, use offdiag instead.\" ) opts [ \"offdiag\" ] = opts [ \"upper\" ] # Prepare diag/upper/lower if type ( opts [ \"diag\" ]) is not list : opts [ \"diag\" ] = [ opts [ \"diag\" ] for _ in range ( len ( samples ))] if type ( opts [ \"offdiag\" ]) is not list : opts [ \"offdiag\" ] = [ opts [ \"offdiag\" ] for _ in range ( len ( samples ))] # if type(opts['lower']) is not list: # opts['lower'] = [opts['lower'] for _ in range(len(samples))] opts [ \"lower\" ] = None diag_func = get_diag_func ( samples , limits , opts , ** kwargs ) def offdiag_func ( row , col , limits , ** kwargs ): if len ( samples ) > 0 : for n , v in enumerate ( samples ): if opts [ \"offdiag\" ][ n ] == \"hist\" or opts [ \"offdiag\" ][ n ] == \"hist2d\" : hist , xedges , yedges = np . histogram2d ( v [:, col ], v [:, row ], range = [ [ limits [ col ][ 0 ], limits [ col ][ 1 ]], [ limits [ row ][ 0 ], limits [ row ][ 1 ]], ], ** opts [ \"hist_offdiag\" ], ) plt . imshow ( hist . T , origin = \"lower\" , extent = ( xedges [ 0 ], xedges [ - 1 ], yedges [ 0 ], yedges [ - 1 ], ), aspect = \"auto\" , ) elif opts [ \"offdiag\" ][ n ] in [ \"kde\" , \"kde2d\" , \"contour\" , \"contourf\" , ]: density = gaussian_kde ( v [:, [ col , row ]] . T , bw_method = opts [ \"kde_offdiag\" ][ \"bw_method\" ], ) X , Y = np . meshgrid ( np . linspace ( limits [ col ][ 0 ], limits [ col ][ 1 ], opts [ \"kde_offdiag\" ][ \"bins\" ], ), np . linspace ( limits [ row ][ 0 ], limits [ row ][ 1 ], opts [ \"kde_offdiag\" ][ \"bins\" ], ), ) positions = np . vstack ([ X . ravel (), Y . ravel ()]) Z = np . reshape ( density ( positions ) . T , X . shape ) if opts [ \"offdiag\" ][ n ] == \"kde\" or opts [ \"offdiag\" ][ n ] == \"kde2d\" : plt . imshow ( Z , extent = ( limits [ col ][ 0 ], limits [ col ][ 1 ], limits [ row ][ 0 ], limits [ row ][ 1 ], ), origin = \"lower\" , aspect = \"auto\" , ) elif opts [ \"offdiag\" ][ n ] == \"contour\" : if opts [ \"contour_offdiag\" ][ \"percentile\" ]: Z = probs2contours ( Z , opts [ \"contour_offdiag\" ][ \"levels\" ]) else : Z = ( Z - Z . min ()) / ( Z . max () - Z . min ()) plt . contour ( X , Y , Z , origin = \"lower\" , extent = [ limits [ col ][ 0 ], limits [ col ][ 1 ], limits [ row ][ 0 ], limits [ row ][ 1 ], ], colors = opts [ \"samples_colors\" ][ n ], levels = opts [ \"contour_offdiag\" ][ \"levels\" ], ) else : pass elif opts [ \"offdiag\" ][ n ] == \"scatter\" : plt . scatter ( v [:, col ], v [:, row ], color = opts [ \"samples_colors\" ][ n ], ** opts [ \"scatter_offdiag\" ], ) elif opts [ \"offdiag\" ][ n ] == \"plot\" : plt . plot ( v [:, col ], v [:, row ], color = opts [ \"samples_colors\" ][ n ], ** opts [ \"plot_offdiag\" ], ) else : pass return _arrange_plots ( diag_func , offdiag_func , dim , limits , points , opts , fig = fig , axes = axes )","title":"pairplot()"},{"location":"reference/#sbi.analysis.plot.marginal_plot","text":"Plot samples in a row showing 1D marginals of selected dimensions. Each of the plots can be interpreted as a 1D-marginal of the distribution that the samples were drawn from. Parameters: Name Type Description Default samples Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] Samples used to build the histogram. required points Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] List of additional points to scatter. None limits Union[List, torch.Tensor] Array containing the plot xlim for each parameter dimension. If None, just use the min and max of the passed samples None subset Optional[List[int]] List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1 st and 3 rd dimension but will discard the 0 th and 2 nd (and, if they exist, the 4 th , 5 th and so on). None diag Optional[str] Plotting style for 1D marginals, {hist, kde cond, None}. 'hist' figsize Tuple Size of the entire figure. (10, 10) labels Optional[List[str]] List of strings specifying the names of the parameters. None ticks Union[List, torch.Tensor] Position of the ticks. [] points_colors Colors of the points . required fig matplotlib figure to plot on. None axes matplotlib axes corresponding to fig. None **kwargs Additional arguments to adjust the plot, e.g., samples_colors , points_colors and many more, see the source code in _get_default_opts() in sbi.analysis.plot for details. {} Returns: figure and axis of posterior distribution plot Source code in sbi/analysis/plot.py def marginal_plot ( samples : Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ], points : Optional [ Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ] ] = None , limits : Optional [ Union [ List , torch . Tensor ]] = None , subset : Optional [ List [ int ]] = None , diag : Optional [ str ] = \"hist\" , figsize : Tuple = ( 10 , 10 ), labels : Optional [ List [ str ]] = None , ticks : Union [ List , torch . Tensor ] = [], fig = None , axes = None , ** kwargs , ): \"\"\" Plot samples in a row showing 1D marginals of selected dimensions. Each of the plots can be interpreted as a 1D-marginal of the distribution that the samples were drawn from. Args: samples: Samples used to build the histogram. points: List of additional points to scatter. limits: Array containing the plot xlim for each parameter dimension. If None, just use the min and max of the passed samples subset: List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1st and 3rd dimension but will discard the 0th and 2nd (and, if they exist, the 4th, 5th and so on). diag: Plotting style for 1D marginals, {hist, kde cond, None}. figsize: Size of the entire figure. labels: List of strings specifying the names of the parameters. ticks: Position of the ticks. points_colors: Colors of the `points`. fig: matplotlib figure to plot on. axes: matplotlib axes corresponding to fig. **kwargs: Additional arguments to adjust the plot, e.g., `samples_colors`, `points_colors` and many more, see the source code in `_get_default_opts()` in `sbi.analysis.plot` for details. Returns: figure and axis of posterior distribution plot \"\"\" opts = _get_default_opts () # update the defaults dictionary by the current values of the variables (passed by # the user) opts = _update ( opts , locals ()) opts = _update ( opts , kwargs ) samples , dim , limits = prepare_for_plot ( samples , limits ) # Prepare diag/upper/lower if type ( opts [ \"diag\" ]) is not list : opts [ \"diag\" ] = [ opts [ \"diag\" ] for _ in range ( len ( samples ))] diag_func = get_diag_func ( samples , limits , opts , ** kwargs ) return _arrange_plots ( diag_func , None , dim , limits , points , opts , fig = fig , axes = axes )","title":"marginal_plot()"},{"location":"reference/#sbi.analysis.plot.conditional_pairplot","text":"Plot conditional distribution given all other parameters. The conditionals can be interpreted as slices through the density at a location given by condition . For example: Say we have a 3D density with parameters \\(\\theta_0\\) , \\(\\theta_1\\) , \\(\\theta_2\\) and a condition \\(c\\) passed by the user in the condition argument. For the plot of \\(\\theta_0\\) on the diagonal, this will plot the conditional \\(p(\\theta_0 | \\theta_1=c[1], \\theta_2=c[2])\\) . For the upper diagonal of \\(\\theta_1\\) and \\(\\theta_2\\) , it will plot \\(p(\\theta_1, \\theta_2 | \\theta_0=c[0])\\) . All other diagonals and upper-diagonals are built in the corresponding way. Parameters: Name Type Description Default density Any Probability density with a log_prob() method. required condition Tensor Condition that all but the one/two regarded parameters are fixed to. The condition should be of shape (1, dim_theta), i.e. it could e.g. be a sample from the posterior distribution. required limits Union[List, torch.Tensor] Limits in between which each parameter will be evaluated. required points Union[List[numpy.ndarray], List[torch.Tensor], numpy.ndarray, torch.Tensor] Additional points to scatter. None subset Optional[List[int]] List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1 st and 3 rd dimension but will discard the 0 th and 2 nd (and, if they exist, the 4 th , 5 th and so on) None resolution int Resolution of the grid at which we evaluate the pdf . 50 figsize Tuple Size of the entire figure. (10, 10) labels Optional[List[str]] List of strings specifying the names of the parameters. None ticks Union[List, torch.Tensor] Position of the ticks. [] points_colors Colors of the points . required fig matplotlib figure to plot on. None axes matplotlib axes corresponding to fig. None **kwargs Additional arguments to adjust the plot, e.g., samples_colors , points_colors and many more, see the source code in _get_default_opts() in sbi.analysis.plot for details. {} Returns: figure and axis of posterior distribution plot Source code in sbi/analysis/plot.py def conditional_pairplot ( density : Any , condition : torch . Tensor , limits : Union [ List , torch . Tensor ], points : Optional [ Union [ List [ np . ndarray ], List [ torch . Tensor ], np . ndarray , torch . Tensor ] ] = None , subset : Optional [ List [ int ]] = None , resolution : int = 50 , figsize : Tuple = ( 10 , 10 ), labels : Optional [ List [ str ]] = None , ticks : Union [ List , torch . Tensor ] = [], fig = None , axes = None , ** kwargs , ): r \"\"\" Plot conditional distribution given all other parameters. The conditionals can be interpreted as slices through the `density` at a location given by `condition`. For example: Say we have a 3D density with parameters $\\theta_0$, $\\theta_1$, $\\theta_2$ and a condition $c$ passed by the user in the `condition` argument. For the plot of $\\theta_0$ on the diagonal, this will plot the conditional $p(\\theta_0 | \\theta_1=c[1], \\theta_2=c[2])$. For the upper diagonal of $\\theta_1$ and $\\theta_2$, it will plot $p(\\theta_1, \\theta_2 | \\theta_0=c[0])$. All other diagonals and upper-diagonals are built in the corresponding way. Args: density: Probability density with a `log_prob()` method. condition: Condition that all but the one/two regarded parameters are fixed to. The condition should be of shape (1, dim_theta), i.e. it could e.g. be a sample from the posterior distribution. limits: Limits in between which each parameter will be evaluated. points: Additional points to scatter. subset: List containing the dimensions to plot. E.g. subset=[1,3] will plot plot only the 1st and 3rd dimension but will discard the 0th and 2nd (and, if they exist, the 4th, 5th and so on) resolution: Resolution of the grid at which we evaluate the `pdf`. figsize: Size of the entire figure. labels: List of strings specifying the names of the parameters. ticks: Position of the ticks. points_colors: Colors of the `points`. fig: matplotlib figure to plot on. axes: matplotlib axes corresponding to fig. **kwargs: Additional arguments to adjust the plot, e.g., `samples_colors`, `points_colors` and many more, see the source code in `_get_default_opts()` in `sbi.analysis.plot` for details. Returns: figure and axis of posterior distribution plot \"\"\" device = density . _device if hasattr ( density , \"_device\" ) else \"cpu\" # Setting these is required because _pairplot_scaffold will check if opts['diag'] is # `None`. This would break if opts has no key 'diag'. Same for 'upper'. diag = \"cond\" offdiag = \"cond\" opts = _get_default_opts () # update the defaults dictionary by the current values of the variables (passed by # the user) opts = _update ( opts , locals ()) opts = _update ( opts , kwargs ) opts [ \"lower\" ] = None dim , limits , eps_margins = prepare_for_conditional_plot ( condition , opts ) diag_func = get_conditional_diag_func ( opts , limits , eps_margins , resolution ) def offdiag_func ( row , col , ** kwargs ): p_image = ( eval_conditional_density ( opts [ \"density\" ], opts [ \"condition\" ] . to ( device ), limits . to ( device ), row , col , resolution = resolution , eps_margins1 = eps_margins [ row ], eps_margins2 = eps_margins [ col ], ) . to ( \"cpu\" ) . numpy () ) plt . imshow ( p_image . T , origin = \"lower\" , extent = ( limits [ col , 0 ] . item (), limits [ col , 1 ] . item (), limits [ row , 0 ] . item (), limits [ row , 1 ] . item (), ), aspect = \"auto\" , ) return _arrange_plots ( diag_func , offdiag_func , dim , limits , points , opts , fig = fig , axes = axes )","title":"conditional_pairplot()"},{"location":"reference/#sbi.analysis.conditional_density.conditional_corrcoeff","text":"Returns the conditional correlation matrix of a distribution. To compute the conditional distribution, we condition all but two parameters to values from condition , and then compute the Pearson correlation coefficient \\(\\rho\\) between the remaining two parameters under the distribution density . We do so for any pair of parameters specified in subset , thus creating a matrix containing conditional correlations between any pair of parameters. If condition is a batch of conditions, this function computes the conditional correlation matrix for each one of them and returns the mean. Parameters: Name Type Description Default density Any Probability density function with .log_prob() function. required limits Tensor Limits within which to evaluate the density . required condition Tensor Values to condition the density on. If a batch of conditions is passed, we compute the conditional correlation matrix for each of them and return the average conditional correlation matrix. required subset Optional[List[int]] Evaluate the conditional distribution only on a subset of dimensions. If None this function uses all dimensions. None resolution int Number of grid points on which the conditional distribution is evaluated. A higher value increases the accuracy of the estimated correlation but also increases the computational cost. 50 Returns: Average conditional correlation matrix of shape either (num_dim, num_dim) or (len(subset), len(subset)) if subset was specified. Source code in sbi/analysis/conditional_density.py def conditional_corrcoeff ( density : Any , limits : Tensor , condition : Tensor , subset : Optional [ List [ int ]] = None , resolution : int = 50 , ) -> Tensor : r \"\"\"Returns the conditional correlation matrix of a distribution. To compute the conditional distribution, we condition all but two parameters to values from `condition`, and then compute the Pearson correlation coefficient $\\rho$ between the remaining two parameters under the distribution `density`. We do so for any pair of parameters specified in `subset`, thus creating a matrix containing conditional correlations between any pair of parameters. If `condition` is a batch of conditions, this function computes the conditional correlation matrix for each one of them and returns the mean. Args: density: Probability density function with `.log_prob()` function. limits: Limits within which to evaluate the `density`. condition: Values to condition the `density` on. If a batch of conditions is passed, we compute the conditional correlation matrix for each of them and return the average conditional correlation matrix. subset: Evaluate the conditional distribution only on a subset of dimensions. If `None` this function uses all dimensions. resolution: Number of grid points on which the conditional distribution is evaluated. A higher value increases the accuracy of the estimated correlation but also increases the computational cost. Returns: Average conditional correlation matrix of shape either `(num_dim, num_dim)` or `(len(subset), len(subset))` if `subset` was specified. \"\"\" device = density . _device if hasattr ( density , \"_device\" ) else \"cpu\" subset_ = subset if subset is not None else range ( condition . shape [ 1 ]) correlation_matrices = [] for cond in condition : correlation_matrices . append ( torch . stack ( [ compute_corrcoeff ( eval_conditional_density ( density , cond . to ( device ), limits . to ( device ), dim1 = dim1 , dim2 = dim2 , resolution = resolution , ), limits [[ dim1 , dim2 ]] . to ( device ), ) for dim1 in subset_ for dim2 in subset_ if dim1 < dim2 ] ) ) average_correlations = torch . mean ( torch . stack ( correlation_matrices ), dim = 0 ) # `average_correlations` is still a vector containing the upper triangular entries. # Below, assemble them into a matrix: av_correlation_matrix = torch . zeros (( len ( subset_ ), len ( subset_ )), device = device ) triu_indices = torch . triu_indices ( row = len ( subset_ ), col = len ( subset_ ), offset = 1 , device = device ) av_correlation_matrix [ triu_indices [ 0 ], triu_indices [ 1 ]] = average_correlations # Make the matrix symmetric by copying upper diagonal to lower diagonal. av_correlation_matrix = torch . triu ( av_correlation_matrix ) + torch . tril ( av_correlation_matrix . T ) av_correlation_matrix . fill_diagonal_ ( 1.0 ) return av_correlation_matrix","title":"conditional_corrcoeff()"},{"location":"examples/00_HH_simulator/","text":"Inference on Hodgkin-Huxley model: tutorial \u00b6 In this tutorial, we use sbi to do inference on a Hodgkin-Huxley model from neuroscience (Hodgkin and Huxley, 1952). We will learn two parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ) based on a current-clamp recording, that we generate synthetically (in practice, this would be an experimental observation). Note, you find the original version of this notebook at https://github.com/sbi-dev/sbi/blob/main/examples/00_HH_simulator.ipynb in the sbi repository. First we are going to import basic packages. import numpy as np import torch # visualization import matplotlib as mpl import matplotlib.pyplot as plt # sbi from sbi import utils as utils from sbi import analysis as analysis from sbi.inference.base import infer # remove top and right axis from plots mpl . rcParams [ \"axes.spines.right\" ] = False mpl . rcParams [ \"axes.spines.top\" ] = False Different required components \u00b6 Before running inference, let us define the different required components: observed data prior over model parameters simulator 1. Observed data \u00b6 Let us assume we current-clamped a neuron and recorded the following voltage trace: In fact, this voltage trace was not measured experimentally but synthetically generated by simulating a Hodgkin-Huxley model with particular parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ). We will come back to this point later in the tutorial. 2. Simulator \u00b6 We would like to infer the posterior over the two parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ) of a Hodgkin-Huxley model, given the observed electrophysiological recording above. The model has channel kinetics as in Pospischil et al. 2008 , and is defined by the following set of differential equations (parameters of interest highlighted in orange): \\[ \\scriptsize \\begin{align} C_m\\frac{dV}{dt}&=g_1\\left(E_1-V\\right)+ \\color{orange}{\\bar{g}_{Na}}m^3h\\left(E_{Na}-V\\right)+ \\color{orange}{\\bar{g}_{K}}n^4\\left(E_K-V\\right)+ \\bar{g}_Mp\\left(E_K-V\\right)+ I_{inj}+ \\sigma\\eta\\left(t\\right)\\\\ \\frac{dq}{dt}&=\\frac{q_\\infty\\left(V\\right)-q}{\\tau_q\\left(V\\right)},\\;q\\in\\{m,h,n,p\\} \\end{align} \\] Above, \\(V\\) represents the membrane potential, \\(C_m\\) is the membrane capacitance, \\(g_{\\text{l}}\\) is the leak conductance, \\(E_{\\text{l}}\\) is the membrane reversal potential, \\(\\bar{g}_c\\) is the density of channels of type \\(c\\) ( \\(\\text{Na}^+\\) , \\(\\text{K}^+\\) , M), \\(E_c\\) is the reversal potential of \\(c\\) , ( \\(m\\) , \\(h\\) , \\(n\\) , \\(p\\) ) are the respective channel gating kinetic variables, and \\(\\sigma \\eta(t)\\) is the intrinsic neural noise. The right hand side of the voltage dynamics is composed of a leak current, a voltage-dependent \\(\\text{Na}^+\\) current, a delayed-rectifier \\(\\text{K}^+\\) current, a slow voltage-dependent \\(\\text{K}^+\\) current responsible for spike-frequency adaptation, and an injected current \\(I_{\\text{inj}}\\) . Channel gating variables \\(q\\) have dynamics fully characterized by the neuron membrane potential \\(V\\) , given the respective steady-state \\(q_{\\infty}(V)\\) and time constant \\(\\tau_{q}(V)\\) (details in Pospischil et al. 2008). The input current \\(I_{\\text{inj}}\\) is defined as from HH_helper_functions import syn_current I , t_on , t_off , dt , t , A_soma = syn_current () The Hodgkin-Huxley simulator is given by: from HH_helper_functions import HHsimulator Putting the input current and the simulator together: def run_HH_model ( params ): params = np . asarray ( params ) # input current, time step I , t_on , t_off , dt , t , A_soma = syn_current () t = np . arange ( 0 , len ( I ), 1 ) * dt # initial voltage V0 = - 70 states = HHsimulator ( V0 , params . reshape ( 1 , - 1 ), dt , t , I ) return dict ( data = states . reshape ( - 1 ), time = t , dt = dt , I = I . reshape ( - 1 )) To get an idea of the output of the Hodgkin-Huxley model, let us generate some voltage traces for different parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), given the input current \\(I_{\\text{inj}}\\) : # three sets of (g_Na, g_K) params = np . array ([[ 50.0 , 1.0 ], [ 4.0 , 1.5 ], [ 20.0 , 15.0 ]]) num_samples = len ( params [:, 0 ]) sim_samples = np . zeros (( num_samples , len ( I ))) for i in range ( num_samples ): sim_samples [ i , :] = run_HH_model ( params = params [ i , :])[ \"data\" ] # colors for traces col_min = 2 num_colors = num_samples + col_min cm1 = mpl . cm . Blues col1 = [ cm1 ( 1.0 * i / num_colors ) for i in range ( col_min , num_colors )] fig = plt . figure ( figsize = ( 7 , 5 )) gs = mpl . gridspec . GridSpec ( 2 , 1 , height_ratios = [ 4 , 1 ]) ax = plt . subplot ( gs [ 0 ]) for i in range ( num_samples ): plt . plot ( t , sim_samples [ i , :], color = col1 [ i ], lw = 2 ) plt . ylabel ( \"voltage (mV)\" ) ax . set_xticks ([]) ax . set_yticks ([ - 80 , - 20 , 40 ]) ax = plt . subplot ( gs [ 1 ]) plt . plot ( t , I * A_soma * 1e3 , \"k\" , lw = 2 ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"input (nA)\" ) ax . set_xticks ([ 0 , max ( t ) / 2 , max ( t )]) ax . set_yticks ([ 0 , 1.1 * np . max ( I * A_soma * 1e3 )]) ax . yaxis . set_major_formatter ( mpl . ticker . FormatStrFormatter ( \" %.2f \" )) plt . show () As can be seen, the voltage traces can be quite diverse for different parameter values. Often, we are not interested in matching the exact trace, but only in matching certain features thereof. In this example of the Hodgkin-Huxley model, the summary features are the number of spikes, the mean resting potential, the standard deviation of the resting potential, and the first four voltage moments: mean, standard deviation, skewness and kurtosis. Using the function calculate_summary_statistics() imported below, we obtain these statistics from the output of the Hodgkin Huxley simulator. from HH_helper_functions import calculate_summary_statistics Lastly, we define a function that performs all of the above steps at once. The function simulation_wrapper takes in conductance values, runs the Hodgkin Huxley model and then returns the summary statistics. def simulation_wrapper ( params ): \"\"\" Returns summary statistics from conductance values in `params`. Summarizes the output of the HH simulator and converts it to `torch.Tensor`. \"\"\" obs = run_HH_model ( params ) summstats = torch . as_tensor ( calculate_summary_statistics ( obs )) return summstats sbi takes any function as simulator. Thus, sbi also has the flexibility to use simulators that utilize external packages, e.g., Brian ( http://briansimulator.org/ ), nest ( https://www.nest-simulator.org/ ), or NEURON ( https://neuron.yale.edu/neuron/ ). External simulators do not even need to be Python-based as long as they store simulation outputs in a format that can be read from Python. All that is necessary is to wrap your external simulator of choice into a Python callable that takes a parameter set and outputs a set of summary statistics we want to fit the parameters to. 3. Prior over model parameters \u00b6 Now that we have the simulator, we need to define a function with the prior over the model parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), which in this case is chosen to be a Uniform distribution: prior_min = [ 0.5 , 1e-4 ] prior_max = [ 80.0 , 15.0 ] prior = utils . torchutils . BoxUniform ( low = torch . as_tensor ( prior_min ), high = torch . as_tensor ( prior_max ) ) Inference \u00b6 Now that we have all the required components, we can run inference with SNPE to identify parameters whose activity matches this trace. posterior = infer ( simulation_wrapper , prior , method = \"SNPE\" , num_simulations = 300 , num_workers = 4 ) HBox(children=(FloatProgress(value=0.0, description='Running 300 simulations in 300 batches.', max=300.0, styl\u2026 Neural network successfully converged after 233 epochs. Note sbi can parallelize your simulator. If you experience problems with parallelization, try setting num_workers=1 and please give us an error report as a GitHub issue . Coming back to the observed data \u00b6 As mentioned at the beginning of the tutorial, the observed data are generated by the Hodgkin-Huxley model with a set of known parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ). To illustrate how to compute the summary statistics of the observed data, let us regenerate the observed data: # true parameters and respective labels true_params = np . array ([ 50.0 , 5.0 ]) labels_params = [ r \"$g_ {Na} $\" , r \"$g_ {K} $\" ] observation_trace = run_HH_model ( true_params ) observation_summary_statistics = calculate_summary_statistics ( observation_trace ) As we already shown above, the observed voltage traces look as follows: fig = plt . figure ( figsize = ( 7 , 5 )) gs = mpl . gridspec . GridSpec ( 2 , 1 , height_ratios = [ 4 , 1 ]) ax = plt . subplot ( gs [ 0 ]) plt . plot ( observation_trace [ \"time\" ], observation_trace [ \"data\" ]) plt . ylabel ( \"voltage (mV)\" ) plt . title ( \"observed data\" ) plt . setp ( ax , xticks = [], yticks = [ - 80 , - 20 , 40 ]) ax = plt . subplot ( gs [ 1 ]) plt . plot ( observation_trace [ \"time\" ], I * A_soma * 1e3 , \"k\" , lw = 2 ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"input (nA)\" ) ax . set_xticks ([ 0 , max ( observation_trace [ \"time\" ]) / 2 , max ( observation_trace [ \"time\" ])]) ax . set_yticks ([ 0 , 1.1 * np . max ( I * A_soma * 1e3 )]) ax . yaxis . set_major_formatter ( mpl . ticker . FormatStrFormatter ( \" %.2f \" )) Analysis of the posterior given the observed data \u00b6 After running the inference algorithm, let us inspect the inferred posterior distribution over the parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), given the observed trace. To do so, we first draw samples (i.e. consistent parameter sets) from the posterior: samples = posterior . sample (( 10000 ,), x = observation_summary_statistics ) HBox(children=(FloatProgress(value=0.0, description='Drawing 10000 posterior samples', max=10000.0, style=Prog\u2026 fig , axes = analysis . pairplot ( samples , limits = [[ 0.5 , 80 ], [ 1e-4 , 15.0 ]], ticks = [[ 0.5 , 80 ], [ 1e-4 , 15.0 ]], figsize = ( 5 , 5 ), points = true_params , points_offdiag = { \"markersize\" : 6 }, points_colors = \"r\" , ); As can be seen, the inferred posterior contains the ground-truth parameters (red) in a high-probability region. Now, let us sample parameters from the posterior distribution, simulate the Hodgkin-Huxley model for this parameter set and compare the simulations with the observed data: # Draw a sample from the posterior and convert to numpy for plotting. posterior_sample = posterior . sample (( 1 ,), x = observation_summary_statistics ) . numpy () HBox(children=(FloatProgress(value=0.0, description='Drawing 1 posterior samples', max=1.0, style=ProgressStyl\u2026 fig = plt . figure ( figsize = ( 7 , 5 )) # plot observation t = observation_trace [ \"time\" ] y_obs = observation_trace [ \"data\" ] plt . plot ( t , y_obs , lw = 2 , label = \"observation\" ) # simulate and plot samples x = run_HH_model ( posterior_sample ) plt . plot ( t , x [ \"data\" ], \"--\" , lw = 2 , label = \"posterior sample\" ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"voltage (mV)\" ) ax = plt . gca () handles , labels = ax . get_legend_handles_labels () ax . legend ( handles [:: - 1 ], labels [:: - 1 ], bbox_to_anchor = ( 1.3 , 1 ), loc = \"upper right\" ) ax . set_xticks ([ 0 , 60 , 120 ]) ax . set_yticks ([ - 80 , - 20 , 40 ]); As can be seen, the sample from the inferred posterior leads to simulations that closely resemble the observed data, confirming that SNPE did a good job at capturing the observed data in this simple case. References \u00b6 A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117(4):500\u2013544, 1952. M. Pospischil, M. Toledo-Rodriguez, C. Monier, Z. Piwkowska, T. Bal, Y. Fr\u00e9gnac, H. Markram, and A. Destexhe. Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biological Cybernetics, 99(4-5), 2008.","title":"Hodgkin-Huxley example"},{"location":"examples/00_HH_simulator/#inference-on-hodgkin-huxley-model-tutorial","text":"In this tutorial, we use sbi to do inference on a Hodgkin-Huxley model from neuroscience (Hodgkin and Huxley, 1952). We will learn two parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ) based on a current-clamp recording, that we generate synthetically (in practice, this would be an experimental observation). Note, you find the original version of this notebook at https://github.com/sbi-dev/sbi/blob/main/examples/00_HH_simulator.ipynb in the sbi repository. First we are going to import basic packages. import numpy as np import torch # visualization import matplotlib as mpl import matplotlib.pyplot as plt # sbi from sbi import utils as utils from sbi import analysis as analysis from sbi.inference.base import infer # remove top and right axis from plots mpl . rcParams [ \"axes.spines.right\" ] = False mpl . rcParams [ \"axes.spines.top\" ] = False","title":"Inference on Hodgkin-Huxley model: tutorial"},{"location":"examples/00_HH_simulator/#different-required-components","text":"Before running inference, let us define the different required components: observed data prior over model parameters simulator","title":"Different required components"},{"location":"examples/00_HH_simulator/#1-observed-data","text":"Let us assume we current-clamped a neuron and recorded the following voltage trace: In fact, this voltage trace was not measured experimentally but synthetically generated by simulating a Hodgkin-Huxley model with particular parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ). We will come back to this point later in the tutorial.","title":"1. Observed data"},{"location":"examples/00_HH_simulator/#2-simulator","text":"We would like to infer the posterior over the two parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ) of a Hodgkin-Huxley model, given the observed electrophysiological recording above. The model has channel kinetics as in Pospischil et al. 2008 , and is defined by the following set of differential equations (parameters of interest highlighted in orange): \\[ \\scriptsize \\begin{align} C_m\\frac{dV}{dt}&=g_1\\left(E_1-V\\right)+ \\color{orange}{\\bar{g}_{Na}}m^3h\\left(E_{Na}-V\\right)+ \\color{orange}{\\bar{g}_{K}}n^4\\left(E_K-V\\right)+ \\bar{g}_Mp\\left(E_K-V\\right)+ I_{inj}+ \\sigma\\eta\\left(t\\right)\\\\ \\frac{dq}{dt}&=\\frac{q_\\infty\\left(V\\right)-q}{\\tau_q\\left(V\\right)},\\;q\\in\\{m,h,n,p\\} \\end{align} \\] Above, \\(V\\) represents the membrane potential, \\(C_m\\) is the membrane capacitance, \\(g_{\\text{l}}\\) is the leak conductance, \\(E_{\\text{l}}\\) is the membrane reversal potential, \\(\\bar{g}_c\\) is the density of channels of type \\(c\\) ( \\(\\text{Na}^+\\) , \\(\\text{K}^+\\) , M), \\(E_c\\) is the reversal potential of \\(c\\) , ( \\(m\\) , \\(h\\) , \\(n\\) , \\(p\\) ) are the respective channel gating kinetic variables, and \\(\\sigma \\eta(t)\\) is the intrinsic neural noise. The right hand side of the voltage dynamics is composed of a leak current, a voltage-dependent \\(\\text{Na}^+\\) current, a delayed-rectifier \\(\\text{K}^+\\) current, a slow voltage-dependent \\(\\text{K}^+\\) current responsible for spike-frequency adaptation, and an injected current \\(I_{\\text{inj}}\\) . Channel gating variables \\(q\\) have dynamics fully characterized by the neuron membrane potential \\(V\\) , given the respective steady-state \\(q_{\\infty}(V)\\) and time constant \\(\\tau_{q}(V)\\) (details in Pospischil et al. 2008). The input current \\(I_{\\text{inj}}\\) is defined as from HH_helper_functions import syn_current I , t_on , t_off , dt , t , A_soma = syn_current () The Hodgkin-Huxley simulator is given by: from HH_helper_functions import HHsimulator Putting the input current and the simulator together: def run_HH_model ( params ): params = np . asarray ( params ) # input current, time step I , t_on , t_off , dt , t , A_soma = syn_current () t = np . arange ( 0 , len ( I ), 1 ) * dt # initial voltage V0 = - 70 states = HHsimulator ( V0 , params . reshape ( 1 , - 1 ), dt , t , I ) return dict ( data = states . reshape ( - 1 ), time = t , dt = dt , I = I . reshape ( - 1 )) To get an idea of the output of the Hodgkin-Huxley model, let us generate some voltage traces for different parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), given the input current \\(I_{\\text{inj}}\\) : # three sets of (g_Na, g_K) params = np . array ([[ 50.0 , 1.0 ], [ 4.0 , 1.5 ], [ 20.0 , 15.0 ]]) num_samples = len ( params [:, 0 ]) sim_samples = np . zeros (( num_samples , len ( I ))) for i in range ( num_samples ): sim_samples [ i , :] = run_HH_model ( params = params [ i , :])[ \"data\" ] # colors for traces col_min = 2 num_colors = num_samples + col_min cm1 = mpl . cm . Blues col1 = [ cm1 ( 1.0 * i / num_colors ) for i in range ( col_min , num_colors )] fig = plt . figure ( figsize = ( 7 , 5 )) gs = mpl . gridspec . GridSpec ( 2 , 1 , height_ratios = [ 4 , 1 ]) ax = plt . subplot ( gs [ 0 ]) for i in range ( num_samples ): plt . plot ( t , sim_samples [ i , :], color = col1 [ i ], lw = 2 ) plt . ylabel ( \"voltage (mV)\" ) ax . set_xticks ([]) ax . set_yticks ([ - 80 , - 20 , 40 ]) ax = plt . subplot ( gs [ 1 ]) plt . plot ( t , I * A_soma * 1e3 , \"k\" , lw = 2 ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"input (nA)\" ) ax . set_xticks ([ 0 , max ( t ) / 2 , max ( t )]) ax . set_yticks ([ 0 , 1.1 * np . max ( I * A_soma * 1e3 )]) ax . yaxis . set_major_formatter ( mpl . ticker . FormatStrFormatter ( \" %.2f \" )) plt . show () As can be seen, the voltage traces can be quite diverse for different parameter values. Often, we are not interested in matching the exact trace, but only in matching certain features thereof. In this example of the Hodgkin-Huxley model, the summary features are the number of spikes, the mean resting potential, the standard deviation of the resting potential, and the first four voltage moments: mean, standard deviation, skewness and kurtosis. Using the function calculate_summary_statistics() imported below, we obtain these statistics from the output of the Hodgkin Huxley simulator. from HH_helper_functions import calculate_summary_statistics Lastly, we define a function that performs all of the above steps at once. The function simulation_wrapper takes in conductance values, runs the Hodgkin Huxley model and then returns the summary statistics. def simulation_wrapper ( params ): \"\"\" Returns summary statistics from conductance values in `params`. Summarizes the output of the HH simulator and converts it to `torch.Tensor`. \"\"\" obs = run_HH_model ( params ) summstats = torch . as_tensor ( calculate_summary_statistics ( obs )) return summstats sbi takes any function as simulator. Thus, sbi also has the flexibility to use simulators that utilize external packages, e.g., Brian ( http://briansimulator.org/ ), nest ( https://www.nest-simulator.org/ ), or NEURON ( https://neuron.yale.edu/neuron/ ). External simulators do not even need to be Python-based as long as they store simulation outputs in a format that can be read from Python. All that is necessary is to wrap your external simulator of choice into a Python callable that takes a parameter set and outputs a set of summary statistics we want to fit the parameters to.","title":"2. Simulator"},{"location":"examples/00_HH_simulator/#3-prior-over-model-parameters","text":"Now that we have the simulator, we need to define a function with the prior over the model parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), which in this case is chosen to be a Uniform distribution: prior_min = [ 0.5 , 1e-4 ] prior_max = [ 80.0 , 15.0 ] prior = utils . torchutils . BoxUniform ( low = torch . as_tensor ( prior_min ), high = torch . as_tensor ( prior_max ) )","title":"3. Prior over model parameters"},{"location":"examples/00_HH_simulator/#inference","text":"Now that we have all the required components, we can run inference with SNPE to identify parameters whose activity matches this trace. posterior = infer ( simulation_wrapper , prior , method = \"SNPE\" , num_simulations = 300 , num_workers = 4 ) HBox(children=(FloatProgress(value=0.0, description='Running 300 simulations in 300 batches.', max=300.0, styl\u2026 Neural network successfully converged after 233 epochs. Note sbi can parallelize your simulator. If you experience problems with parallelization, try setting num_workers=1 and please give us an error report as a GitHub issue .","title":"Inference"},{"location":"examples/00_HH_simulator/#coming-back-to-the-observed-data","text":"As mentioned at the beginning of the tutorial, the observed data are generated by the Hodgkin-Huxley model with a set of known parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ). To illustrate how to compute the summary statistics of the observed data, let us regenerate the observed data: # true parameters and respective labels true_params = np . array ([ 50.0 , 5.0 ]) labels_params = [ r \"$g_ {Na} $\" , r \"$g_ {K} $\" ] observation_trace = run_HH_model ( true_params ) observation_summary_statistics = calculate_summary_statistics ( observation_trace ) As we already shown above, the observed voltage traces look as follows: fig = plt . figure ( figsize = ( 7 , 5 )) gs = mpl . gridspec . GridSpec ( 2 , 1 , height_ratios = [ 4 , 1 ]) ax = plt . subplot ( gs [ 0 ]) plt . plot ( observation_trace [ \"time\" ], observation_trace [ \"data\" ]) plt . ylabel ( \"voltage (mV)\" ) plt . title ( \"observed data\" ) plt . setp ( ax , xticks = [], yticks = [ - 80 , - 20 , 40 ]) ax = plt . subplot ( gs [ 1 ]) plt . plot ( observation_trace [ \"time\" ], I * A_soma * 1e3 , \"k\" , lw = 2 ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"input (nA)\" ) ax . set_xticks ([ 0 , max ( observation_trace [ \"time\" ]) / 2 , max ( observation_trace [ \"time\" ])]) ax . set_yticks ([ 0 , 1.1 * np . max ( I * A_soma * 1e3 )]) ax . yaxis . set_major_formatter ( mpl . ticker . FormatStrFormatter ( \" %.2f \" ))","title":"Coming back to the observed data"},{"location":"examples/00_HH_simulator/#analysis-of-the-posterior-given-the-observed-data","text":"After running the inference algorithm, let us inspect the inferred posterior distribution over the parameters ( \\(\\bar g_{Na}\\) , \\(\\bar g_K\\) ), given the observed trace. To do so, we first draw samples (i.e. consistent parameter sets) from the posterior: samples = posterior . sample (( 10000 ,), x = observation_summary_statistics ) HBox(children=(FloatProgress(value=0.0, description='Drawing 10000 posterior samples', max=10000.0, style=Prog\u2026 fig , axes = analysis . pairplot ( samples , limits = [[ 0.5 , 80 ], [ 1e-4 , 15.0 ]], ticks = [[ 0.5 , 80 ], [ 1e-4 , 15.0 ]], figsize = ( 5 , 5 ), points = true_params , points_offdiag = { \"markersize\" : 6 }, points_colors = \"r\" , ); As can be seen, the inferred posterior contains the ground-truth parameters (red) in a high-probability region. Now, let us sample parameters from the posterior distribution, simulate the Hodgkin-Huxley model for this parameter set and compare the simulations with the observed data: # Draw a sample from the posterior and convert to numpy for plotting. posterior_sample = posterior . sample (( 1 ,), x = observation_summary_statistics ) . numpy () HBox(children=(FloatProgress(value=0.0, description='Drawing 1 posterior samples', max=1.0, style=ProgressStyl\u2026 fig = plt . figure ( figsize = ( 7 , 5 )) # plot observation t = observation_trace [ \"time\" ] y_obs = observation_trace [ \"data\" ] plt . plot ( t , y_obs , lw = 2 , label = \"observation\" ) # simulate and plot samples x = run_HH_model ( posterior_sample ) plt . plot ( t , x [ \"data\" ], \"--\" , lw = 2 , label = \"posterior sample\" ) plt . xlabel ( \"time (ms)\" ) plt . ylabel ( \"voltage (mV)\" ) ax = plt . gca () handles , labels = ax . get_legend_handles_labels () ax . legend ( handles [:: - 1 ], labels [:: - 1 ], bbox_to_anchor = ( 1.3 , 1 ), loc = \"upper right\" ) ax . set_xticks ([ 0 , 60 , 120 ]) ax . set_yticks ([ - 80 , - 20 , 40 ]); As can be seen, the sample from the inferred posterior leads to simulations that closely resemble the observed data, confirming that SNPE did a good job at capturing the observed data in this simple case.","title":"Analysis of the posterior given the observed data"},{"location":"examples/00_HH_simulator/#references","text":"A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117(4):500\u2013544, 1952. M. Pospischil, M. Toledo-Rodriguez, C. Monier, Z. Piwkowska, T. Bal, Y. Fr\u00e9gnac, H. Markram, and A. Destexhe. Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biological Cybernetics, 99(4-5), 2008.","title":"References"},{"location":"examples/01_decision_making_model/","text":"SBI for decision-making models \u00b6 In a previous tutorial , we showed how to use SBI with trial-based iid data. Such scenarios can arise, for example, in models of perceptual decision making. In addition to trial-based iid data points, these models often come with mixed data types and varying experimental conditions. Here, we show how sbi can be used to perform inference in such models with the MNLE method. Trial-based SBI with mixed data types \u00b6 In some cases, models with trial-based data additionally return data with mixed data types, e.g., continous and discrete data. For example, most computational models of decision-making have continuous reaction times and discrete choices as output. This can induce a problem when performing trial-based SBI that relies on learning a neural likelihood: It is challenging for most density estimators to handle both, continuous and discrete data at the same time. However, there is a recent SBI method for solving this problem, it\u2019s called Mixed Neural Likelihood Estimation (MNLE). It works just like NLE, but with mixed data types. The trick is that it learns two separate density estimators, one for the discrete part of the data, and one for the continuous part, and combines the two to obtain the final neural likelihood. Crucially, the continuous density estimator is trained conditioned on the output of the discrete one, such that statistical dependencies between the discrete and continuous data (e.g., between choices and reaction times) are modeled as well. The interested reader is referred to the original paper available here . MNLE was recently added to sbi (see this PR and also issue ) and follows the same API as SNLE . In this tutorial we will show how to apply MNLE to mixed data, and how to deal with varying experimental conditions. Toy problem for MNLE \u00b6 To illustrate MNLE we set up a toy simulator that outputs mixed data and for which we know the likelihood such we can obtain reference posterior samples via MCMC. Simulator : To simulate mixed data we do the following Sample reaction time from inverse Gamma Sample choices from Binomial Return reaction time \\(rt \\in (0, \\infty)\\) and choice index \\(c \\in \\{0, 1\\}\\) \\[ c \\sim \\text{Binomial}(\\rho) \\\\ rt \\sim \\text{InverseGamma}(\\alpha=2, \\beta) \\\\ \\] Prior : The priors of the two parameters \\(\\rho\\) and \\(\\beta\\) are independent. We define a Beta prior over the probabilty parameter of the Binomial used in the simulator and a Gamma prior over the shape-parameter of the inverse Gamma used in the simulator: \\[ p(\\beta, \\rho) = p(\\beta) \\; p(\\rho) ; \\\\ p(\\beta) = \\text{Gamma}(1, 0.5) \\\\ p(\\text{probs}) = \\text{Beta}(2, 2) \\] Because the InverseGamma and the Binomial likelihoods are well-defined we can perform MCMC on this problem and obtain reference-posterior samples. import matplotlib.pyplot as plt import torch from torch import Tensor from sbi.inference import MNLE from pyro.distributions import InverseGamma from torch.distributions import Beta , Binomial , Categorical , Gamma from sbi.utils import MultipleIndependent from sbi.utils.metrics import c2st from sbi.analysis import pairplot from sbi.inference import MCMCPosterior from sbi.utils.torchutils import atleast_2d from sbi.inference.potentials.likelihood_based_potential import ( MixedLikelihoodBasedPotential , ) from sbi.utils.conditional_density_utils import ConditionedPotential from sbi.utils import mcmc_transform from sbi.inference.potentials.base_potential import BasePotential # Toy simulator for mixed data def mixed_simulator ( theta : Tensor , concentration_scaling : float = 1.0 ): \"\"\"Returns a sample from a mixed distribution given parameters theta. Args: theta: batch of parameters, shape (batch_size, 2) concentration_scaling: scaling factor for the concentration parameter of the InverseGamma distribution, mimics an experimental condition. \"\"\" beta , ps = theta [:, : 1 ], theta [:, 1 :] choices = Binomial ( probs = ps ) . sample () rts = InverseGamma ( concentration = concentration_scaling * torch . ones_like ( beta ), rate = beta ) . sample () return torch . cat (( rts , choices ), dim = 1 ) # The potential function defines the ground truth likelihood and allows us to obtain reference posterior samples via MCMC. class PotentialFunctionProvider ( BasePotential ): allow_iid_x = True # type: ignore def __init__ ( self , prior , x_o , concentration_scaling = 1.0 , device = \"cpu\" ): super () . __init__ ( prior , x_o , device ) self . concentration_scaling = concentration_scaling def __call__ ( self , theta , track_gradients : bool = True ): theta = atleast_2d ( theta ) with torch . set_grad_enabled ( track_gradients ): iid_ll = self . iid_likelihood ( theta ) return iid_ll + self . prior . log_prob ( theta ) def iid_likelihood ( self , theta ): lp_choices = torch . stack ( [ Binomial ( probs = th . reshape ( 1 , - 1 )) . log_prob ( self . x_o [:, 1 :]) for th in theta [:, 1 :] ], dim = 1 , ) lp_rts = torch . stack ( [ InverseGamma ( concentration = self . concentration_scaling * torch . ones_like ( beta_i ), rate = beta_i , ) . log_prob ( self . x_o [:, : 1 ]) for beta_i in theta [:, : 1 ] ], dim = 1 , ) joint_likelihood = ( lp_choices + lp_rts ) . squeeze () assert joint_likelihood . shape == torch . Size ([ self . x_o . shape [ 0 ], theta . shape [ 0 ]]) return joint_likelihood . sum ( 0 ) # Define independent prior. prior = MultipleIndependent ( [ Gamma ( torch . tensor ([ 1.0 ]), torch . tensor ([ 0.5 ])), Beta ( torch . tensor ([ 2.0 ]), torch . tensor ([ 2.0 ])), ], validate_args = False , ) Obtain reference-posterior samples via analytical likelihood and MCMC \u00b6 torch . manual_seed ( 42 ) num_trials = 10 num_samples = 1000 theta_o = prior . sample (( 1 ,)) x_o = mixed_simulator ( theta_o . repeat ( num_trials , 1 )) mcmc_kwargs = dict ( num_chains = 20 , warmup_steps = 50 , method = \"slice_np_vectorized\" , init_strategy = \"proposal\" , ) true_posterior = MCMCPosterior ( potential_fn = PotentialFunctionProvider ( prior , x_o ), proposal = prior , theta_transform = mcmc_transform ( prior , enable_transform = True ), ** mcmc_kwargs , ) true_samples = true_posterior . sample (( num_samples ,)) /Users/janbolts/qode/sbi/sbi/utils/sbiutils.py:342: UserWarning: An x with a batch size of 10 was passed. It will be interpreted as a batch of independent and identically distributed data X={x_1, ..., x_n}, i.e., data generated based on the same underlying (unknown) parameter. The resulting posterior will be with respect to entire batch, i.e,. p(theta | X). warnings.warn( Running vectorized MCMC with 20 chains: 0%| | 0/20000 [00:00 1 , you might experience an error that a certain object from your simulator could not be pickled (an example can be found here ). This can be fixed by forcing sbi to pickle with dill instead of the default cloudpickle . To do so, adjust your code as follows: Install dill : pip install dill At the very beginning of your python script, set the pickler to dill : from joblib.externals.loky import set_loky_pickler set_loky_pickler ( \"dill\" ) Move all imports required by your simulator into the simulator: # Imports specified outside of the simulator will break dill: import torch def my_simulator ( parameters ): return torch . ones ( 1 , 10 ) # Therefore, move the imports into the simulator: def my_simulator ( parameters ): import torch return torch . ones ( 1 , 10 ) Alternative: parallelize yourself \u00b6 You can also write your own code to parallelize simulations with whatever multiprocessing framework you prefer. You can then simulate your data outside of sbi and pass the simulated data as shown in the flexible interface : Some more background \u00b6 sbi uses joblib to parallelize simulations, which in turn uses pickle or cloudpickle to serialize the simulator. Almost all simulators will be picklable with cloudpickle , but we have experienced issues e.g. with neuron simulators, see here .","title":"When using multiple workers, I get a pickling error. Can I still use multiprocessing?"},{"location":"faq/question_03/#when-using-multiple-workers-i-get-a-pickling-error-can-i-still-use-multiprocessing","text":"Yes, but you will have to make a few adjustments to your code. Some background: when using num_workers > 1 , you might experience an error that a certain object from your simulator could not be pickled (an example can be found here ). This can be fixed by forcing sbi to pickle with dill instead of the default cloudpickle . To do so, adjust your code as follows: Install dill : pip install dill At the very beginning of your python script, set the pickler to dill : from joblib.externals.loky import set_loky_pickler set_loky_pickler ( \"dill\" ) Move all imports required by your simulator into the simulator: # Imports specified outside of the simulator will break dill: import torch def my_simulator ( parameters ): return torch . ones ( 1 , 10 ) # Therefore, move the imports into the simulator: def my_simulator ( parameters ): import torch return torch . ones ( 1 , 10 )","title":"When using multiple workers, I get a pickling error. Can I still use multiprocessing?"},{"location":"faq/question_03/#alternative-parallelize-yourself","text":"You can also write your own code to parallelize simulations with whatever multiprocessing framework you prefer. You can then simulate your data outside of sbi and pass the simulated data as shown in the flexible interface :","title":"Alternative: parallelize yourself"},{"location":"faq/question_03/#some-more-background","text":"sbi uses joblib to parallelize simulations, which in turn uses pickle or cloudpickle to serialize the simulator. Almost all simulators will be picklable with cloudpickle , but we have experienced issues e.g. with neuron simulators, see here .","title":"Some more background"},{"location":"faq/question_04/","text":"Can I use the GPU for training the density estimator? \u00b6 TLDR; Yes, by passing device=\"cuda\" and by passing a prior that lives on the device name your passed. But no speed-ups for default density estimators. Yes. When creating the inference object in the flexible interface, you can pass the device as an argument, e.g., inference = SNPE ( prior , device = \"cuda\" , density_estimator = \"maf\" ) The device is set to \"cpu\" by default, and it can be set to anything, as long as it maps to an existing PyTorch CUDA device. sbi will take care of copying the net and the training data to and from the device . Note that the prior must be on the training device already, e.g., when passing device=\"cuda:0\" , make sure to pass a prior object that was created on that device, e.g., prior = torch.distributions.MultivariateNormal(loc=torch.zeros(2, device=\"cuda:0\"), covariance_matrix=torch.eye(2, device=\"cuda:0\")) . Performance \u00b6 Whether or not you reduce your training time when training on a GPU depends on the problem at hand. We provide a couple of default density estimators for SNPE , SNLE and SNRE , e.g., a mixture density network ( density_estimator=\"mdn\" ) or a Masked Autoregressive Flow ( density_estimator=\"maf\" ). For those default density estimators we do not expect a speed up. This is because the underlying neural networks are quite shallow and not tall, e.g., they do not have many parameters or matrix operations that profit a lot from being executed on the GPU. A speed up through training on the GPU will most likely become visible when you are using convolutional modules in your neural networks. E.g., when passing an embedding net for image processing like in this example: https://github.com/sbi-dev/sbi/blob/main/tutorials/05_embedding_net.ipynb .","title":"Can I use the GPU for training the density estimator?"},{"location":"faq/question_04/#can-i-use-the-gpu-for-training-the-density-estimator","text":"TLDR; Yes, by passing device=\"cuda\" and by passing a prior that lives on the device name your passed. But no speed-ups for default density estimators. Yes. When creating the inference object in the flexible interface, you can pass the device as an argument, e.g., inference = SNPE ( prior , device = \"cuda\" , density_estimator = \"maf\" ) The device is set to \"cpu\" by default, and it can be set to anything, as long as it maps to an existing PyTorch CUDA device. sbi will take care of copying the net and the training data to and from the device . Note that the prior must be on the training device already, e.g., when passing device=\"cuda:0\" , make sure to pass a prior object that was created on that device, e.g., prior = torch.distributions.MultivariateNormal(loc=torch.zeros(2, device=\"cuda:0\"), covariance_matrix=torch.eye(2, device=\"cuda:0\")) .","title":"Can I use the GPU for training the density estimator?"},{"location":"faq/question_04/#performance","text":"Whether or not you reduce your training time when training on a GPU depends on the problem at hand. We provide a couple of default density estimators for SNPE , SNLE and SNRE , e.g., a mixture density network ( density_estimator=\"mdn\" ) or a Masked Autoregressive Flow ( density_estimator=\"maf\" ). For those default density estimators we do not expect a speed up. This is because the underlying neural networks are quite shallow and not tall, e.g., they do not have many parameters or matrix operations that profit a lot from being executed on the GPU. A speed up through training on the GPU will most likely become visible when you are using convolutional modules in your neural networks. E.g., when passing an embedding net for image processing like in this example: https://github.com/sbi-dev/sbi/blob/main/tutorials/05_embedding_net.ipynb .","title":"Performance"},{"location":"faq/question_05/","text":"How should I save and load objects in sbi ? \u00b6 NeuralPosterior objects are picklable. import pickle # ... run inference posterior = inference . build_posterior () with open ( \"/path/to/my_posterior.pkl\" , \"wb\" ) as handle : pickle . dump ( posterior , handle ) Note: posterior objects that were saved under sbi v0.17.2 or older can not be loaded under sbi v0.18.0 or newer. Note: if you try to load a posterior that was saved under sbi v0.14.x or earlier under sbi v0.15.x until sbi v0.17.x , you have to add: import sys from sbi.utils import user_input_checks_utils sys . modules [ \"sbi.user_input.user_input_checks_utils\" ] = user_input_checks_utils to your script before loading the posterior. As of sbi v0.18.0 , NeuralInference objects are also picklable. import pickle # ... run inference posterior = inference . build_posterior () with open ( \"/path/to/my_inference.pkl\" , \"wb\" ) as handle : pickle . dump ( inference , handle ) However, saving and loading the inference object will slightly modify the object (in order to make it serializable). These modifications lead to the following two changes in behaviour: 1) Retraining from scratch is not supported, i.e. .train(..., retrain_from_scratch=True) does not work. 2) When the loaded object calls the .train() method, it generates a new tensorboard summary writer (instead of appending to the current one).","title":"How should I save and load objects in sbi?"},{"location":"faq/question_05/#how-should-i-save-and-load-objects-in-sbi","text":"NeuralPosterior objects are picklable. import pickle # ... run inference posterior = inference . build_posterior () with open ( \"/path/to/my_posterior.pkl\" , \"wb\" ) as handle : pickle . dump ( posterior , handle ) Note: posterior objects that were saved under sbi v0.17.2 or older can not be loaded under sbi v0.18.0 or newer. Note: if you try to load a posterior that was saved under sbi v0.14.x or earlier under sbi v0.15.x until sbi v0.17.x , you have to add: import sys from sbi.utils import user_input_checks_utils sys . modules [ \"sbi.user_input.user_input_checks_utils\" ] = user_input_checks_utils to your script before loading the posterior. As of sbi v0.18.0 , NeuralInference objects are also picklable. import pickle # ... run inference posterior = inference . build_posterior () with open ( \"/path/to/my_inference.pkl\" , \"wb\" ) as handle : pickle . dump ( inference , handle ) However, saving and loading the inference object will slightly modify the object (in order to make it serializable). These modifications lead to the following two changes in behaviour: 1) Retraining from scratch is not supported, i.e. .train(..., retrain_from_scratch=True) does not work. 2) When the loaded object calls the .train() method, it generates a new tensorboard summary writer (instead of appending to the current one).","title":"How should I save and load objects in sbi?"},{"location":"faq/question_06/","text":"Can I stop neural network training and resume it later? \u00b6 Many clusters have a time limit and sbi might exceed this limit. You can circumvent this problem by using the flexible interface . After simulations are finished, sbi trains a neural network. If this process takes too long, you can stop training and resume it later. The syntax is: inference = SNPE ( prior = prior ) inference = inference . append_simulations ( theta , x ) inference . train ( max_num_epochs = 300 ) # Pick `max_num_epochs` such that it does not exceed the runtime. with open ( \"path/to/my/inference.pkl\" , \"wb\" ) as handle : pickle . dump ( inference , handle ) # To resume training: with open ( \"path/to/my/inference.pkl\" , \"rb\" ) as handle : inference_from_disk = pickle . load ( handle ) inference_from_disk . train ( resume_training = True , max_num_epochs = 600 ) # Run epochs 301 until 600 (or stop early). posterior = inference_from_disk . build_posterior ()","title":"Can I stop neural network training and resume it later?"},{"location":"faq/question_06/#can-i-stop-neural-network-training-and-resume-it-later","text":"Many clusters have a time limit and sbi might exceed this limit. You can circumvent this problem by using the flexible interface . After simulations are finished, sbi trains a neural network. If this process takes too long, you can stop training and resume it later. The syntax is: inference = SNPE ( prior = prior ) inference = inference . append_simulations ( theta , x ) inference . train ( max_num_epochs = 300 ) # Pick `max_num_epochs` such that it does not exceed the runtime. with open ( \"path/to/my/inference.pkl\" , \"wb\" ) as handle : pickle . dump ( inference , handle ) # To resume training: with open ( \"path/to/my/inference.pkl\" , \"rb\" ) as handle : inference_from_disk = pickle . load ( handle ) inference_from_disk . train ( resume_training = True , max_num_epochs = 600 ) # Run epochs 301 until 600 (or stop early). posterior = inference_from_disk . build_posterior ()","title":"Can I stop neural network training and resume it later?"},{"location":"faq/question_07/","text":"Can I use a custom prior with sbi? \u00b6 sbi works with torch distributions only so we recommend to use those whenever possible. For example, when you are used to using scipy.stats distributions as priors then we recommend using the corresponding torch.distributions , most common distributions are implemented there. In case you want to use a custom prior that is not in the set of common distributions that\u2019s possible as well: You need to write a prior class that mimicks the behaviour of a torch.distributions.Distribution class. Then sbi will wrap this class to make it a fully functional torch Distribution . Essentially, the class needs two methods: .sample(sample_shape) , where sample_shape is a shape tuple, e.g., (n,) , and returns a batch of n samples, e.g., of shape (n, 2)` for a two dimenional prior. .log_prob(value) method that returns the \u201clog probs\u201d of parameters under the prior, e.g., for a batches of n parameters with shape (n, ndims) it should return a log probs array of shape (n,) . For sbi > 0.17.2 this could look like the following: class CustomUniformPrior : \"\"\"User defined numpy uniform prior. Custom prior with user-defined valid .sample and .log_prob methods. \"\"\" def __init__ ( self , lower : Tensor , upper : Tensor , return_numpy : bool = False ): self . lower = lower self . upper = upper self . dist = BoxUniform ( lower , upper ) self . return_numpy = return_numpy def sample ( self , sample_shape = torch . Size ([])): samples = self . dist . sample ( sample_shape ) return samples . numpy () if self . return_numpy else samples def log_prob ( self , values ): if self . return_numpy : values = torch . as_tensor ( values ) log_probs = self . dist . log_prob ( values ) return log_probs . numpy () if self . return_numpy else log_probs Once you have such a class you can wrap into a Distribution using the process_prior function sbi provides: from sbi.utils import process_prior custom_prior = CustomUniformPrior ( torch . zeros ( 2 ), torch . ones ( 2 )) prior , * _ = process_prior ( custom_prior ) # Keeping only the first return. # use this wrapped prior in sbi... In sbi it is sometimes required to check the support of the prior, e.g., when the prior support is bounded and one wants to reject samples from the posterior density estimator that lie outside the prior support. In torch Distributions this is handled automatically, however, when using a custom prior it is not. Thus, if your prior has bounded support (like the one above) it makes sense to pass the bounds to the wrapper function such that sbi can pass them to torch Distributions : from sbi.utils import process_prior custom_prior = CustomUniformPrior ( torch . zeros ( 2 ), torch . ones ( 2 )) prior = process_prior ( custom_prior , custom_prior_wrapper_kwargs = dict ( lower_bound = torch . zeros ( 2 ), upper_bound = torch . ones ( 2 ))) # use this wrapped prior in sbi... Note that in custom_prior_wrapper_kwargs you can pass additinal arguments for the wrapper, e.g., validate_args or arg_constraints see the Distribution documentation for more details. If you are running sbi < 0.17.2 and use SNLE the code above will produce a NotImplementedError (see #581 ). In this case you need to update to a newer version of sbi or use SNPE instead.","title":"Can I use a custom prior with sbi?"},{"location":"faq/question_07/#can-i-use-a-custom-prior-with-sbi","text":"sbi works with torch distributions only so we recommend to use those whenever possible. For example, when you are used to using scipy.stats distributions as priors then we recommend using the corresponding torch.distributions , most common distributions are implemented there. In case you want to use a custom prior that is not in the set of common distributions that\u2019s possible as well: You need to write a prior class that mimicks the behaviour of a torch.distributions.Distribution class. Then sbi will wrap this class to make it a fully functional torch Distribution . Essentially, the class needs two methods: .sample(sample_shape) , where sample_shape is a shape tuple, e.g., (n,) , and returns a batch of n samples, e.g., of shape (n, 2)` for a two dimenional prior. .log_prob(value) method that returns the \u201clog probs\u201d of parameters under the prior, e.g., for a batches of n parameters with shape (n, ndims) it should return a log probs array of shape (n,) . For sbi > 0.17.2 this could look like the following: class CustomUniformPrior : \"\"\"User defined numpy uniform prior. Custom prior with user-defined valid .sample and .log_prob methods. \"\"\" def __init__ ( self , lower : Tensor , upper : Tensor , return_numpy : bool = False ): self . lower = lower self . upper = upper self . dist = BoxUniform ( lower , upper ) self . return_numpy = return_numpy def sample ( self , sample_shape = torch . Size ([])): samples = self . dist . sample ( sample_shape ) return samples . numpy () if self . return_numpy else samples def log_prob ( self , values ): if self . return_numpy : values = torch . as_tensor ( values ) log_probs = self . dist . log_prob ( values ) return log_probs . numpy () if self . return_numpy else log_probs Once you have such a class you can wrap into a Distribution using the process_prior function sbi provides: from sbi.utils import process_prior custom_prior = CustomUniformPrior ( torch . zeros ( 2 ), torch . ones ( 2 )) prior , * _ = process_prior ( custom_prior ) # Keeping only the first return. # use this wrapped prior in sbi... In sbi it is sometimes required to check the support of the prior, e.g., when the prior support is bounded and one wants to reject samples from the posterior density estimator that lie outside the prior support. In torch Distributions this is handled automatically, however, when using a custom prior it is not. Thus, if your prior has bounded support (like the one above) it makes sense to pass the bounds to the wrapper function such that sbi can pass them to torch Distributions : from sbi.utils import process_prior custom_prior = CustomUniformPrior ( torch . zeros ( 2 ), torch . ones ( 2 )) prior = process_prior ( custom_prior , custom_prior_wrapper_kwargs = dict ( lower_bound = torch . zeros ( 2 ), upper_bound = torch . ones ( 2 ))) # use this wrapped prior in sbi... Note that in custom_prior_wrapper_kwargs you can pass additinal arguments for the wrapper, e.g., validate_args or arg_constraints see the Distribution documentation for more details. If you are running sbi < 0.17.2 and use SNLE the code above will produce a NotImplementedError (see #581 ). In this case you need to update to a newer version of sbi or use SNPE instead.","title":"Can I use a custom prior with sbi?"},{"location":"tutorial/00_getting_started/","text":"Getting started with sbi \u00b6 Note, you can find the original version of this notebook at https://github.com/sbi-dev/sbi/blob/main/tutorials/00_getting_started.ipynb in the sbi repository. import torch from sbi import utils as utils from sbi import analysis as analysis from sbi.inference.base import infer Running the inference procedure \u00b6 sbi provides a simple interface to run state-of-the-art algorithms for simulation-based inference. For inference, you need to provide two ingredients: 1) a prior distribution that allows to sample parameter sets. 2) a simulator that takes parameter sets and produces simulation outputs. For example, we can have a 3-dimensional parameter space with a uniform prior between [-1,1] and a simple simulator that for the sake of example adds 1.0 and some Gaussian noise to the parameter set: num_dim = 3 prior = utils . BoxUniform ( low =- 2 * torch . ones ( num_dim ), high = 2 * torch . ones ( num_dim )) def simulator ( parameter_set ): return 1.0 + parameter_set + torch . randn ( parameter_set . shape ) * 0.1 sbi can then run inference: # Other methods are \"SNLE\" or \"SNRE\". posterior = infer ( simulator , prior , method = \"SNPE\" , num_simulations = 1000 ) Running 1000 simulations.: 0%| | 0/1000 [00:001 , the posterior is no longer amortized: it will give good results when sampled around x=observation , but possibly bad results for other x . Once we have obtained the posterior, we can .sample() , .log_prob() , or .pairplot() in the same way as for the simple interface. posterior_samples = posterior . sample (( 10000 ,), x = x_o ) # plot posterior samples _ = analysis . pairplot ( posterior_samples , limits = [[ - 2 , 2 ], [ - 2 , 2 ], [ - 2 , 2 ]], figsize = ( 5 , 5 ) ) Drawing 10000 posterior samples: 0%| | 0/10000 [00:001 , the posterior is no longer amortized: it will give good results when sampled around x=observation , but possibly bad results for other x . Once we have obtained the posterior, we can .sample() , .log_prob() , or .pairplot() in the same way as for the simple interface. posterior_samples = posterior . sample (( 10000 ,), x = x_o ) # plot posterior samples _ = analysis . pairplot ( posterior_samples , limits = [[ - 2 , 2 ], [ - 2 , 2 ], [ - 2 , 2 ]], figsize = ( 5 , 5 ) ) Drawing 10000 posterior samples: 0%| | 0/10000 [00:00 The simulator model \u00b6 The simulator model that we consider has two parameters: \\(r\\) and \\(\\theta\\) . On each run, it generates 100 two-dimensional points centered around \\((r \\cos(\\theta), r \\sin(\\theta))\\) and perturbed by a Gaussian noise with variance 0.01. Instead of simply outputting the \\((x,y)\\) coordinates of each data point, the model generates a grayscale image of the scattered points with dimensions 32 by 32. This image is further perturbed by an uniform noise with values betweeen 0 and 0.2. The code below defines such model. def simulator_model ( parameter , return_points = False ): \"\"\"Simulator model with two-dimensional input parameter and 1024-dimensional output This simulator serves as a basic example for using a neural net for learning summary features. It has only two input parameters but generates high-dimensional output vectors. The data is generated as follows: (-) Input: parameter = [r, theta] (1) Generate 100 two-dimensional points centered around (r cos(theta),r sin(theta)) and perturbed by a Gaussian noise with variance 0.01 (2) Create a grayscale image I of the scattered points with dimensions 32 by 32 (3) Perturb I with an uniform noise with values betweeen 0 and 0.2 (-) Output: I Parameters ---------- parameter : array-like, shape (2) The two input parameters of the model, ordered as [r, theta] return_points : bool (default: False) Whether the simulator should return the coordinates of the simulated data points as well Returns ------- I: torch tensor, shape (1, 1024) Output flattened image (optional) points: array-like, shape (100, 2) Coordinates of the 2D simulated data points \"\"\" r = parameter [ 0 ] theta = parameter [ 1 ] sigma_points = 0.10 npoints = 100 points = [] for _ in range ( npoints ): x = r * torch . cos ( theta ) + sigma_points * torch . randn ( 1 ) y = r * torch . sin ( theta ) + sigma_points * torch . randn ( 1 ) points . append ([ x , y ]) points = torch . as_tensor ( points ) nx = 32 ny = 32 sigma_image = 0.20 I = torch . zeros ( nx , ny ) for point in points : pi = int (( point [ 0 ] - ( - 1 )) / (( + 1 ) - ( - 1 )) * nx ) pj = int (( point [ 1 ] - ( - 1 )) / (( + 1 ) - ( - 1 )) * ny ) if ( pi < nx ) and ( pj < ny ): I [ pi , pj ] = 1 I = I + sigma_image * torch . rand ( nx , ny ) I = I . T I = I . reshape ( 1 , - 1 ) if return_points : return I , points else : return I The figure below shows an example of the output of the simulator when \\(r = 0.70\\) and \\(\\theta = \\pi/4\\) # simulate samples true_parameter = torch . tensor ([ 0.70 , torch . pi / 4 ]) x_observed , x_points = simulator_model ( true_parameter , return_points = True ) # plot the observation fig , ax = plt . subplots ( facecolor = \"white\" , figsize = ( 11.15 , 5.61 ), ncols = 2 , constrained_layout = True ) circle = plt . Circle (( 0 , 0 ), 1.0 , color = \"k\" , ls = \"--\" , lw = 0.8 , fill = False ) ax [ 0 ] . add_artist ( circle ) ax [ 0 ] . scatter ( x_points [:, 0 ], x_points [:, 1 ], s = 20 ) ax [ 0 ] . set_xlabel ( \"x\" ) ax [ 0 ] . set_ylabel ( \"y\" ) ax [ 0 ] . set_xlim ( - 1 , + 1 ) ax [ 0 ] . set_xticks ([ - 1 , 0.0 , + 1.0 ]) ax [ 0 ] . set_ylim ( - 1 , + 1 ) ax [ 0 ] . set_yticks ([ - 1 , 0.0 , + 1.0 ]) ax [ 0 ] . set_title ( r \"original simulated points with $r = 0.70$ and $\\theta = \\pi/4$\" ) ax [ 1 ] . imshow ( x_observed . view ( 32 , 32 ), origin = \"lower\" , cmap = \"gray\" ) ax [ 1 ] . set_xticks ([]) ax [ 1 ] . set_yticks ([]) ax [ 1 ] . set_title ( \"noisy observed data (gray image with 32 x 32 pixels)\" ) Text(0.5, 1.0, 'noisy observed data (gray image with 32 x 32 pixels)') Defining an embedding_net \u00b6 An inference procedure applied to the output data from this simulator model determines the posterior distribution of \\(r\\) and \\(\\theta\\) given an observation of \\(x\\) , which lives in a 1024 dimensional space (32 x 32 = 1024). To avoid working directly on these high-dimensional vectors, one can use a convolutional neural network (CNN) that takes the 32x32 images as input and encodes them into 8-dimensional feature vectors. This CNN is trained along with the neural density estimator of the inference procedure and serves as an automatic summary statistics extractor. We define and instantiate the CNN as follows: class SummaryNet ( nn . Module ): def __init__ ( self ): super () . __init__ () # 2D convolutional layer self . conv1 = nn . Conv2d ( in_channels = 1 , out_channels = 6 , kernel_size = 5 , padding = 2 ) # Maxpool layer that reduces 32x32 image to 4x4 self . pool = nn . MaxPool2d ( kernel_size = 8 , stride = 8 ) # Fully connected layer taking as input the 6 flattened output arrays from the maxpooling layer self . fc = nn . Linear ( in_features = 6 * 4 * 4 , out_features = 8 ) def forward ( self , x ): x = x . view ( - 1 , 1 , 32 , 32 ) x = self . pool ( F . relu ( self . conv1 ( x ))) x = x . view ( - 1 , 6 * 4 * 4 ) x = F . relu ( self . fc ( x )) return x embedding_net = SummaryNet () The inference procedure \u00b6 With the embedding_net defined and instantiated, we can follow the usual workflow of an inference procedure in sbi . The embedding_net object appears as an input argument when instantiating the neural density estimator with utils.posterior_nn . # set prior distribution for the parameters prior = utils . BoxUniform ( low = torch . tensor ([ 0.0 , 0.0 ]), high = torch . tensor ([ 1.0 , 2 * torch . pi ]) ) # make a SBI-wrapper on the simulator object for compatibility simulator_wrapper , prior = prepare_for_sbi ( simulator_model , prior ) # instantiate the neural density estimator neural_posterior = utils . posterior_nn ( model = \"maf\" , embedding_net = embedding_net , hidden_features = 10 , num_transforms = 2 ) # setup the inference procedure with the SNPE-C procedure inference = SNPE ( prior = prior , density_estimator = neural_posterior ) # run the inference procedure on one round and 10000 simulated data points theta , x = simulate_for_sbi ( simulator_wrapper , prior , num_simulations = 10000 ) Running 10000 simulations.: 0%| | 0/10000 [00:00","title":"Learning summary statistics with a neural net"},{"location":"tutorial/05_embedding_net/#the-simulator-model","text":"The simulator model that we consider has two parameters: \\(r\\) and \\(\\theta\\) . On each run, it generates 100 two-dimensional points centered around \\((r \\cos(\\theta), r \\sin(\\theta))\\) and perturbed by a Gaussian noise with variance 0.01. Instead of simply outputting the \\((x,y)\\) coordinates of each data point, the model generates a grayscale image of the scattered points with dimensions 32 by 32. This image is further perturbed by an uniform noise with values betweeen 0 and 0.2. The code below defines such model. def simulator_model ( parameter , return_points = False ): \"\"\"Simulator model with two-dimensional input parameter and 1024-dimensional output This simulator serves as a basic example for using a neural net for learning summary features. It has only two input parameters but generates high-dimensional output vectors. The data is generated as follows: (-) Input: parameter = [r, theta] (1) Generate 100 two-dimensional points centered around (r cos(theta),r sin(theta)) and perturbed by a Gaussian noise with variance 0.01 (2) Create a grayscale image I of the scattered points with dimensions 32 by 32 (3) Perturb I with an uniform noise with values betweeen 0 and 0.2 (-) Output: I Parameters ---------- parameter : array-like, shape (2) The two input parameters of the model, ordered as [r, theta] return_points : bool (default: False) Whether the simulator should return the coordinates of the simulated data points as well Returns ------- I: torch tensor, shape (1, 1024) Output flattened image (optional) points: array-like, shape (100, 2) Coordinates of the 2D simulated data points \"\"\" r = parameter [ 0 ] theta = parameter [ 1 ] sigma_points = 0.10 npoints = 100 points = [] for _ in range ( npoints ): x = r * torch . cos ( theta ) + sigma_points * torch . randn ( 1 ) y = r * torch . sin ( theta ) + sigma_points * torch . randn ( 1 ) points . append ([ x , y ]) points = torch . as_tensor ( points ) nx = 32 ny = 32 sigma_image = 0.20 I = torch . zeros ( nx , ny ) for point in points : pi = int (( point [ 0 ] - ( - 1 )) / (( + 1 ) - ( - 1 )) * nx ) pj = int (( point [ 1 ] - ( - 1 )) / (( + 1 ) - ( - 1 )) * ny ) if ( pi < nx ) and ( pj < ny ): I [ pi , pj ] = 1 I = I + sigma_image * torch . rand ( nx , ny ) I = I . T I = I . reshape ( 1 , - 1 ) if return_points : return I , points else : return I The figure below shows an example of the output of the simulator when \\(r = 0.70\\) and \\(\\theta = \\pi/4\\) # simulate samples true_parameter = torch . tensor ([ 0.70 , torch . pi / 4 ]) x_observed , x_points = simulator_model ( true_parameter , return_points = True ) # plot the observation fig , ax = plt . subplots ( facecolor = \"white\" , figsize = ( 11.15 , 5.61 ), ncols = 2 , constrained_layout = True ) circle = plt . Circle (( 0 , 0 ), 1.0 , color = \"k\" , ls = \"--\" , lw = 0.8 , fill = False ) ax [ 0 ] . add_artist ( circle ) ax [ 0 ] . scatter ( x_points [:, 0 ], x_points [:, 1 ], s = 20 ) ax [ 0 ] . set_xlabel ( \"x\" ) ax [ 0 ] . set_ylabel ( \"y\" ) ax [ 0 ] . set_xlim ( - 1 , + 1 ) ax [ 0 ] . set_xticks ([ - 1 , 0.0 , + 1.0 ]) ax [ 0 ] . set_ylim ( - 1 , + 1 ) ax [ 0 ] . set_yticks ([ - 1 , 0.0 , + 1.0 ]) ax [ 0 ] . set_title ( r \"original simulated points with $r = 0.70$ and $\\theta = \\pi/4$\" ) ax [ 1 ] . imshow ( x_observed . view ( 32 , 32 ), origin = \"lower\" , cmap = \"gray\" ) ax [ 1 ] . set_xticks ([]) ax [ 1 ] . set_yticks ([]) ax [ 1 ] . set_title ( \"noisy observed data (gray image with 32 x 32 pixels)\" ) Text(0.5, 1.0, 'noisy observed data (gray image with 32 x 32 pixels)')","title":"The simulator model"},{"location":"tutorial/05_embedding_net/#defining-an-embedding_net","text":"An inference procedure applied to the output data from this simulator model determines the posterior distribution of \\(r\\) and \\(\\theta\\) given an observation of \\(x\\) , which lives in a 1024 dimensional space (32 x 32 = 1024). To avoid working directly on these high-dimensional vectors, one can use a convolutional neural network (CNN) that takes the 32x32 images as input and encodes them into 8-dimensional feature vectors. This CNN is trained along with the neural density estimator of the inference procedure and serves as an automatic summary statistics extractor. We define and instantiate the CNN as follows: class SummaryNet ( nn . Module ): def __init__ ( self ): super () . __init__ () # 2D convolutional layer self . conv1 = nn . Conv2d ( in_channels = 1 , out_channels = 6 , kernel_size = 5 , padding = 2 ) # Maxpool layer that reduces 32x32 image to 4x4 self . pool = nn . MaxPool2d ( kernel_size = 8 , stride = 8 ) # Fully connected layer taking as input the 6 flattened output arrays from the maxpooling layer self . fc = nn . Linear ( in_features = 6 * 4 * 4 , out_features = 8 ) def forward ( self , x ): x = x . view ( - 1 , 1 , 32 , 32 ) x = self . pool ( F . relu ( self . conv1 ( x ))) x = x . view ( - 1 , 6 * 4 * 4 ) x = F . relu ( self . fc ( x )) return x embedding_net = SummaryNet ()","title":"Defining an embedding_net"},{"location":"tutorial/05_embedding_net/#the-inference-procedure","text":"With the embedding_net defined and instantiated, we can follow the usual workflow of an inference procedure in sbi . The embedding_net object appears as an input argument when instantiating the neural density estimator with utils.posterior_nn . # set prior distribution for the parameters prior = utils . BoxUniform ( low = torch . tensor ([ 0.0 , 0.0 ]), high = torch . tensor ([ 1.0 , 2 * torch . pi ]) ) # make a SBI-wrapper on the simulator object for compatibility simulator_wrapper , prior = prepare_for_sbi ( simulator_model , prior ) # instantiate the neural density estimator neural_posterior = utils . posterior_nn ( model = \"maf\" , embedding_net = embedding_net , hidden_features = 10 , num_transforms = 2 ) # setup the inference procedure with the SNPE-C procedure inference = SNPE ( prior = prior , density_estimator = neural_posterior ) # run the inference procedure on one round and 10000 simulated data points theta , x = simulate_for_sbi ( simulator_wrapper , prior , num_simulations = 10000 ) Running 10000 simulations.: 0%| | 0/10000 [00:00] 1.3 Summary statistics \u00b6 We will compare two methods for defining summary statistics. One method uses three summary statistics which are function evaluations at three points in time. The other method uses a single summary statistic: the mean squared error between the observed and the simulated trace. In the second case, one then tries to obtain the posterior \\(p(\\theta | 0)\\) , i.e. the error being zero. These two methods are implemented below: \\(\\textbf{get_3_values()}\\) returns 3 function evaluations at \\(x=-0.5, x=0\\) and \\(x=0.75\\) . \\(\\textbf{get_MSE()}\\) returns the mean squared error between true and a quadratic function corresponding to a prior distributions sample. def get_3_values ( theta , seed = None ): \"\"\" Return 3 'y' values corresponding to x=-0.5,0,0.75 as summary statistic vector \"\"\" return np . array ( [ eval ( theta , - 0.5 , seed = seed ), eval ( theta , 0 , seed = seed ), eval ( theta , 0.75 , seed = seed ), ] ) . T def get_MSE ( theta , theta_o , seed = None ): \"\"\" Return the mean-squared error (MSE) i.e. Euclidean distance from the observation function \"\"\" _ , y = create_t_x ( theta_o , seed = seed ) # truth _ , y_ = create_t_x ( theta , seed = seed ) # simulations return np . mean ( np . square ( y_ - y ), axis = 0 , keepdims = True ) . T # MSE Let\u2019s try a couple of samples from our prior and see their summary statistics. Notice that these indeed change in small amounts every time you rerun it due to the noise, except if you set the seed. 1.4 Simulating data \u00b6 Let us see various plots of prior samples and their summary statistics versus the truth, i.e. our artificial observation. t , x_truth = create_t_x ( theta_o ) plt . plot ( t , x_truth , \"k\" , zorder = 1 , label = \"truth\" ) n_samples = 100 theta = prior . sample (( n_samples ,)) t , x = create_t_x ( theta . numpy ()) plt . plot ( t , x , \"grey\" , zorder = 0 ) plt . legend ()
Running the inference procedurenum_dim = 3
prior = utils.BoxUniform(low=-2 * torch.ones(num_dim), high=2 * torch.ones(num_dim))
-
def simulator(parameter_set):
return 1.0 + parameter_set + torch.randn(parameter_set.shape) * 0.1
sbi can then run inference:
-posterior = infer(simulator, prior, method="SNPE", num_simulations=1000)
+# Other methods are "SNLE" or "SNRE".
+posterior = infer(simulator, prior, method="SNPE", num_simulations=1000)
Running 1000 simulations.: 0%| | 0/1000 [00:00<?, ?it/s]
@@ -1043,32 +947,7 @@ Running the inference procedure
Next steps¶
-The single-line interface described above provides an easy entry for using sbi. However, if you are working on a larger project or need additional features, we strongly recommend using the flexible interface.
-Requirements for the simulator, prior, and observation¶
-In the interface described above, you need to provide a prior and a simulator for training. Let’s talk about what requirements they need to satisfy.
-Prior¶
-A prior is a distribution object that allows to sample parameter sets. Any class for the prior is allowed as long as it allows to call prior.sample() and prior.log_prob().
-Simulator¶
-The simulator is a Python callable that takes in a parameter set and outputs data with some (even if very small) stochasticity.
-Allowed data types and shapes for input and output:
-
-- the input parameter set and the output have to be either a
np.ndarray or a torch.Tensor.
-- the input parameter set should have either shape
(1,N) or (N), and the output must have shape (1,M) or (M).
-
-You can call simulators not written in Python as long as you wrap them in a Python function.
-Observation¶
-Once you have a trained posterior, you will want to evaluate or sample the posterior \(p(\theta|x_o)\) at certain observed values \(x_o\):
-
-- The allowable data types are either Numpy
np.ndarray or a torch torch.Tensor.
-- The shape must be either
(1,M) or just (M).
-
-Running different algorithms¶
-sbi implements three classes of algorithms that can be used to obtain the posterior distribution: SNPE, SNLE, and SNRE. You can try the different algorithms by simply swapping out the method:
-posterior = infer(simulator, prior, method="SNPE", num_simulations=1000)
-posterior = infer(simulator, prior, method="SNLE", num_simulations=1000)
-posterior = infer(simulator, prior, method="SNRE", num_simulations=1000)
-
-You can then infer, sample, evaluate, and plot the posterior as described above.
+The single-line interface described above provides an easy entry for using sbi. However, on almost any real-world problem that goes beyond a simple demonstration, we strongly recommend using the flexible interface.
@@ -1104,13 +983,13 @@ Running different algorithms
+
diff --git a/tutorial/09_sensitivity_analysis/index.html b/tutorial/09_sensitivity_analysis/index.html
index b0e691433..b7a80d3cf 100644
--- a/tutorial/09_sensitivity_analysis/index.html
+++ b/tutorial/09_sensitivity_analysis/index.html
@@ -311,20 +311,6 @@
-
-
- Amortized inference
-
-
-
-
-
-
-
-
-
-
-
Flexible interface
@@ -340,8 +326,8 @@
-
- Sampler interface
+
+ Amortized inference
@@ -415,8 +401,8 @@
-
- Using Variational Inference for Building Posteriors
+
+ Sampling algorithms in sbi
@@ -457,8 +443,8 @@
-
- Handling invalid simulations
+
+ SBI with trial-based data
@@ -471,8 +457,8 @@
-
- Crafting summary statistics
+
+ Handling invalid simulations
@@ -485,8 +471,8 @@
-
- SBI with trial-based data
+
+ Crafting summary statistics
diff --git a/tutorial/10_crafting_summary_statistics/index.html b/tutorial/10_crafting_summary_statistics/index.html
index d842e288a..f69c3a25c 100644
--- a/tutorial/10_crafting_summary_statistics/index.html
+++ b/tutorial/10_crafting_summary_statistics/index.html
@@ -311,20 +311,6 @@
-
-
- Amortized inference
-
-
-
-
-
-
-
-
-
-
-
Flexible interface
@@ -340,8 +326,8 @@
-
- Sampler interface
+
+ Amortized inference
@@ -417,8 +403,8 @@
-
- Using Variational Inference for Building Posteriors
+
+ Sampling algorithms in sbi
@@ -458,6 +444,20 @@
+
+
+ SBI with trial-based data
+
+
+
+
+
+
+
+
+
+
+
Handling invalid simulations
@@ -491,20 +491,6 @@
-
-
-
-
-
-
-
- SBI with trial-based data
-
-
-
-
-
-
@@ -1167,13 +1153,13 @@ 1.7 Explicit recommendations
+
@@ -933,7 +844,7 @@
- The sampler interface¶
+ Sampling algorithms in sbi¶
Note: this tutorial requires that the user is already familiar with the flexible interface.
sbi implements three methods: SNPE, SNLE, and SNRE. When using SNPE, the trained neural network directly approximates the posterior. Thus, sampling from the posterior can be done by sampling from the trained neural network. The neural networks trained in SNLE and SNRE approximate the likelihood(-ratio). Thus, in order to draw samples from the posterior, one has to perform additional sampling steps, e.g. Markov-chain Monte-Carlo (MCMC). In sbi, the implemented samplers are:
@@ -947,12 +858,44 @@ The sampler interface, mcmc_method="slice_np", and mcmc_parameters={}. However, for full flexibility in customizing the sampler, we recommend using the sampler interface. This interface is described here. Further details can be found here.
-Main syntax for SNLE¶
+
Below, we will demonstrate how these samplers can be used in sbi. First, we train the neural network as always:
+import torch
+from sbi.inference import SNLE
+
+# dummy Gaussian simulator for demonstration
+num_dim = 2
+prior = torch.distributions.MultivariateNormal(torch.zeros(num_dim), torch.eye(num_dim))
+theta = prior.sample((1000,))
+x = theta + torch.randn((1000, num_dim))
+x_o = torch.randn((1, num_dim))
+
+inference = SNLE(prior=prior, show_progress_bars=False)
+likelihood_estimator = inference.append_simulations(theta, x).train()
+
+And then we pass the options for which sampling method to use to the build_posterior() method:
+# Sampling with MCMC
+sampling_algorithm = "mcmc"
+mcmc_method = "slice_np" # or nuts, or hmc
+posterior = inference.build_posterior(sample_with=sampling_algorithm, mcmc_method=mcmc_method)
+
+# Sampling with variational inference
+sampling_algorithm = "vi"
+vi_method = "rKL" # or fKL
+posterior = inference.build_posterior(sample_with=sampling_algorithm, vi_method=vi_method)
+# Unlike other methods, vi needs a training step for every observation.
+posterior = posterior.set_default_x(x_o).train()
+
+# Sampling with rejection sampling
+sampling_algorithm = "rejection"
+posterior = inference.build_posterior(sample_with=sampling_algorithm)
+
+More flexibility in adjusting the sampler¶
+With the above syntax, you can easily try out different sampling algorithms. However, in many cases, you might want to customize your sampler. Below, we demonstrate how you can change hyperparameters of the samplers (e.g. number of warm-up steps of MCMC) or how you can write your own sampler from scratch.
+Main syntax (for SNLE and SNRE)¶
+As above, we begin by training the neural network as always:
import torch
from sbi.inference import SNLE
-from sbi.inference import likelihood_estimator_based_potential, MCMCPosterior
# dummy Gaussian simulator for demonstration
num_dim = 2
@@ -963,17 +906,31 @@ Main syntax for SNLEinference = SNLE(show_progress_bars=False)
likelihood_estimator = inference.append_simulations(theta, x).train()
+
+ Neural network successfully converged after 52 epochs.
+
+
+Then, for full flexibility on using the sampler, we do not use the .build_posterior() method, but instead we explicitly define the potential function and the sampling algorithm (see below for explanation):
+from sbi.inference import likelihood_estimator_based_potential, MCMCPosterior
potential_fn, parameter_transform = likelihood_estimator_based_potential(
likelihood_estimator, prior, x_o
)
posterior = MCMCPosterior(
- potential_fn, proposal=prior, theta_transform=parameter_transform
+ potential_fn, proposal=prior, theta_transform=parameter_transform, warmup_steps=10
)
- Neural network successfully converged after 52 epochs.
-
+If you want to use variational inference or rejection sampling, you have to replace the last line with VIPosterior or RejectionPosterior:
+# For VI, we have to train.
+posterior = VIPosterior(
+ potential_fn, proposal=prior, theta_transform=parameter_transform
+).train()
+posterior = RejectionPosterior(
+ potential_fn, proposal=prior, theta_transform=parameter_transform
+)
+
+At this point, you could also plug the potential_fn into any sampler of your choice and not rely on any of the in-built sbi-samplers.
Further explanation¶
The first lines are the same as for the flexible interface:
-
+
diff --git a/tutorial/13_diagnostics_simulation_based_calibration/index.html b/tutorial/13_diagnostics_simulation_based_calibration/index.html
index ade247774..8bba7e91b 100644
--- a/tutorial/13_diagnostics_simulation_based_calibration/index.html
+++ b/tutorial/13_diagnostics_simulation_based_calibration/index.html
@@ -311,20 +311,6 @@
-
num_dim = 3
prior = utils.BoxUniform(low=-2 * torch.ones(num_dim), high=2 * torch.ones(num_dim))
-
def simulator(parameter_set):
return 1.0 + parameter_set + torch.randn(parameter_set.shape) * 0.1
sbi can then run inference:posterior = infer(simulator, prior, method="SNPE", num_simulations=1000)
+# Other methods are "SNLE" or "SNRE".
+posterior = infer(simulator, prior, method="SNPE", num_simulations=1000)
Running 1000 simulations.: 0%| | 0/1000 [00:00<?, ?it/s]
@@ -1043,32 +947,7 @@ Running the inference procedure
Next steps¶
-The single-line interface described above provides an easy entry for using sbi. However, if you are working on a larger project or need additional features, we strongly recommend using the flexible interface.
-Requirements for the simulator, prior, and observation¶
-In the interface described above, you need to provide a prior and a simulator for training. Let’s talk about what requirements they need to satisfy.
-Prior¶
-A prior is a distribution object that allows to sample parameter sets. Any class for the prior is allowed as long as it allows to call prior.sample() and prior.log_prob().
-Simulator¶
-The simulator is a Python callable that takes in a parameter set and outputs data with some (even if very small) stochasticity.
-Allowed data types and shapes for input and output:
-
-- the input parameter set and the output have to be either a
np.ndarray or a torch.Tensor.
-- the input parameter set should have either shape
(1,N) or (N), and the output must have shape (1,M) or (M).
-
-You can call simulators not written in Python as long as you wrap them in a Python function.
-Observation¶
-Once you have a trained posterior, you will want to evaluate or sample the posterior \(p(\theta|x_o)\) at certain observed values \(x_o\):
-
-- The allowable data types are either Numpy
np.ndarray or a torch torch.Tensor.
-- The shape must be either
(1,M) or just (M).
-
-Running different algorithms¶
-sbi implements three classes of algorithms that can be used to obtain the posterior distribution: SNPE, SNLE, and SNRE. You can try the different algorithms by simply swapping out the method:
-posterior = infer(simulator, prior, method="SNPE", num_simulations=1000)
-posterior = infer(simulator, prior, method="SNLE", num_simulations=1000)
-posterior = infer(simulator, prior, method="SNRE", num_simulations=1000)
-
-You can then infer, sample, evaluate, and plot the posterior as described above.
+The single-line interface described above provides an easy entry for using sbi. However, on almost any real-world problem that goes beyond a simple demonstration, we strongly recommend using the flexible interface.
@@ -1104,13 +983,13 @@ Running different algorithms
+
diff --git a/tutorial/09_sensitivity_analysis/index.html b/tutorial/09_sensitivity_analysis/index.html
index b0e691433..b7a80d3cf 100644
--- a/tutorial/09_sensitivity_analysis/index.html
+++ b/tutorial/09_sensitivity_analysis/index.html
@@ -311,20 +311,6 @@
-
-
- Amortized inference
-
-
-
-
-
-
-
-
-
-
-
Flexible interface
@@ -340,8 +326,8 @@
-
- Sampler interface
+
+ Amortized inference
@@ -415,8 +401,8 @@
-
- Using Variational Inference for Building Posteriors
+
+ Sampling algorithms in sbi
@@ -457,8 +443,8 @@
-
- Handling invalid simulations
+
+ SBI with trial-based data
@@ -471,8 +457,8 @@
-
- Crafting summary statistics
+
+ Handling invalid simulations
@@ -485,8 +471,8 @@
-
- SBI with trial-based data
+
+ Crafting summary statistics
diff --git a/tutorial/10_crafting_summary_statistics/index.html b/tutorial/10_crafting_summary_statistics/index.html
index d842e288a..f69c3a25c 100644
--- a/tutorial/10_crafting_summary_statistics/index.html
+++ b/tutorial/10_crafting_summary_statistics/index.html
@@ -311,20 +311,6 @@
-
-
- Amortized inference
-
-
-
-
-
-
-
-
-
-
-
Flexible interface
@@ -340,8 +326,8 @@
-
- Sampler interface
+
+ Amortized inference
@@ -417,8 +403,8 @@
-
- Using Variational Inference for Building Posteriors
+
+ Sampling algorithms in sbi
@@ -458,6 +444,20 @@
+
+
+ SBI with trial-based data
+
+
+
+
+
+
+
+
+
+
+
Handling invalid simulations
@@ -491,20 +491,6 @@
-
-
-
-
-
-
-
- SBI with trial-based data
-
-
-
-
-
-
@@ -1167,13 +1153,13 @@ 1.7 Explicit recommendations
+
@@ -933,7 +844,7 @@
- The sampler interface¶
+ Sampling algorithms in sbi¶
Note: this tutorial requires that the user is already familiar with the flexible interface.
sbi implements three methods: SNPE, SNLE, and SNRE. When using SNPE, the trained neural network directly approximates the posterior. Thus, sampling from the posterior can be done by sampling from the trained neural network. The neural networks trained in SNLE and SNRE approximate the likelihood(-ratio). Thus, in order to draw samples from the posterior, one has to perform additional sampling steps, e.g. Markov-chain Monte-Carlo (MCMC). In sbi, the implemented samplers are:
@@ -947,12 +858,44 @@ The sampler interface, mcmc_method="slice_np", and mcmc_parameters={}. However, for full flexibility in customizing the sampler, we recommend using the sampler interface. This interface is described here. Further details can be found here.
-Main syntax for SNLE¶
+
Below, we will demonstrate how these samplers can be used in sbi. First, we train the neural network as always:
+import torch
+from sbi.inference import SNLE
+
+# dummy Gaussian simulator for demonstration
+num_dim = 2
+prior = torch.distributions.MultivariateNormal(torch.zeros(num_dim), torch.eye(num_dim))
+theta = prior.sample((1000,))
+x = theta + torch.randn((1000, num_dim))
+x_o = torch.randn((1, num_dim))
+
+inference = SNLE(prior=prior, show_progress_bars=False)
+likelihood_estimator = inference.append_simulations(theta, x).train()
+
+And then we pass the options for which sampling method to use to the build_posterior() method:
+# Sampling with MCMC
+sampling_algorithm = "mcmc"
+mcmc_method = "slice_np" # or nuts, or hmc
+posterior = inference.build_posterior(sample_with=sampling_algorithm, mcmc_method=mcmc_method)
+
+# Sampling with variational inference
+sampling_algorithm = "vi"
+vi_method = "rKL" # or fKL
+posterior = inference.build_posterior(sample_with=sampling_algorithm, vi_method=vi_method)
+# Unlike other methods, vi needs a training step for every observation.
+posterior = posterior.set_default_x(x_o).train()
+
+# Sampling with rejection sampling
+sampling_algorithm = "rejection"
+posterior = inference.build_posterior(sample_with=sampling_algorithm)
+
+More flexibility in adjusting the sampler¶
+With the above syntax, you can easily try out different sampling algorithms. However, in many cases, you might want to customize your sampler. Below, we demonstrate how you can change hyperparameters of the samplers (e.g. number of warm-up steps of MCMC) or how you can write your own sampler from scratch.
+Main syntax (for SNLE and SNRE)¶
+As above, we begin by training the neural network as always:
import torch
from sbi.inference import SNLE
-from sbi.inference import likelihood_estimator_based_potential, MCMCPosterior
# dummy Gaussian simulator for demonstration
num_dim = 2
@@ -963,17 +906,31 @@ Main syntax for SNLEinference = SNLE(show_progress_bars=False)
likelihood_estimator = inference.append_simulations(theta, x).train()
+
+ Neural network successfully converged after 52 epochs.
+
+
+Then, for full flexibility on using the sampler, we do not use the .build_posterior() method, but instead we explicitly define the potential function and the sampling algorithm (see below for explanation):
+from sbi.inference import likelihood_estimator_based_potential, MCMCPosterior
potential_fn, parameter_transform = likelihood_estimator_based_potential(
likelihood_estimator, prior, x_o
)
posterior = MCMCPosterior(
- potential_fn, proposal=prior, theta_transform=parameter_transform
+ potential_fn, proposal=prior, theta_transform=parameter_transform, warmup_steps=10
)
- Neural network successfully converged after 52 epochs.
-
+If you want to use variational inference or rejection sampling, you have to replace the last line with VIPosterior or RejectionPosterior:
+# For VI, we have to train.
+posterior = VIPosterior(
+ potential_fn, proposal=prior, theta_transform=parameter_transform
+).train()
+posterior = RejectionPosterior(
+ potential_fn, proposal=prior, theta_transform=parameter_transform
+)
+
+At this point, you could also plug the potential_fn into any sampler of your choice and not rely on any of the in-built sbi-samplers.
Further explanation¶
The first lines are the same as for the flexible interface:
-
+
SBI with iid data and permutation-invariant embeddings¶
There are scenarios in which we observe multiple data points per experiment and we can assume that they are independent and identically distributed (iid, i.e., they are assumed to have the same underlying model parameters). -For example, in a decision-making experiments, the experiment is often repeated in trials with the same experimental settings and conditions. The corresponding set of trials is then assumed to be “iid”. +For example, in decision-making experiments, the experiment is often repeated in trials with the same experimental settings and conditions. The corresponding set of trials is then assumed to be “iid” given a single parameter set. In such a scenario, we may want to obtain the posterior given a set of observation \(p(\theta | X=\{x_i\}_i^N)\).
Amortization of neural network training: iid-inference with NLE / NRE¶
-For some SBI variants the iid assumption can be exploited: when using a likelihood-based SBI method (SNLE, SNRE) one can train the density or ratio estimator on single-trial data, and then perform inference with MCMC. Crucially, because the data is iid and the estimator is trained on single-trial data, one can repeat the inference with a different x_o (a different set of trials, or different number of trials) without having to retrain the density estimator. One can interpet this as amortization of the SBI training: we can obtain a neural likelihood, or likelihood-ratio estimate for new x_os without retraining, but we still have to run MCMC or VI to do inference.
In addition, one can not only change the number of trials of a new x_o, but also the entire inference setting.
-For example, one can apply hierarchical inference scenarios with changing hierarchical denpendencies between the model parameters–all without having to retrain the density estimator because that is based on estimating single-trail likelihoods.
For some SBI variants the iid assumption can be exploited: when using a likelihood-based SBI method (SNLE, SNRE) one can train the density or ratio estimator on single-trial data, and then perform inference with MCMC or variational inference (VI). Crucially, because the data is iid and the estimator is trained on single-trial data, one can repeat the inference with a different x_o (a different set of trials, or different number of trials) without having to retrain the density estimator. One can interpet this as amortization of the SBI training: we can obtain a neural likelihood, or likelihood-ratio estimate for new x_os without retraining, but we still have to run MCMC or VI to do inference.
In addition, one cannot only change the number of trials of a new x_o, but also the entire inference setting.
+For example, one can apply hierarchical inference with changing hierarchical denpendencies between the model parameters–all without having to retrain the density estimator because it estimates single-trail likelihoods.
Full amortization: iid-inference with NPE and permutation-invariant embedding nets¶
-When performing neural posterior estimation (SNPE) we cannot exploit the iid assumption directly because we are learning a density estimator in theta.
+
When performing neural posterior estimation (SNPE) we cannot exploit the iid assumption directly.
Thus, the underlying neural network takes x as input and predicts the parameters of the density estimator.
-As a consequence, if x is a set of iid observations \(X=\{x_i\}_i^N\) then the neural network has to be invariant to permutations of this set, i.e., it has to be permutation invariant.
-Overall, this means that we can use SNPE for inference with iid data, however, we need to provide a corresponding embedding network that handles the iid-data and is permutation invariant.
-This will likely require some hyperparameter tuning and more training data for the inference to work accurately. But once we have this, the inference is fully amortized, i.e., we can get new posterior samples basically instantly without retraining and without running MCMC or VI.
Let us first have a look how trial-based inference works in SBI before we discuss models with “mixed data types”.
x is a set of iid observations \(X=\{x_i\}_i^N\) then the neural network has to be invariant to permutations of this set, i.e., it has to be permutation invariant. In addition, the neural network has to be able to consume a varying number of iid datapoints in order to be amortized over the number of trials.
+Therefore, in order to use SNPE for inference on iid data, we need to provide a corresponding embedding network that handles the iid-data.
+This will likely require some hyperparameter tuning and more training data for inference to work accurately. But once we have this, inference is fully amortized, i.e., we can get new posterior samples almost instantly without retraining and without running MCMC or VI.
SBI with trial-based data¶
-For illustration we use a simple linear Gaussian simulator, as in previous tutorials. The simulator takes a single parameter (vector), the mean of the Gaussian, and its variance is set to one.
-We define a Gaussian prior over the mean and perform inference.
-The observed data is again a from a Gaussian with some fixed “ground-truth” parameter \(\theta_o\).
-Crucially, the observed data x_o can consist of multiple samples given the same ground-truth parameters and these samples are then iid:
For illustration, we use a simple linear Gaussian simulator, as in previous tutorials. The simulator takes a single parameter (vector) which is the mean of a Gaussian. The simulator then adds noise with a fixed variance (set to one). +We define a Gaussian prior over the mean and perform inference.
+The observed data is also sampled from a Gaussian with some fixed “ground-truth” parameter \(\theta_o\).
+Crucially, the observed data x_o can consist of multiple samples given the same ground-truth parameters and these samples are iid given \(\theta\):
\[
\theta \sim \mathcal{N}(\mu_0,\; \Sigma_0) \\
x | \theta \sim \mathcal{N}(\theta,\; \Sigma=I) \\
\mathbf{x_o} = \{x_o^i\}_{i=1}^N \sim \mathcal{N}(\theta_o,\; \Sigma=I)
\]
-For this toy problem the ground-truth posterior is well defined, it is again a Gaussian, centered on the mean of \(\mathbf{x_o}\) and with variance scaled by the number of trials \(N\), i.e., the more trials we observe, the more information about the underlying \(\theta_o\) we have and the more concentrated the posteriors becomes.
+For this toy problem, the ground-truth posterior is well defined, it is again a Gaussian, centered on the mean of \(\mathbf{x_o}\) and with variance scaled by the number of trials \(N\), i.e., the more trials we observe, the more information about the underlying \(\theta_o\) we have and the more concentrated the posteriors becomes.
We will illustrate this below:
import torch
import matplotlib.pyplot as plt
@@ -1156,8 +1141,7 @@
Indeed, with increasing number of trials the posterior density concentrates around the true underlying parameter.
IID inference with NLE¶
-(S)NLE can easily perform inference given multiple IID x because it is based on learning the likelihood. Once the likelihood is learned on single trials, i.e., a neural network that given a single observation and a parameter predicts the likelihood of that observation given the parameter, one can perform MCMC to obtain posterior samples.
-MCMC relies on evaluating ratios of likelihoods of candidate parameters to either accept or reject them to be posterior samples. When inferring the posterior given multiple IID observation, these likelihoods are just the joint likelihoods of each IID observation given the current parameter candidate. Thus, given a neural likelihood from SNLE, we can calculate these joint likelihoods and perform MCMC given IID data, we just have to multiply together (or add in log-space) the individual trial-likelihoods (sbi takes care of that).
+(S)NLE and (S)NRE can perform inference given multiple IID obserations by using only single-trial training data (i.e., for training, we run the simulator only once per parameter set). Once the likelihood is learned on single trials (i.e., a neural network that predicts the likelihood of a single observation given a parameter set), one can sample the posterior for any number of trials. This works because, given a single-trial neural likelihood from (S)NLE or (S)NRE, we can calculate the joint likelihoods of all trials by multiplying them together (or adding them in log-space). The joint likelihood can then be plugged into MCMC or VI. sbi takes care of all of these steps, so you do not have to implement anything yourself:
# Train SNLE.
inferer = SNLE(prior, show_progress_bars=True, density_estimator="mdn")
theta, x = simulate_for_sbi(simulator, prior, 10000, simulation_batch_size=1000)
@@ -1265,10 +1249,11 @@ IID inference with NLE
IID inference with NPE using permutation-invariant embedding nets¶
-For NPE we need to define an embedding net that handles the set-like structure of iid-data, i.e., that it permutation invariant and can handle different number of trials.
+For NPE we need to define an embedding net that handles the set-like structure of iid-data, i.e., that it permutation invariant and can handle different number of trials.
We implemented several embedding net classes that allow to construct such a permutation- and number-of-trials invariant embedding net.
To become permutation invariant, the neural net first learns embeddings for single trials and then performs a permutation invariant operation on those embeddings, e.g., by taking the sum or the mean (Chen et al. 2018, Radev et al. 2021).
-To become invariant w.r.t. the number-of-trials, we train the net with varying number of trials for each parameter setting. As it is difficult to handle tensors of varying lengths in the SBI training loop, we construct a training data set in which “unobserved” trials are mask by NaNs (and ignore the resulting SBI warning about NaNs in the training data).
+To become invariant w.r.t. the number-of-trials, we train the net with varying number of trials for each parameter setting. This means that, unlike for (S)NLE and (S)NRE, (S)NPE requires to run the simulator multiple times for individual parameter sets to generate the training data.
+In order to implement this in sbi, “unobserved” trials in the training dataset have to be masked by NaNs (and ignore the resulting SBI warning about NaNs in the training data).
Construct training data set.¶
Using Variational Inference for Building Posteriors¶
-In the previous tutorial, we saw how to build the posterior and how to specialize on one specific observation x_o. If one uses SNPE, then the posterior can be sampled from directly, yet this comes at the expense of necessary correction terms during training, since the samples are obtained from the “wrong” prior for num_rounds > 1. For SNLE or SNRE, MCMC sampling is required, which is computationally expensive. With SNVI (sequential neural variational inference), it is possible to directly sample from the posterior without any corrections during training or without expensive MCMC for sampling. This is possible by learning the posterior with variational inference techniques. For this, an additional network (one for the likelihood or likelihood-to-evidence-ratio) must be trained first.
If one uses SNPE, then the posterior can be sampled from directly (without MCMC). Contrary to that, for SNLE or SNRE, MCMC sampling is required, which is computationally expensive. With SNVI (sequential neural variational inference), it is possible to directly sample from the posterior without any corrections during training or without expensive MCMC for sampling. This is possible by learning the posterior with variational inference techniques. For this, an additional network (one for the likelihood or likelihood-to-evidence-ratio) must be trained first.
Main syntax¶
inference = SNLE(prior)
@@ -1047,41 +980,6 @@ Linear Gaussian example
-
-