diff --git a/.travis.yml b/.travis.yml index 92de1003..259ef3bb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,6 +11,7 @@ addons: - gfortran - libboost-random-dev - python3.7-dev + - python3.9-dev - python3-numpy - swig - libmpich-dev @@ -31,6 +32,10 @@ jobs: python: "3.8" env: - UNIT_TEST=true + - name: "Python 3.9" + python: "3.9" + env: + - UNIT_TEST=true # Test coverage and without Pytorch for a single version - name: "Coverage" python: "3.8" @@ -61,6 +66,7 @@ jobs: install: - sudo apt-get install -y r-base # install R for testing example +- pip install cython - pip install rpy2 # install the rpy2 library for testing example with R - pip install -r requirements.txt - pip install -r requirements/backend-spark.txt diff --git a/Makefile b/Makefile index f92a5f65..74e85769 100644 --- a/Makefile +++ b/Makefile @@ -4,10 +4,10 @@ MAKEDIRS=$(shell find examples -name Makefile -exec dirname {} \;) whl_file = abcpy-${VERSION}-py3-none-any.whl .DEFAULT: help -.PHONY: help clean doc doctest exampletest package test uninstall unittest unittest_mpi install reinstall $(MAKEDIRS) +.PHONY: help clean doc doctest exampletest exampletest_mpi package test uninstall unittest unittest_mpi install reinstall $(MAKEDIRS) help: - @echo Targets are: clean, doc, doctest, exampletest, package, uninstall, unittest, unittest_mpi , test + @echo Targets are: clean, doc, doctest, exampletest, exampletest_mpi, package, uninstall, unittest, unittest_mpi, test clean: find . -name "*.pyc" -type f -delete @@ -26,6 +26,11 @@ test: unittest unittest_mpi exampletest exampletest_mpi doctest unittest: @echo "Running standard unit tests.." python3 -m unittest discover -s tests -v -p "*_tests.py" || (echo "Error in standard unit tests."; exit 1) + @# remove temporary files created during testing + @if test -f net.pth; then rm net.pth; fi + @if test -f scaler.pkl; then rm scaler.pkl; fi + @if test -f tmp.jnl; then rm tmp.jnl; fi + @if test -f journal_tests_testfile.pkl; then rm journal_tests_testfile.pkl; fi unittest_mpi: @echo "Running MPI backend unit tests.." diff --git a/README.md b/README.md index 7734297e..fe375d96 100644 --- a/README.md +++ b/README.md @@ -2,19 +2,50 @@ ABCpy is a scientific library written in Python for Bayesian uncertainty quantification in absence of likelihood function, which parallelizes existing approximate Bayesian computation (ABC) -algorithms and other likelihood-free inference schemes. It presently includes: - -* RejectionABC -* PMCABC (Population Monte Carlo ABC) -* SMCABC (Sequential Monte Carlo ABC) -* RSMCABC (Replenishment SMC-ABC) -* APMCABC (Adaptive Population Monte Carlo ABC) -* SABC (Simulated Annealing ABC) -* ABCsubsim (ABC using subset simulation) -* PMC (Population Monte Carlo) using approximations of likelihood functions -* Random Forest Model Selection Scheme -* Semi-automatic summary selection (with Neural networks) -* summary selection using distance learning (with Neural networks) +algorithms and other likelihood-free inference schemes. + +# Content + +ABCpy presently includes the following **ABC algorithms**: + +* [RejectionABC](https://www.genetics.org/content/145/2/505) +* [PMCABC (Population Monte Carlo ABC)](https://www.annualreviews.org/doi/abs/10.1146/annurev-ecolsys-102209-144621) +* [SMCABC (Sequential Monte Carlo ABC)](https://link.springer.com/article/10.1007/s11222-011-9271-y) +* [RSMCABC (Replenishment SMC-ABC)](https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1541-0420.2010.01410.x) +* [APMCABC (Adaptive Population Monte Carlo ABC)](https://link.springer.com/article/10.1007/s00180-013-0428-3) +* [SABC (Simulated Annealing ABC)](https://link.springer.com/article/10.1007/s11222-014-9507-8) +* [ABCsubsim (ABC using subset simulation)](https://epubs.siam.org/doi/10.1137/130932831) + +The above can be used with the following **distances**: + +* Euclidean Distance +* [Logistic Regression and Penalised Logistic Regression (classification accuracy)](https://link.springer.com/article/10.1007/s11222-017-9738-6) +* Divergences between datasets: + * [Wasserstein Distance](https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/rssb.12312) + * [Sliced Wasserstein Distance](https://ieeexplore.ieee.org/abstract/document/9054735) + * [Gamma Divergence](http://proceedings.mlr.press/v130/fujisawa21a/fujisawa21a.pdf) + * [Kullback Liebler Divergence](http://proceedings.mlr.press/v84/jiang18a/jiang18a.pdf) + * [Maximum Mean Discrepancy](http://proceedings.mlr.press/v51/park16.pdf) + * [Energy Distance](https://arxiv.org/abs/1905.05884) + * [Squared Hellinger Distance](https://arxiv.org/pdf/2006.14126.pdf) + +Moreover, we provide the following methods for directly **approximating the likelihood functions**: +* [Bayesian Synthetic Likelihood](https://www.tandfonline.com/doi/abs/10.1080/10618600.2017.1302882?journalCode=ucgs20) +* [Semiparametric Bayesian Synthetic Likelihood](https://link.springer.com/article/10.1007/s11222-019-09904-x) +* [Penalised Logistic Regression for Ratio Estimation](https://projecteuclid.org/journals/bayesian-analysis/advance-publication/Likelihood-Free-Inference-by-Ratio-Estimation/10.1214/20-BA1238.full) + +The above likelihood approximation methods can be used with the following samplers: + +* [PMC (Population Monte Carlo)](https://www.tandfonline.com/doi/abs/10.1198/106186004X12803) +* Metropolis-Hastings MCMC (Markov Chain Monte Carlo) + +Additional **features** are: +* plotting utilities for the obtained posterior +* several methods for summary selection: + * [Semi-automatic summary selection (with Neural networks)](http://proceedings.mlr.press/v97/wiqvist19a/wiqvist19a.pdf) + * [summary selection using distance learning (with Neural networks)](https://link.springer.com/article/10.1007/s13571-019-00208-8) +* [Random Forest Model Selection Scheme](https://academic.oup.com/bioinformatics/article/32/6/859/1744513) + ABCpy addresses the needs of domain scientists and data scientists by providing @@ -27,9 +58,9 @@ scientists by providing # Documentation For more information, check out the -* [Documentation](http://abcpy.readthedocs.io/en/v0.6.1) -* [Examples](https://github.com/eth-cscs/abcpy/tree/v0.6.1/examples) directory and -* [Reference](http://abcpy.readthedocs.io/en/v0.6.1/abcpy.html) +* [Documentation](http://abcpy.readthedocs.io/en/v0.6.2) +* [Examples](https://github.com/eth-cscs/abcpy/tree/v0.6.2/examples) directory and +* [Reference](http://abcpy.readthedocs.io/en/v0.6.2/abcpy.html) Further, we provide a @@ -93,33 +124,25 @@ ABCpy for your publication, we would appreciate a citation. You can use Publications in which ABCpy was applied: +* L. Pacchiardi, R. Dutta. "Generalized Bayesian Likelihood-Free Inference Using Scoring Rules Estimators", 2021, arXiv:2104.03889. + * L. Pacchiardi, R. Dutta. "Score Matched Conditional Exponential Families for Likelihood-Free Inference", 2020, arXiv:2012.10903. -* R. Dutta, K. Zouaoui-Boudjeltia, C. Kotsalos, A. Rousseau, D. Ribeiro de Sousa, J. M. Desmet, -A. Van Meerhaeghe, A. Mira, and B. Chopard. "Interpretable pathological test for Cardio-vascular -disease: Approximate Bayesian computation with distance learning.", 2020, arXiv:2010.06465. +* R. Dutta, K. Zouaoui-Boudjeltia, C. Kotsalos, A. Rousseau, D. Ribeiro de Sousa, J. M. Desmet, A. Van Meerhaeghe, A. Mira, and B. Chopard. "Interpretable pathological test for Cardio-vascular disease: Approximate Bayesian computation with distance learning.", 2020, arXiv:2010.06465. -* R. Dutta, S. Gomes, D. Kalise, L. Pacchiardi. "Using mobility data in the design of optimal -lockdown strategies for the COVID-19 pandemic in England.", 2020, arXiv:2006.16059. +* R. Dutta, S. Gomes, D. Kalise, L. Pacchiardi. "Using mobility data in the design of optimal lockdown strategies for the COVID-19 pandemic in England.", 2020, arXiv:2006.16059. -* L. Pacchiardi, P. Künzli, M. Schöngens, B. Chopard, R. Dutta, "Distance-Learning for -Approximate Bayesian Computation to Model a Volcanic Eruption", 2020, Sankhya B, ISSN 0976-8394, - [DOI: 10.1007/s13571-019-00208-8](https://doi.org/10.1007/s13571-019-00208-8). +* L. Pacchiardi, P. Künzli, M. Schöngens, B. Chopard, R. Dutta, "Distance-Learning for Approximate Bayesian Computation to Model a Volcanic Eruption", 2020, Sankhya B, 1-30. -* R. Dutta, J. P. Onnela, A. Mira, "Bayesian Inference of Spreading Processes - on Networks", 2018, Proc. R. Soc. A, 474(2215), 20180129. +* R. Dutta, J. P. Onnela, A. Mira, "Bayesian Inference of Spreading Processes on Networks", 2018, Proceedings of Royal Society A, 474(2215), 20180129. -* R. Dutta, Z. Faidon Brotzakis and A. Mira, "Bayesian Calibration of - Force-fields from Experimental Data: TIP4P Water", 2018, Journal of Chemical Physics 149, 154110. +* R. Dutta, Z. Faidon Brotzakis and A. Mira, "Bayesian Calibration of Force-fields from Experimental Data: TIP4P Water", 2018, Journal of Chemical Physics 149, 154110. -* R. Dutta, B. Chopard, J. Lätt, F. Dubois, K. Zouaoui Boudjeltia and A. Mira, - "Parameter Estimation of Platelets Deposition: Approximate Bayesian - Computation with High Performance Computing", 2018, Frontiers in physiology, 9. +* R. Dutta, B. Chopard, J. Lätt, F. Dubois, K. Zouaoui Boudjeltia and A. Mira, "Parameter Estimation of Platelets Deposition: Approximate Bayesian Computation with High Performance Computing", 2018, Frontiers in physiology, 9. -* A. Ebert, R. Dutta, P. Wu, K. Mengersen and A. Mira, "Likelihood-Free - Parameter Estimation for Dynamic Queueing Networks", 2018, arXiv:1804.02526. +* A. Ebert, R. Dutta, P. Wu, K. Mengersen and A. Mira, "Likelihood-Free Parameter Estimation for Dynamic Queueing Networks", 2018, arXiv:1804.02526, To apear in Jouranl of Royal Statistical Scoiety: Series C. -* R. Dutta, M. Schoengens, L. Pacchiardi, A. Ummadisingu, N. Widerman, J. P. Onnela, A. Mira, "ABCpy: A High-Performance Computing Perspective to Approximate Bayesian Computation", 2020, arXiv:1711.04694. +* R. Dutta, M. Schoengens, L. Pacchiardi, A. Ummadisingu, N. Widerman, P. Künzli, J. P. Onnela, A. Mira, "ABCpy: A High-Performance Computing Perspective to Approximate Bayesian Computation", 2018, arXiv:1711.04694, To appear in Journal of Statistical Software. ## License ABCpy is published under the BSD 3-clause license, see [here](LICENSE). diff --git a/VERSION b/VERSION index ee6cdce3..b1d7abc0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.6.1 +0.6.2 \ No newline at end of file diff --git a/abcpy/NN_utilities/algorithms.py b/abcpy/NN_utilities/algorithms.py index dd150772..4a17b0e5 100644 --- a/abcpy/NN_utilities/algorithms.py +++ b/abcpy/NN_utilities/algorithms.py @@ -19,7 +19,7 @@ def contrastive_training(samples, similarity_set, embedding_net, cuda, batch_siz samples_val=None, similarity_set_val=None, early_stopping=False, epochs_early_stopping_interval=1, start_epoch_early_stopping=10, positive_weight=None, load_all_data_GPU=False, margin=1., - lr=None, optimizer=None, scheduler=None, start_epoch_training=0, + lr=None, optimizer=None, scheduler=None, start_epoch_training=0, use_tqdm=True, optimizer_kwargs={}, scheduler_kwargs={}, loader_kwargs={}): """ Implements the algorithm for the contrastive distance learning training of a neural network; need to be provided with a set of samples and the corresponding similarity matrix""" @@ -84,7 +84,8 @@ def contrastive_training(samples, similarity_set, embedding_net, cuda, batch_siz fit(pairs_train_loader, model_contrastive, loss_fn, optimizer, scheduler, n_epochs, cuda, val_loader=pairs_train_loader_val, early_stopping=early_stopping, start_epoch_early_stopping=start_epoch_early_stopping, - epochs_early_stopping_interval=epochs_early_stopping_interval, start_epoch_training=start_epoch_training) + epochs_early_stopping_interval=epochs_early_stopping_interval, start_epoch_training=start_epoch_training, + use_tqdm=use_tqdm) return embedding_net @@ -93,7 +94,7 @@ def triplet_training(samples, similarity_set, embedding_net, cuda, batch_size=16 samples_val=None, similarity_set_val=None, early_stopping=False, epochs_early_stopping_interval=1, start_epoch_early_stopping=10, load_all_data_GPU=False, margin=1., lr=None, optimizer=None, scheduler=None, - start_epoch_training=0, + start_epoch_training=0, use_tqdm=True, optimizer_kwargs={}, scheduler_kwargs={}, loader_kwargs={}): """ Implements the algorithm for the triplet distance learning training of a neural network; need to be provided with a set of samples and the corresponding similarity matrix""" @@ -157,7 +158,7 @@ def triplet_training(samples, similarity_set, embedding_net, cuda, batch_size=16 fit(triplets_train_loader, model_triplet, loss_fn, optimizer, scheduler, n_epochs, cuda, val_loader=triplets_train_loader_val, early_stopping=early_stopping, start_epoch_early_stopping=start_epoch_early_stopping, - epochs_early_stopping_interval=epochs_early_stopping_interval, start_epoch_training=start_epoch_training) + epochs_early_stopping_interval=epochs_early_stopping_interval, start_epoch_training=start_epoch_training, use_tqdm=use_tqdm) return embedding_net @@ -165,7 +166,7 @@ def triplet_training(samples, similarity_set, embedding_net, cuda, batch_size=16 def FP_nn_training(samples, target, embedding_net, cuda, batch_size=1, n_epochs=50, samples_val=None, target_val=None, early_stopping=False, epochs_early_stopping_interval=1, start_epoch_early_stopping=10, load_all_data_GPU=False, - lr=1e-3, optimizer=None, scheduler=None, start_epoch_training=0, optimizer_kwargs={}, + lr=1e-3, optimizer=None, scheduler=None, start_epoch_training=0, use_tqdm=True, optimizer_kwargs={}, scheduler_kwargs={}, loader_kwargs={}): """ Implements the algorithm for the training of a neural network based on regressing the values of the parameters on the corresponding simulation outcomes; it is effectively a training with a mean squared error loss. Needs to be @@ -224,6 +225,6 @@ def FP_nn_training(samples, target, embedding_net, cuda, batch_size=1, n_epochs= fit(data_loader_FP_nn, embedding_net, loss_fn, optimizer, scheduler, n_epochs, cuda, val_loader=data_loader_FP_nn_val, early_stopping=early_stopping, start_epoch_early_stopping=start_epoch_early_stopping, - epochs_early_stopping_interval=epochs_early_stopping_interval, start_epoch_training=start_epoch_training) + epochs_early_stopping_interval=epochs_early_stopping_interval, start_epoch_training=start_epoch_training, use_tqdm=use_tqdm) return embedding_net diff --git a/abcpy/NN_utilities/trainer.py b/abcpy/NN_utilities/trainer.py index a408a16a..12382efe 100644 --- a/abcpy/NN_utilities/trainer.py +++ b/abcpy/NN_utilities/trainer.py @@ -5,7 +5,7 @@ def fit(train_loader, model, loss_fn, optimizer, scheduler, n_epochs, cuda, val_loader=None, early_stopping=False, - epochs_early_stopping_interval=1, start_epoch_early_stopping=10, start_epoch_training=0): + epochs_early_stopping_interval=1, start_epoch_early_stopping=10, start_epoch_training=0, use_tqdm=True): """ Basic function to train a neural network given a train_loader, a loss function and an optimizer. @@ -26,7 +26,7 @@ def fit(train_loader, model, loss_fn, optimizer, scheduler, n_epochs, cuda, val_ for epoch in range(0, start_epoch_training): scheduler.step() - for epoch in tqdm(range(start_epoch_training, n_epochs)): + for epoch in tqdm(range(start_epoch_training, n_epochs), disable=not use_tqdm): # Train stage train_loss = train_epoch(train_loader, model, loss_fn, optimizer, cuda) diff --git a/abcpy/approx_lhd.py b/abcpy/approx_lhd.py index 4d3ddc34..99af55d6 100644 --- a/abcpy/approx_lhd.py +++ b/abcpy/approx_lhd.py @@ -1,14 +1,14 @@ -from abc import ABCMeta, abstractmethod - import numpy as np +from abc import ABCMeta, abstractmethod from glmnet import LogitNet +from scipy.stats import gaussian_kde, rankdata, norm from sklearn.covariance import ledoit_wolf from abcpy.graphtools import GraphTools class Approx_likelihood(metaclass=ABCMeta): - """This abstract base class defines the approximate likelihood + """This abstract base class defines the approximate likelihood function. """ @@ -16,15 +16,23 @@ class Approx_likelihood(metaclass=ABCMeta): def __init__(self, statistics_calc): """ The constructor of a sub-class must accept a non-optional statistics - calculator, which is stored to self.statistics_calc. + calculator; then, it must call the __init__ method of the parent class. This ensures that the + object is initialized correctly so that the _calculate_summary_stat private method can be called when computing + the distances. + Parameters ---------- - statistics_calc : abcpy.stasistics.Statistics + statistics_calc : abcpy.statistics.Statistics Statistics extractor object that conforms to the Statistics class. """ + self.statistics_calc = statistics_calc - raise NotImplemented + # Since the observations do always stay the same, we can save the + # summary statistics of them and not recalculate it each time + self.stat_obs = None + self.data_set = None + self.dataSame = False @abstractmethod def loglikelihood(self, y_obs, y_sim): @@ -48,32 +56,42 @@ def loglikelihood(self, y_obs, y_sim): raise NotImplemented def likelihood(self, y_obs, y_sim): - return np.exp(self.loglikelihood(y_obs, y_sim)) + """Computes the likelihood by taking the exponential of the loglikelihood method. + Parameters + ---------- + y_obs: Python list + Observed data set. + y_sim: Python list + Simulated data set from model at the parameter value. -class SynLikelihood(Approx_likelihood): - """This class implements the approximate likelihood function which computes the approximate - likelihood using the synthetic likelihood approach described in Wood [1]. - For synthetic likelihood approximation, we compute the robust precision matrix using Ledoit and Wolf's [2] - method. - - [1] S. N. Wood. Statistical inference for noisy nonlinear ecological - dynamic systems. Nature, 466(7310):1102–1104, Aug. 2010. - - [2] O. Ledoit and M. Wolf, A Well-Conditioned Estimator for Large-Dimensional Covariance Matrices, - Journal of Multivariate Analysis, Volume 88, Issue 2, pages 365-411, February 2004. - """ + Returns + ------- + float + Computed approximate likelihood. - def __init__(self, statistics_calc): - self.stat_obs = None - self.data_set = None - self.statistics_calc = statistics_calc + """ + return np.exp(self.loglikelihood(y_obs, y_sim)) - def loglikelihood(self, y_obs, y_sim): - # print("DEBUG: SynLikelihood.likelihood().") + def _calculate_summary_stat(self, y_obs, y_sim): + """Helper function that extracts the summary statistics s_obs and s_sim from y_obs and + y_y_sim using the statistics object stored in self.statistics_calc. This stores s_obs for the purpose of checking + whether that is repeated in next calls to the function, and avoiding computing the statitistics for the same + dataset several times. + + Parameters + ---------- + y_obs : array-like + d1 contains n_obs data sets. + y_sim : array-like + d2 contains n_sim data sets. + + Returns + ------- + tuple + Tuple containing numpy.ndarray's with the summary statistics extracted from d1 and d2. + """ if not isinstance(y_obs, list): - # print("type(y_obs) : ", type(y_obs), " , type(y_sim) : ", type(y_sim)) - # print("y_obs : ", y_obs) raise TypeError('Observed data is not of allowed types') if not isinstance(y_sim, list): @@ -81,6 +99,7 @@ def loglikelihood(self, y_obs, y_sim): # Check whether y_obs is same as the stored dataset. if self.data_set is not None: + # check that the the observations have the same length; if not, they can't be the same: if len(y_obs) != len(self.data_set): self.dataSame = False elif len(np.array(y_obs[0]).reshape(-1, )) == 1: @@ -96,6 +115,55 @@ def loglikelihood(self, y_obs, y_sim): # Extract summary statistics from the simulated data stat_sim = self.statistics_calc.statistics(y_sim) + if self.stat_obs.shape[1] != stat_sim.shape[1]: + raise ValueError("The dimension of summaries in the two datasets is different; check the dimension of the" + " provided observations and simulations.") + + return self.stat_obs, stat_sim + + +class SynLikelihood(Approx_likelihood): + + def __init__(self, statistics_calc): + """This class implements the approximate likelihood function which computes the approximate + likelihood using the synthetic likelihood approach described in Wood [1]. + For synthetic likelihood approximation, we compute the robust precision matrix using Ledoit and Wolf's [2] + method. + + [1] S. N. Wood. Statistical inference for noisy nonlinear ecological + dynamic systems. Nature, 466(7310):1102–1104, Aug. 2010. + + [2] O. Ledoit and M. Wolf, A Well-Conditioned Estimator for Large-Dimensional Covariance Matrices, + Journal of Multivariate Analysis, Volume 88, Issue 2, pages 365-411, February 2004. + + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + """ + + super(SynLikelihood, self).__init__(statistics_calc) + + def loglikelihood(self, y_obs, y_sim): + """Computes the loglikelihood. + + Parameters + ---------- + y_obs: Python list + Observed data set. + y_sim: Python list + Simulated data set from model at the parameter value. + + Returns + ------- + float + Computed approximate loglikelihood. + + """ + + stat_obs, stat_sim = self._calculate_summary_stat(y_obs, y_sim) + # Compute the mean, robust precision matrix and determinant of precision matrix mean_sim = np.mean(stat_sim, 0) lw_cov_, _ = ledoit_wolf(stat_sim) @@ -103,61 +171,181 @@ def loglikelihood(self, y_obs, y_sim): sign_logdet, robust_precision_sim_logdet = np.linalg.slogdet(robust_precision_sim) # we do not need sign # print("DEBUG: combining.") # we may have different observation; loop on those now: - # likelihoods = np.zeros(self.stat_obs.shape[0]) - # for i, single_stat_obs in enumerate(self.stat_obs): + # likelihoods = np.zeros(stat_obs.shape[0]) + # for i, single_stat_obs in enumerate(stat_obs): # x_new = np.einsum('i,ij,j->', single_stat_obs - mean_sim, robust_precision_sim, single_stat_obs - mean_sim) # likelihoods[i] = np.exp(-0.5 * x_new) # do without for loop: - diff = self.stat_obs - mean_sim.reshape(1, -1) + diff = stat_obs - mean_sim.reshape(1, -1) x_news = np.einsum('bi,ij,bj->b', diff, robust_precision_sim, diff) logliks = -0.5 * x_news - # looks like we are exponentiating the determinant as well, which is wrong; - # this is however a constant which should not change the algorithms afterwards. - logfactor = 0.5 * self.stat_obs.shape[0] * (np.log(1 / (2 * np.pi)) + robust_precision_sim_logdet) + logfactor = 0.5 * self.stat_obs.shape[0] * robust_precision_sim_logdet return np.sum(logliks) + logfactor # compute the sum of the different loglikelihoods for each observation +class SemiParametricSynLikelihood(Approx_likelihood): + + def __init__(self, statistics_calc, bw_method_marginals="silverman"): + """ + This class implements the approximate likelihood function which computes the approximate + likelihood using the semiparametric Synthetic Likelihood (semiBSL) approach described in [1]. Specifically, this + represents the likelihood as a product of univariate marginals and the copula components (exploiting Sklar's + theorem). + The marginals are approximated from simulations using a Gaussian KDE, while the copula is assumed to be a Gaussian + copula, whose parameters are estimated from data as well. + + This does not yet include shrinkage strategies for the correlation matrix. + + [1] An, Z., Nott, D. J., & Drovandi, C. (2020). Robust Bayesian synthetic likelihood via a semi-parametric approach. + Statistics and Computing, 30(3), 543-557. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + bw_method_marginals : str, scalar or callable, optional + The method used to calculate the estimator bandwidth, passed to `scipy.stats.gaussian_kde`. Following the docs + of that method, this can be 'scott', 'silverman', a scalar constant or a callable. If a scalar, this will be + used directly as `kde.factor`. If a callable, it should take a `gaussian_kde` instance as only parameter + and return a scalar. If None (default), 'silverman' is used. See the Notes in `scipy.stats.gaussian_kde` + for more details. + """ + super(SemiParametricSynLikelihood, self).__init__(statistics_calc) + # create a dict in which store the denominator of the correlation matrix for the different n values; + # this saves from repeating computations: + self.corr_matrix_denominator = {} + self.bw_method_marginals = bw_method_marginals # the bw method to use in the gaussian_kde + + def loglikelihood(self, y_obs, y_sim): + """Computes the loglikelihood. This implementation aims to be equivalent to the `BSL` R package, + but the results are slightly different due to small differences in the way the KDE is performed + + Parameters + ---------- + y_obs: Python list + Observed data set. + y_sim: Python list + Simulated data set from model at the parameter value. + + Returns + ------- + float + Computed approximate loglikelihood. + """ + + stat_obs, stat_sim = self._calculate_summary_stat(y_obs, y_sim) + n_obs, d = stat_obs.shape + if d < 2: + raise RuntimeError("The dimension of the statistics need to be at least 2 in order to apply semiBSL.") + + # first: estimate the marginal KDEs for each coordinate + logpdf_obs = np.zeros_like(stat_obs) # this will contain the estimated pdf at the various observation points + u_obs = np.zeros_like(stat_obs) # this instead will contain the transformed u's using the estimated CDF + for j in range(d): + # estimate the KDE using the data in stat_sim for coordinate j. This leads to slightly different results + # from the R package implementation due to slightly different way to estimate the factor as well as + # different way to evaluate the kernel (they use a weird interpolation there). + kde = gaussian_kde(stat_sim[:, j], bw_method=self.bw_method_marginals) + logpdf_obs[:, j] = kde.logpdf(stat_obs[:, j]) + for i in range(n_obs): # loop over the different observations + u_obs[i, j] = kde.integrate_box_1d(-np.infty, stat_obs[i, j]) # compute the CDF + etas_obs = norm.ppf(u_obs) + + # second: estimate the correlation matrix for the gaussian copula using gaussian rank correlation + R_hat = self._estimate_gaussian_correlation(stat_sim) + R_hat_inv = np.linalg.inv(R_hat) + R_sign_det, R_inv_logdet = np.linalg.slogdet(R_hat_inv) # sign not used + + # third: combine the two to compute the loglikelihood; + # for each observation: + # logliks = np.zeros(n_obs) + # for i in range(n_obs): + # logliks[i] = np.sum(logpdf_obs[i]) # sum along marginals along dimensions + # # add the copula density: + # logliks[i] += 0.5 * R_inv_logdet + # logliks[i] -= 0.5 * np.einsum("i,ij,j->", etas_obs[i], R_hat_inv - np.eye(d), etas_obs[i]) + + # do jointly: + loglik = np.sum(logpdf_obs) # sum along marginals along dimensions + # add the copula density: + copula_density = -0.5 * np.einsum("bi,ij,bj->b", etas_obs, R_hat_inv - np.eye(d), etas_obs) + loglik += np.sum(copula_density) + 0.5 * n_obs * R_inv_logdet + + return loglik + + def _estimate_gaussian_correlation(self, x): + """Estimates the correlation matrix using data in `x` in the way described in [1]. This implementation + gives the same results as the `BSL` R package. + + Parameters + ---------- + x: np.ndarray + Data set. + + Returns + ------- + np.ndarray + Estimated correlation matrix for the gaussian copula. + """ + n, d = x.shape + r = np.zeros_like(x) + for j in range(d): + r[:, j] = rankdata(x[:, j]) + + rqnorm = norm.ppf(r / (n + 1)) + + if n not in self.corr_matrix_denominator.keys(): + # compute the denominator: + self.corr_matrix_denominator[n] = np.sum(norm.ppf((np.arange(n) + 1) / (n + 1)) ** 2) + denominator = self.corr_matrix_denominator[n] + + R_hat = np.einsum('ki,kj->ij', rqnorm, rqnorm) / denominator + + return R_hat + + class PenLogReg(Approx_likelihood, GraphTools): - """This class implements the approximate likelihood function which computes the approximate - likelihood up to a constant using penalized logistic regression described in - Dutta et. al. [1]. It takes one additional function handler defining the - true model and two additional parameters n_folds and n_simulate correspondingly defining number - of folds used to estimate prediction error using cross-validation and the number - of simulated dataset sampled from each parameter to approximate the likelihood - function. For lasso penalized logistic regression we use glmnet of Friedman et. - al. [2]. - - [1] Thomas, O., Dutta, R., Corander, J., Kaski, S., & Gutmann, M. U. (2020). - Likelihood-free inference by ratio estimation. Bayesian Analysis. - - [2] Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization - paths for generalized linear models via coordinate descent. Journal of Statistical - Software, 33(1), 1–22. - - Parameters - ---------- - statistics_calc : abcpy.statistics.Statistics - Statistics extractor object that conforms to the Statistics class. - model : abcpy.models.Model - Model object that conforms to the Model class. - n_simulate : int - Number of data points to simulate for the reference data set; this has to be the same as n_samples_per_param - when calling the sampler. The reference data set is generated by drawing parameters from the prior and samples - from the model when PenLogReg is instantiated. - n_folds: int, optional - Number of folds for cross-validation. The default value is 10. - max_iter: int, optional - Maximum passes over the data. The default is 100000. - seed: int, optional - Seed for the random number generator. The used glmnet solver is not - deterministic, this seed is used for determining the cv folds. The default value is - None. - """ def __init__(self, statistics_calc, model, n_simulate, n_folds=10, max_iter=100000, seed=None): + """This class implements the approximate likelihood function which computes the approximate + likelihood up to a constant using penalized logistic regression described in + Dutta et. al. [1]. It takes one additional function handler defining the + true model and two additional parameters n_folds and n_simulate correspondingly defining number + of folds used to estimate prediction error using cross-validation and the number + of simulated dataset sampled from each parameter to approximate the likelihood + function. For lasso penalized logistic regression we use glmnet of Friedman et. + al. [2]. + + [1] Thomas, O., Dutta, R., Corander, J., Kaski, S., & Gutmann, M. U. (2020). + Likelihood-free inference by ratio estimation. Bayesian Analysis. + + [2] Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization + paths for generalized linear models via coordinate descent. Journal of Statistical + Software, 33(1), 1–22. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + model : abcpy.models.Model + Model object that conforms to the Model class. + n_simulate : int + Number of data points to simulate for the reference data set; this has to be the same as n_samples_per_param + when calling the sampler. The reference data set is generated by drawing parameters from the prior and + samples from the model when PenLogReg is instantiated. + n_folds: int, optional + Number of folds for cross-validation. The default value is 10. + max_iter: int, optional + Maximum passes over the data. The default is 100000. + seed: int, optional + Seed for the random number generator. The used glmnet solver is not + deterministic, this seed is used for determining the cv folds. The default value is + None. + """ + + super(PenLogReg, self).__init__(statistics_calc) # call the super init to initialize correctly self.model = model - self.statistics_calc = statistics_calc self.n_folds = n_folds self.n_simulate = n_simulate self.seed = seed @@ -166,33 +354,23 @@ def __init__(self, statistics_calc, model, n_simulate, n_folds=10, max_iter=1000 # Simulate reference data and extract summary statistics from the reference data self.ref_data_stat = self._simulate_ref_data(rng=self.rng)[0] - self.stat_obs = None - self.data_set = None - def loglikelihood(self, y_obs, y_sim): - if not isinstance(y_obs, list): - raise TypeError('Observed data is not of allowed types') - - if not isinstance(y_sim, list): - raise TypeError('simulated data is not of allowed types') + """Computes the loglikelihood. - # Check whether y_obs is same as the stored dataset. - if self.data_set is not None: - # check that the the observations have the same length; if not, they can't be the same: - if len(y_obs) != len(self.data_set): - self.dataSame = False - elif len(np.array(y_obs[0]).reshape(-1, )) == 1: - self.dataSame = self.data_set == y_obs - else: # otherwise it fails when y_obs[0] is array - self.dataSame = all( - [(np.array(self.data_set[i]) == np.array(y_obs[i])).all() for i in range(len(y_obs))]) + Parameters + ---------- + y_obs: Python list + Observed data set. + y_sim: Python list + Simulated data set from model at the parameter value. - if self.stat_obs is None or self.dataSame is False: - self.stat_obs = self.statistics_calc.statistics(y_obs) - self.data_set = y_obs + Returns + ------- + float + Computed approximate loglikelihood. + """ + stat_obs, stat_sim = self._calculate_summary_stat(y_obs, y_sim) - # Extract summary statistics from the simulated data - stat_sim = self.statistics_calc.statistics(y_sim) if not stat_sim.shape[0] == self.n_simulate: raise RuntimeError("The number of samples in the reference data set is not the same as the number of " "samples in the generated data. Please check that `n_samples` in the `sample()` method" @@ -207,7 +385,7 @@ def loglikelihood(self, y_obs, y_sim): groups += groups # duplicate it as groups need to be defined for both datasets m = LogitNet(alpha=1, n_splits=self.n_folds, max_iter=self.max_iter, random_state=self.seed, scoring="log_loss") m = m.fit(X, y, groups=groups) - result = -np.sum((m.intercept_ + np.sum(np.multiply(m.coef_, self.stat_obs), axis=1)), axis=0) + result = -np.sum((m.intercept_ + np.sum(np.multiply(m.coef_, stat_obs), axis=1)), axis=0) return result @@ -215,6 +393,17 @@ def _simulate_ref_data(self, rng=np.random.RandomState()): """ Simulate the reference data set. This code is run at the initialization of Penlogreg + + Parameters + ---------- + rng: Random number generator, optional + Defines the random number generator to be used. If None, a newly initialized one is used + + Returns + ------- + list + The simulated list of datasets. + """ ref_data_stat = [[None] * self.n_simulate for i in range(len(self.model))] diff --git a/abcpy/continuousmodels.py b/abcpy/continuousmodels.py index 898fb592..b1b6e2bb 100644 --- a/abcpy/continuousmodels.py +++ b/abcpy/continuousmodels.py @@ -1,6 +1,6 @@ import numpy as np from scipy.special import gamma -from scipy.stats import multivariate_normal, norm +from scipy.stats import multivariate_normal, norm, lognorm, expon from abcpy.probabilisticmodels import ProbabilisticModel, Continuous, InputConnector @@ -582,3 +582,195 @@ def pdf(self, input_values, x): density = normalizing_const * pow(tmp, -((df + p) / 2.)) self.calculated_pdf = density return density + + +class LogNormal(ProbabilisticModel, Continuous): + def __init__(self, parameters, name='LogNormal'): + """ + This class implements a probabilistic model following a Lognormal distribution with mean mu and variance sigma. + + Parameters + ---------- + parameters: list + Contains the probabilistic models and hyperparameters from which the model derives. + The list has two entries: from the first entry mean of the underlying normal distribution and from the second entry variance of the underlying normal + distribution is derived. + Note that the second value of the list is strictly greater than 0. + + name: string + The name that should be given to the probabilistic model in the journal file. + """ + + if not isinstance(parameters, list): + raise TypeError('Input for LogNormal has to be of type list.') + if len(parameters) < 2: + raise ValueError('Input for LogNormal has to be of length 2.') + + input_parameters = InputConnector.from_list(parameters) + super(LogNormal, self).__init__(input_parameters, name) + self.visited = False + + def _check_input(self, input_values): + """ + Returns True if the standard deviation is negative. + """ + if len(input_values) != 2: + return False + + if input_values[1] <= 0: + return False + return True + + def _check_output(self, parameters): + """ + Checks parameter values that are given as fixed values. + """ + return True + + def forward_simulate(self, input_values, k, rng=np.random.RandomState(), mpi_comm=None): + """ + Samples from a normal distribution using the current values for each probabilistic model from which the model derives. + + Parameters + ---------- + input_values: list + List of input parameters, in the same order as specified in the InputConnector passed to the init function + k: integer + The number of samples that should be drawn. + rng: Random number generator + Defines the random number generator to be used. The default value uses a random seed to initialize the generator. + + Returns + ------- + list: [np.ndarray] + A list containing the sampled values as np-array. + """ + + mu = input_values[0] + sigma = input_values[1] + result = np.array(rng.lognormal(mu, sigma, k)) + return [np.array([x]).reshape(-1, ) for x in result] + + def get_output_dimension(self): + return 1 + # Why does the following not work here? + # return self._dimension + + def pdf(self, input_values, x): + """ + Calculates the probability density function at point x. + Commonly used to determine whether perturbed parameters are still valid according to the pdf. + + Parameters + ---------- + input_values: list + List of input parameters of the from [mu, sigma] + x: list + The point at which the pdf should be evaluated. + + Returns + ------- + Float: + The evaluated pdf at point x. + """ + + mu = input_values[0] + sigma = input_values[1] + pdf = lognorm(scale=np.exp(mu), s=sigma).pdf(x) + self.calculated_pdf = pdf + return pdf + + +class Exponential(ProbabilisticModel, Continuous): + def __init__(self, parameters, name='Exponential'): + """ + This class implements a probabilistic model following a normal distribution with mean mu and variance sigma. + + Parameters + ---------- + parameters: list + Contains the probabilistic models and hyperparameters from which the model derives. + The list has one entry: the rate $\lambda$ of the exponential distribution, that has therefore pdf: + f(x; \lambda) = \lambda \exp(-\lambda x ), + name: string + The name that should be given to the probabilistic model in the journal file. + """ + + if not isinstance(parameters, list): + raise TypeError('Input for Exponential has to be of type list.') + if len(parameters) != 1: + raise ValueError('Input for Exponential has to be of length 1.') + + input_parameters = InputConnector.from_list(parameters) + super(Exponential, self).__init__(input_parameters, name) + self.visited = False + + def _check_input(self, input_values): + """ + Returns True if the standard deviation is negative. + """ + if len(input_values) != 1: + return False + + if input_values[0] <= 0: + return False + return True + + def _check_output(self, parameters): + """ + Checks parameter values that are given as fixed values. + """ + return True + + def forward_simulate(self, input_values, k, rng=np.random.RandomState(), mpi_comm=None): + """ + Samples from a normal distribution using the current values for each probabilistic model from which the model derives. + + Parameters + ---------- + input_values: list + List of input parameters, in the same order as specified in the InputConnector passed to the init function + k: integer + The number of samples that should be drawn. + rng: Random number generator + Defines the random number generator to be used. The default value uses a random seed to initialize the generator. + + Returns + ------- + list: [np.ndarray] + A list containing the sampled values as np-array. + """ + + rate = input_values[0] + scale = 1 / rate + result = np.array(rng.exponential(scale, k)) + return [np.array([x]).reshape(-1, ) for x in result] + + def get_output_dimension(self): + return 1 + # Why does the following not work here? + # return self._dimension + + def pdf(self, input_values, x): + """ + Calculates the probability density function at point x. + Commonly used to determine whether perturbed parameters are still valid according to the pdf. + + Parameters + ---------- + input_values: list + List of input parameters of the from [rate] + x: list + The point at which the pdf should be evaluated. + + Returns + ------- + Float: + The evaluated pdf at point x. + """ + + rate = input_values[0] + scale = 1 / rate + pdf = expon(scale=scale).pdf(x) + self.calculated_pdf = pdf + return pdf diff --git a/abcpy/distances.py b/abcpy/distances.py index 19a9e608..9d3b25a5 100644 --- a/abcpy/distances.py +++ b/abcpy/distances.py @@ -1,8 +1,9 @@ -from abc import ABCMeta, abstractmethod - import numpy as np +import warnings +from abc import ABCMeta, abstractmethod from glmnet import LogitNet from sklearn import linear_model +from sklearn.neighbors import NearestNeighbors from abcpy.utils import wass_dist @@ -130,24 +131,40 @@ def _calculate_summary_stat(self, d1, d2): s2 = self.statistics_calc.statistics(d2) + if self.s1.shape[1] != s2.shape[1]: + raise ValueError("The dimension of summaries in the two datasets is different; check the dimension of the" + " provided observations and simulations.") + return self.s1, s2 +class Divergence(Distance, metaclass=ABCMeta): + """This is an abstract class which subclasses Distance, and is used as a parent class for all divergence + estimators; more specifically, it is used for all Distances which compare the empirical distribution of simulations + and observations.""" + + @abstractmethod + def _estimate_always_positive(self): + """This returns whether the implemented divergence always returns positive values or not. In fact, some + estimators may return negative values, which may break some inference algorithms""" + raise NotImplementedError + + class Euclidean(Distance): """ This class implements the Euclidean distance between two vectors. The maximum value of the distance is np.inf. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + """ - def __init__(self, statistics): - """ - Parameters - ---------- - statistics_calc : abcpy.statistics.Statistics - Statistics extractor object that conforms to the Statistics class. - """ - super(Euclidean, self).__init__(statistics) + def __init__(self, statistics_calc): + super(Euclidean, self).__init__(statistics_calc) def distance(self, d1, d2): """Calculates the distance between two datasets, by computing Euclidean distance between each element of d1 and @@ -185,7 +202,7 @@ def dist_max(self): return np.inf -class PenLogReg(Distance): +class PenLogReg(Divergence): """ This class implements a distance measure based on the classification accuracy. @@ -202,16 +219,15 @@ class PenLogReg(Distance): [2] Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1), 1–22. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. """ - def __init__(self, statistics): - """ - Parameters - ---------- - statistics_calc : abcpy.statistics.Statistics - Statistics extractor object that conforms to the Statistics class. - """ - super(PenLogReg, self).__init__(statistics) + def __init__(self, statistics_calc): + super(PenLogReg, self).__init__(statistics_calc) self.n_folds = 10 # for cross validation in PenLogReg @@ -263,8 +279,11 @@ def dist_max(self): """ return 1.0 + def _estimate_always_positive(self): + return False + -class LogReg(Distance): +class LogReg(Divergence): """This class implements a distance measure based on the classification accuracy [1]. The classification accuracy is calculated between two dataset d1 and d2 using logistics regression and return it as a distance. The maximum value of the distance is 1.0. @@ -273,20 +292,18 @@ class LogReg(Distance): [1] Gutmann, M. U., Dutta, R., Kaski, S., & Corander, J. (2018). Likelihood-free inference via classification. Statistics and Computing, 28(2), 411-425. - """ - def __init__(self, statistics, seed=None): - """ - Parameters - ---------- - statistics_calc : abcpy.statistics.Statistics - Statistics extractor object that conforms to the Statistics class. - seed : integer, optionl - Seed used to initialize the Random Numbers Generator used to determine the (random) cross validation split - in the Logistic Regression classifier. - """ + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + seed : integer, optionl + Seed used to initialize the Random Numbers Generator used to determine the (random) cross validation split + in the Logistic Regression classifier. + """ - super(LogReg, self).__init__(statistics) + def __init__(self, statistics_calc, seed=None): + super(LogReg, self).__init__(statistics_calc) # seed is used for a RandomState for the random split in the LogisticRegression classifier: self.rng = np.random.RandomState(seed=seed) @@ -331,35 +348,101 @@ def dist_max(self): """ return 1.0 + def _estimate_always_positive(self): + return False -class Wasserstein(Distance): - """This class implements a distance measure based on the 2-Wasserstein distance, as used in [1]. This considers the - several simulations/observations in the datasets as iid samples from the model for a fixed parameter value/from the - data generating model, and computes the 2-Wasserstein distance between the empirical distributions those - simulations/observations define. + +class Wasserstein(Divergence): + """This class implements a distance measure based on the 2-Wasserstein distance, as used in [1]. This considers the + several simulations/observations in the datasets as iid samples from the model for a fixed parameter value/from the + data generating model, and computes the 2-Wasserstein distance between the empirical distributions those + simulations/observations define. [1] Bernton, E., Jacob, P.E., Gerber, M. and Robert, C.P. (2019), Approximate Bayesian computation with the Wasserstein distance. J. R. Stat. Soc. B, 81: 235-269. doi:10.1111/rssb.12312 + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + num_iter_max : integer, optional + The maximum number of iterations in the linear programming algorithm to estimate the Wasserstein distance. + Default to 100000. + """ - def __init__(self, statistics, num_iter_max=100000): - """ + def __init__(self, statistics_calc, num_iter_max=100000): + super(Wasserstein, self).__init__(statistics_calc) + + self.num_iter_max = num_iter_max + + def distance(self, d1, d2): + """Calculates the distance between two datasets. + Parameters ---------- - statistics_calc : abcpy.statistics.Statistics - Statistics extractor object that conforms to the Statistics class. - num_iter_max : integer, optional - The maximum number of iterations in the linear programming algorithm to estimate the Wasserstein distance. - Default to 100000. + d1: Python list + Contains n1 data points. + d2: Python list + Contains n2 data points. + + Returns + ------- + numpy.float + The distance between the two input data sets. + """ + s1, s2 = self._calculate_summary_stat(d1, d2) + + # compute the Wasserstein distance between the empirical distributions: + return wass_dist(samples_1=s1, samples_2=s2, num_iter_max=self.num_iter_max) + + def dist_max(self): + """ + Returns + ------- + numpy.float + The maximal possible value of the desired distance function. """ - super(Wasserstein, self).__init__(statistics) + # As the statistics are positive, the max possible value is 1 + return np.inf - self.num_iter_max = num_iter_max + def _estimate_always_positive(self): + return True + + +class SlicedWasserstein(Divergence): + """This class implements a distance measure based on the sliced 2-Wasserstein distance, as used in [1]. + This considers the several simulations/observations in the datasets as iid samples from the model for a fixed + parameter value/from the data generating model, and computes the sliced 2-Wasserstein distance between the + empirical distributions those simulations/observations define. Specifically, the sliced Wasserstein distance + is a cheaper version of the Wasserstein distance which consists of projecting the multivariate data on 1d directions + and computing the 1d Wasserstein distance, which is computationally cheap. The resulting sliced Wasserstein + distance is obtained by averaging over a given number of projections. + + [1] Nadjahi, K., De Bortoli, V., Durmus, A., Badeau, R., & Şimşekli, U. (2020, May). Approximate bayesian + computation with the sliced-wasserstein distance. In ICASSP 2020-2020 IEEE International Conference on Acoustics, + Speech and Signal Processing (ICASSP) (pp. 5470-5474). IEEE. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + n_projections : int, optional + Number of 1d projections used for estimating the sliced Wasserstein distance. Default value is 50. + rng : np.random.RandomState, optional + random number generators used to generate the projections. If not provided, a new one is instantiated. + """ + + def __init__(self, statistics_calc, n_projections=50, rng=np.random.RandomState()): + super(SlicedWasserstein, self).__init__(statistics_calc) + + self.n_projections = n_projections + self.rng = rng def distance(self, d1, d2): """Calculates the distance between two datasets. - + Parameters ---------- d1: Python list @@ -375,7 +458,592 @@ def distance(self, d1, d2): s1, s2 = self._calculate_summary_stat(d1, d2) # compute the Wasserstein distance between the empirical distributions: - return wass_dist(samples_1=s1, samples_2=s2, num_iter_max=self.num_iter_max) + return self.sliced_wasserstein_distance(X_s=s1, X_t=s2, n_projections=self.n_projections, seed=self.rng) + + def dist_max(self): + """ + Returns + ------- + numpy.float + The maximal possible value of the desired distance function. + """ + + # As the statistics are positive, the max possible value is 1 + return np.inf + + def _estimate_always_positive(self): + return True + + # the following two functions are taken from + # https://github.com/PythonOT/POT/blob/78b44af2434f494c8f9e4c8c91003fbc0e1d4415/ot/sliced.py + # Author: Adrien Corenflos + # + # License: MIT License + @staticmethod + def get_random_projections(n_projections, d, seed=None): + r""" + Generates n_projections samples from the uniform on the unit sphere of dimension d-1: :math:`\mathcal{U}(\mathcal{S}^{d-1})` + Parameters + ---------- + n_projections : int + number of samples requested + d : int + dimension of the space + seed: int or RandomState, optional + Seed used for numpy random number generator + Returns + ------- + out: ndarray, shape (n_projections, d) + The uniform unit vectors on the sphere + Examples + -------- + >>> n_projections = 100 + >>> d = 5 + >>> projs = get_random_projections(n_projections, d) + >>> np.allclose(np.sum(np.square(projs), 1), 1.) # doctest: +NORMALIZE_WHITESPACE + True + """ + + if not isinstance(seed, np.random.RandomState): + random_state = np.random.RandomState(seed) + else: + random_state = seed + + projections = random_state.normal(0., 1., [n_projections, d]) + norm = np.linalg.norm(projections, ord=2, axis=1, keepdims=True) + projections = projections / norm + return projections + + def sliced_wasserstein_distance(self, X_s, X_t, a=None, b=None, n_projections=50, seed=None, log=False): + r""" + Computes a Monte-Carlo approximation of the 2-Sliced Wasserstein distance + .. math:: + \mathcal{SWD}_2(\mu, \nu) = \underset{\theta \sim \mathcal{U}(\mathbb{S}^{d-1})}{\mathbb{E}}[\mathcal{W}_2^2(\theta_\# \mu, \theta_\# \nu)]^{\frac{1}{2}} + where : + - :math:`\theta_\# \mu` stands for the pushforwars of the projection :math:`\mathbb{R}^d \ni X \mapsto \langle \theta, X \rangle` + Parameters + ---------- + X_s : ndarray, shape (n_samples_a, dim) + samples in the source domain + X_t : ndarray, shape (n_samples_b, dim) + samples in the target domain + a : ndarray, shape (n_samples_a,), optional + samples weights in the source domain + b : ndarray, shape (n_samples_b,), optional + samples weights in the target domain + n_projections : int, optional + Number of projections used for the Monte-Carlo approximation + seed: int or RandomState or None, optional + Seed used for numpy random number generator + log: bool, optional + if True, sliced_wasserstein_distance returns the projections used and their associated EMD. + Returns + ------- + cost: float + Sliced Wasserstein Cost + log : dict, optional + log dictionary return only if log==True in parameters + Examples + -------- + >>> n_samples_a = 20 + >>> reg = 0.1 + >>> X = np.random.normal(0., 1., (n_samples_a, 5)) + >>> sliced_wasserstein_distance(X, X, seed=0) # doctest: +NORMALIZE_WHITESPACE + 0.0 + References + ---------- + .. [31] Bonneel, Nicolas, et al. "Sliced and radon wasserstein barycenters of measures." Journal of Mathematical Imaging and Vision 51.1 (2015): 22-45 + """ + from ot.lp import emd2_1d + + X_s = np.asanyarray(X_s) + X_t = np.asanyarray(X_t) + + n = X_s.shape[0] + m = X_t.shape[0] + + if X_s.shape[1] != X_t.shape[1]: + raise ValueError( + "X_s and X_t must have the same number of dimensions {} and {} respectively given".format(X_s.shape[1], + X_t.shape[1])) + + if a is None: + a = np.full(n, 1 / n) + if b is None: + b = np.full(m, 1 / m) + + d = X_s.shape[1] + + projections = self.get_random_projections(n_projections, d, seed) + + X_s_projections = np.dot(projections, X_s.T) + X_t_projections = np.dot(projections, X_t.T) + + if log: + projected_emd = np.empty(n_projections) + else: + projected_emd = None + + res = 0. + + for i, (X_s_proj, X_t_proj) in enumerate(zip(X_s_projections, X_t_projections)): + emd = emd2_1d(X_s_proj, X_t_proj, a, b, log=False, dense=False) + if projected_emd is not None: + projected_emd[i] = emd + res += emd + + res = (res / n_projections) ** 0.5 + if log: + return res, {"projections": projections, "projected_emds": projected_emd} + return res + + +class GammaDivergence(Divergence): + """ + This implements an empirical estimator of the gamma-divergence for ABC as suggested in [1]. In [1], the + gamma-divergence was proposed as a divergence which is robust to outliers. The estimator is based on a nearest + neighbor density estimate. + Specifically, this considers the + several simulations/observations in the datasets as iid samples from the model for a fixed parameter value/from the + data generating model, and estimates the divergence between the empirical distributions those + simulations/observations define. + + [1] Fujisawa, M., Teshima, T., Sato, I., & Sugiyama, M. + γ-ABC: Outlier-robust approximate Bayesian computation based on a + robust divergence estimator. + In A. Banerjee and K. Fukumizu (Eds.), Proceedings of 24th + International Conference on Artificial Intelligence and Statistics + (AISTATS2021), Proceedings of Machine Learning Research, vol.130, + pp.1783-1791, online, Apr. 13-15, 2021. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + k : int, optional + nearest neighbor number for the density estimate. Default value is 1 + gam : float, optional + the gamma parameter in the definition of the divergence. Default value is 0.1 + + """ + + def __init__(self, statistics_calc, k=1, gam=0.1): + super(GammaDivergence, self).__init__(statistics_calc) + + self.k = k # number of nearest neighbors used in the estimation algorithm + self.gam = gam + + def distance(self, d1, d2): + """Calculates the distance between two datasets. + + Parameters + ---------- + d1: Python list + Contains n1 data points. + d2: Python list + Contains n2 data points. + + Returns + ------- + numpy.float + The distance between the two input data sets. + """ + s1, s2 = self._calculate_summary_stat(d1, d2) + + if s1.shape[0] > self.k or s2.shape[0] > self.k: + assert ValueError(f"The provided value of k ({self.k}) is smaller or equal than the number of samples " + f"in one of the two datasets; that should instead be larger") + + # estimate the gamma divergence using the empirical distributions + return self.skl_estimator_gamma_q(s1=s1, s2=s2, k=self.k, gam=self.gam) + + def dist_max(self): + """ + Returns + ------- + numpy.float + The maximal possible value of the desired distance function. + """ + + # As the statistics are positive, the max possible value is 1 + return np.inf + + @staticmethod + def skl_estimator_gamma_q(s1, s2, k=1, gam=0.1): + """ Gamma-Divergence estimator using scikit-learn's NearestNeighbours + s1: (N_1,D) Sample drawn from distribution P + s2: (N_2,D) Sample drawn from distribution Q + k: Number of neighbours considered (default 1) + return: estimated D(P|Q) + + Adapted from code provided by Masahiro Fujisawa (University of Tokyo / RIKEN AIP) + """ + + n, m = len(s1), len(s2) # NOTE: here different convention of n, m wrt MMD and EnergyDistance + d = float(s1.shape[1]) + + radius = 10 # this is not used at all... + s1_neighbourhood = NearestNeighbors(n_neighbors=k + 1, radius=radius, algorithm='kd_tree').fit(s1) + s2_neighbourhood = NearestNeighbors(n_neighbors=k, radius=radius, algorithm='kd_tree').fit(s2) + s3_neighbourhood = NearestNeighbors(n_neighbors=k + 1, radius=radius, algorithm='kd_tree').fit(s2) + + d_gam = d * gam + + s1_distances, indices = s1_neighbourhood.kneighbors(s1, k + 1) + s2_distances, indices = s2_neighbourhood.kneighbors(s1, k) + rho = s1_distances[:, -1] + nu = s2_distances[:, -1] + if np.any(rho == 0): + warnings.warn( + f"The distance between an element of the first dataset and its {k}-th NN in the same dataset " + f"is 0; this causes divergences in the code, and it is due to elements which are repeated " + f"{k + 1} times in the first dataset. Increasing the value of k usually solves this.", + RuntimeWarning) + # notice: the one below becomes 0 when one element in the s1 dataset is equal to one in the s2 dataset + # and k=1 (as the distance between those two would be 0, which gives infinity when dividing) + if np.any(nu == 0): + warnings.warn(f"The distance between an element of the first dataset and its {k}-th NN in the second " + f"dataset is 0; this causes divergences in the code, and it is usually due to equal " + f"elements" + f" in the two datasets. Increasing the value of k usually solves this.", RuntimeWarning) + second_term = np.sum(1 / (rho ** d_gam)) / (n * (n - 1) ** gam) + fourth_term = np.sum(1 / (nu ** d_gam)) / (n * m ** gam) + + s3_distances, indices = s3_neighbourhood.kneighbors(s2, k + 1) + rho_q = s3_distances[:, -1] + + if np.any(rho_q == 0): + warnings.warn( + f"The distance between an element of the second dataset and its {k}-th NN in the same dataset " + f"is 0; this causes divergences in the code, and it is due to elements which are repeated " + f"{k + 1} times in the second dataset. Increasing the value of k usually solves this.", + RuntimeWarning) + + third_term = np.sum(1 / (rho_q ** d_gam)) + # third_term /= m * (m ** gam) # original code: I think the second term here should be m - 1 + third_term /= m * (m - 1) ** gam # corrected version + + third_term = third_term ** gam + fourth_term = fourth_term ** (1 + gam) + D = (1 / (gam * (gam + 1))) * (np.log((second_term * third_term) / fourth_term)) + return D + + def _estimate_always_positive(self): + return False + + +class KLDivergence(Divergence): + """ + This implements an empirical estimator of the KL divergence for ABC as suggested in [1]. The estimator is based + on a nearest neighbor density estimate. + Specifically, this considers the + several simulations/observations in the datasets as iid samples from the model for a fixed parameter value/from the + data generating model, and estimates the divergence between the empirical distributions those + simulations/observations define. + + [1] Jiang, B. (2018, March). Approximate Bayesian computation with Kullback-Leibler divergence as data discrepancy. + In International Conference on Artificial Intelligence and Statistics (pp. 1711-1721). PMLR. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + k : int, optional + nearest neighbor number for the density estimate. Default value is 1 + """ + + def __init__(self, statistics_calc, k=1): + super(KLDivergence, self).__init__(statistics_calc) + + self.k = k # number of nearest neighbors used in the estimation algorithm + + def distance(self, d1, d2): + """Calculates the distance between two datasets. + + Parameters + ---------- + d1: Python list + Contains n1 data points. + d2: Python list + Contains n2 data points. + + Returns + ------- + numpy.float + The distance between the two input data sets. + """ + s1, s2 = self._calculate_summary_stat(d1, d2) + + if s1.shape[0] > self.k or s2.shape[0] > self.k: + assert ValueError(f"The provided value of k ({self.k}) is smaller or equal than the number of samples " + f"in one of the two datasets; that should instead be larger") + + # estimate the KL divergence using the empirical distributions + return self.skl_estimator_KL_div(s1=s1, s2=s2, k=self.k) + + def dist_max(self): + """ + Returns + ------- + numpy.float + The maximal possible value of the desired distance function. + """ + + # As the statistics are positive, the max possible value is 1 + return np.inf + + @staticmethod + def skl_estimator_KL_div(s1, s2, k=1): + """ + Adapted from https://github.com/nhartland/KL-divergence-estimators/blob/5473a23f5f13d7557100504611c57c9225b1a6eb/src/knn_divergence.py + + MIT license + + KL-Divergence estimator using scikit-learn's NearestNeighbours + s1: (N_1,D) Sample drawn from distribution P + s2: (N_2,D) Sample drawn from distribution Q + k: Number of neighbours considered (default 1) + return: estimated D(P|Q) + """ + + n, m = len(s1), len(s2) # NOTE: here different convention of n, m wrt MMD and EnergyDistance + d = float(s1.shape[1]) + + radius = 10 # this is useless + s1_neighbourhood = NearestNeighbors(n_neighbors=k + 1, radius=radius, algorithm='kd_tree').fit(s1) + s2_neighbourhood = NearestNeighbors(n_neighbors=k, radius=radius, algorithm='kd_tree').fit(s2) + + s1_distances, indices = s1_neighbourhood.kneighbors(s1, k + 1) + s2_distances, indices = s2_neighbourhood.kneighbors(s1, k) + rho = s1_distances[:, -1] + nu = s2_distances[:, -1] + if np.any(rho == 0): + warnings.warn( + f"The distance between an element of the first dataset and its {k}-th NN in the same dataset " + f"is 0; this causes divergences in the code, and it is due to elements which are repeated " + f"{k + 1} times in the first dataset. Increasing the value of k usually solves this.", + RuntimeWarning) + D = np.sum(np.log(nu / rho)) + + return (d / n) * D + np.log(m / (n - 1)) # this second term should be enough for it to be valid for m \neq n + + def _estimate_always_positive(self): + return False + + +class MMD(Divergence): + """ + This implements an empirical estimator of the MMD for ABC as suggested in [1]. This class implements a gaussian + kernel by default but allows specifying different kernel functions. Notice that the original version in [1] + suggested an unbiased estimate, which however can return negative values. We also provide a biased but provably + positive estimator following the remarks in [2]. + Specifically, this considers the + several simulations/observations in the datasets as iid samples from the model for a fixed parameter value/from the + data generating model, and estimates the MMD between the empirical distributions those + simulations/observations define. + + [1] Park, M., Jitkrittum, W., & Sejdinovic, D. (2016, May). K2-ABC: Approximate Bayesian computation with + kernel embeddings. In Artificial Intelligence and Statistics (pp. 398-407). PMLR. + [2] Nguyen, H. D., Arbel, J., Lü, H., & Forbes, F. (2020). Approximate Bayesian computation via the energy + statistic. IEEE Access, 8, 131683-131698. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + kernel : str or callable + Can be a string denoting the kernel, or a function. If a string, only gaussian is implemented for now; in + that case, you can also provide an additional keyword parameter 'sigma' which is used as the sigma in the + kernel. Default is the gaussian kernel. + biased_estimator : boolean, optional + Whether to use the biased (but always positive) or unbiased estimator; by default, it uses the biased one. + kernel_kwargs + Additional keyword arguments to be passed to the distance calculator. + """ + + def __init__(self, statistics_calc, kernel="gaussian", biased_estimator=False, **kernel_kwargs): + super(MMD, self).__init__(statistics_calc) + + self.kernel_vectorized = False + if not isinstance(kernel, str) and not callable(kernel): + raise RuntimeError("'kernel' must be either a string or a function.") + if isinstance(kernel, str): + if kernel == "gaussian": + self.kernel = self.def_gaussian_kernel(**kernel_kwargs) + self.kernel_vectorized = True # the gaussian kernel is vectorized + else: + raise NotImplementedError("The required kernel is not implemented.") + else: + self.kernel = kernel # if kernel is a callable already + + self.biased_estimator = biased_estimator + + def distance(self, d1, d2): + """Calculates the distance between two datasets. + + Parameters + ---------- + d1: Python list + Contains n1 data points. + d2: Python list + Contains n2 data points. + + Returns + ------- + numpy.float + The distance between the two input data sets. + """ + s1, s2 = self._calculate_summary_stat(d1, d2) + + # compute the Gram matrix + K11, K22, K12 = self.compute_Gram_matrix(s1, s2) + + # Estimate MMD + if self.biased_estimator: + return self.MMD_V_estimator(K11, K22, K12) + else: + return self.MMD_unbiased(K11, K22, K12) + + def dist_max(self): + """ + Returns + ------- + numpy.float + The maximal possible value of the desired distance function. + """ + + # As the statistics are positive, the max possible value is 1 + return np.inf + + @staticmethod + def def_gaussian_kernel(sigma=1): + # notice in the MMD paper they set sigma to a median value over the observation; check that. + sigma_2 = 2 * sigma ** 2 + + # def Gaussian_kernel(x, y): + # xy = x - y + # # assert np.allclose(np.dot(xy, xy), np.linalg.norm(xy) ** 2) + # return np.exp(- np.dot(xy, xy) / sigma_2) + def Gaussian_kernel_vectorized(X, Y): + """Here X and Y have shape (n_samples_x, n_features) and (n_samples_y, n_features); + this directly computes the kernel for all pairwise components""" + XY = X.reshape(X.shape[0], 1, -1) - Y.reshape(1, Y.shape[0], -1) # pairwise differences + return np.exp(- np.einsum('xyi,xyi->xy', XY, XY) / sigma_2) + + return Gaussian_kernel_vectorized + + def compute_Gram_matrix(self, s1, s2): + + if self.kernel_vectorized: + K11 = self.kernel(s1, s1) + K22 = self.kernel(s2, s2) + K12 = self.kernel(s1, s2) + else: + m = s1.shape[0] + n = s2.shape[0] + + K11 = np.zeros((m, m)) + K22 = np.zeros((n, n)) + K12 = np.zeros((m, n)) + + for i in range(m): + # we assume the function to be symmetric; this saves some steps: + for j in range(i, m): + K11[j, i] = K11[i, j] = self.kernel(s1[i], s1[j]) + + for i in range(n): + # we assume the function to be symmetric; this saves some steps: + for j in range(i, n): + K22[j, i] = K22[i, j] = self.kernel(s2[i], s2[j]) + + for i in range(m): + for j in range(n): + K12[i, j] = self.kernel(s1[i], s2[j]) + + # can we improve the above? Could use map but would not change too much likely. + return K11, K22, K12 + + @staticmethod + def MMD_unbiased(Kxx, Kyy, Kxy): + # from https://github.com/eugenium/MMD/blob/2fe67cbc7378f10f3b273cfd8d8bbd2135db5798/mmd.py + # The estimate when distribution of x is not equal to y + m = Kxx.shape[0] + n = Kyy.shape[0] + + t1 = (1. / (m * (m - 1))) * np.sum(Kxx - np.diag(np.diagonal(Kxx))) + t2 = (2. / (m * n)) * np.sum(Kxy) + t3 = (1. / (n * (n - 1))) * np.sum(Kyy - np.diag(np.diagonal(Kyy))) + + MMDsquared = (t1 - t2 + t3) + + return MMDsquared + + @staticmethod + def MMD_V_estimator(Kxx, Kyy, Kxy): + # The estimate when distribution of x is not equal to y + m = Kxx.shape[0] + n = Kyy.shape[0] + + t1 = (1. / (m * m)) * np.sum(Kxx) + t2 = (2. / (m * n)) * np.sum(Kxy) + t3 = (1. / (n * n)) * np.sum(Kyy) + + MMDsquared = (t1 - t2 + t3) + + return MMDsquared + + def _estimate_always_positive(self): + return self.biased_estimator + + +class EnergyDistance(MMD): + """ + This implements an empirical estimator of the Energy Distance for ABC as suggested in [1]. + This class uses the Euclidean distance by default as a base distance, but allows to pass different distances. + Moreover, when the Euclidean distance is specified, it is possible to pass an additional keyword argument `beta` + which denotes the power of the distance to consider. + In [1], the authors suggest to use a biased but provably positive estimator; we also provide an unbiased estimate, + which however can return negative values. + Specifically, this considers the + several simulations/observations in the datasets as iid samples from the model for a fixed parameter value/from the + data generating model, and estimates the MMD between the empirical distributions those + simulations/observations define. + + [1] Nguyen, H. D., Arbel, J., Lü, H., & Forbes, F. (2020). Approximate Bayesian computation via the energy + statistic. IEEE Access, 8, 131683-131698. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + base_distance : str or callable + Can be a string denoting the kernel, or a function. If a string, only Euclidean distance is implemented + for now; in that case, you can also provide an additional keyword parameter 'beta' which is the power + of the distance to consider. By default, this uses the Euclidean distance. + biased_estimator : boolean, optional + Whether to use the biased (but always positive) or unbiased estimator; by default, it uses the biased one. + base_distance_kwargs + Additional keyword arguments to be passed to the distance calculator. + """ + + def __init__(self, statistics_calc, base_distance="Euclidean", biased_estimator=True, **base_distance_kwargs): + if not isinstance(base_distance, str) and not callable(base_distance): + raise RuntimeError("'base_distance' must be either a string or a function.") + if isinstance(base_distance, str): + if base_distance == "Euclidean": + self.base_distance = self.def_Euclidean_distance(**base_distance_kwargs) + else: + raise NotImplementedError("The required kernel is not implemented.") + else: + self.base_distance = base_distance # if base_distance is a callable already + + self.biased_estimator = biased_estimator + + def negative_distance(*args): + return - self.base_distance(*args) + + super(EnergyDistance, self).__init__(statistics_calc, kernel=negative_distance, + biased_estimator=self.biased_estimator) def dist_max(self): """ @@ -387,3 +1055,147 @@ def dist_max(self): # As the statistics are positive, the max possible value is 1 return np.inf + + @staticmethod + def def_Euclidean_distance(beta=1): + if beta <= 0 or beta > 2: + raise RuntimeError("'beta' not in the right range (0,2]") + + if beta == 1: + def Euclidean_distance(x, y): + return np.linalg.norm(x - y) + else: + def Euclidean_distance(x, y): + return np.linalg.norm(x - y) ** beta + + return Euclidean_distance + + +class SquaredHellingerDistance(Divergence): + """ + This implements an empirical estimator of the squared Hellinger distance for ABC. Using the Hellinger distance was + suggested originally in [1], but as that work did not provide originally any implementation details, this + implementation is original. The estimator is based on a nearest neighbor density estimate. + Specifically, this considers the + several simulations/observations in the datasets as iid samples from the model for a fixed parameter value/from the + data generating model, and estimates the divergence between the empirical distributions those + simulations/observations define. + + [1] Frazier, D. T. (2020). Robust and Efficient Approximate Bayesian Computation: A Minimum Distance Approach. + arXiv preprint arXiv:2006.14126. + + Parameters + ---------- + statistics_calc : abcpy.statistics.Statistics + Statistics extractor object that conforms to the Statistics class. + k : int, optional + nearest neighbor number for the density estimate. Default value is 1 + """ + + def __init__(self, statistics_calc, k=1): + super(SquaredHellingerDistance, self).__init__(statistics_calc) + + self.k = k # number of nearest neighbors used in the estimation algorithm + + def distance(self, d1, d2): + """Calculates the distance between two datasets. + + Parameters + ---------- + d1: Python list + Contains n1 data points. + d2: Python list + Contains n2 data points. + + Returns + ------- + numpy.float + The distance between the two input data sets. + """ + s1, s2 = self._calculate_summary_stat(d1, d2) + + if s1.shape[0] > self.k or s2.shape[0] > self.k: + assert ValueError(f"The provided value of k ({self.k}) is smaller or equal than the number of samples " + f"in one of the two datasets; that should instead be larger") + + # estimate the gamma divergence using the empirical distributions + return self.skl_estimator_squared_Hellinger_distance(s1=s1, s2=s2, k=self.k) + + def dist_max(self): + """ + Returns + ------- + numpy.float + The maximal possible value of the desired distance function. + """ + + return 2 + + @staticmethod + def skl_estimator_squared_Hellinger_distance(s1, s2, k=1): + """ Squared Hellinger distance estimator using scikit-learn's NearestNeighbours + s1: (N_1,D) Sample drawn from distribution P + s2: (N_2,D) Sample drawn from distribution Q + k: Number of neighbours considered (default 1) + return: estimated D(P|Q) + + """ + + n, m = len(s1), len(s2) # NOTE: here different convention of n, m wrt MMD and EnergyDistance + d = float(s1.shape[1]) + d_2 = d / 2 + + radius = 10 # this is not used at all... + s1_neighbourhood_k1 = NearestNeighbors(n_neighbors=k + 1, radius=radius, algorithm='kd_tree').fit(s1) + s1_neighbourhood_k = NearestNeighbors(n_neighbors=k, radius=radius, algorithm='kd_tree').fit(s1) + s2_neighbourhood_k1 = NearestNeighbors(n_neighbors=k + 1, radius=radius, algorithm='kd_tree').fit(s2) + s2_neighbourhood_k = NearestNeighbors(n_neighbors=k, radius=radius, algorithm='kd_tree').fit(s2) + + s1_distances, indices = s1_neighbourhood_k1.kneighbors(s1, k + 1) + s2_distances, indices = s2_neighbourhood_k.kneighbors(s1, k) + rho = s1_distances[:, -1] + nu = s2_distances[:, -1] + if np.any(rho == 0): + warnings.warn( + f"The distance between an element of the first dataset and its {k}-th NN in the same dataset " + f"is 0; this is due to elements which are repeated " + f"{k + 1} times in the first dataset, and may lead to a poor estimate of the distance. " + f"Increasing the value of k usually solves this.", + RuntimeWarning) + + if np.any(nu == 0): + warnings.warn(f"The distance between an element of the first dataset and its {k}-th NN in the second " + f"dataset is 0; this causes divergences in the code, and it is usually due to equal " + f"elements" + f" in the two datasets. Increasing the value of k usually solves this.", RuntimeWarning) + first_estimator = np.sum((rho / nu) ** d_2) + first_estimator = 2 - 2 * np.sqrt((n - 1) / m) * first_estimator + + s2_distances, indices = s2_neighbourhood_k1.kneighbors(s2, k + 1) + s1_distances, indices = s1_neighbourhood_k.kneighbors(s2, k) + rho = s2_distances[:, -1] + nu = s1_distances[:, -1] + if np.any(rho == 0): + warnings.warn( + f"The distance between an element of the second dataset and its {k}-th NN in the same dataset " + f"is 0; this is due to elements which are repeated " + f"{k + 1} times in the second dataset, and may lead to a poor estimate of the distance. " + f"Increasing the value of k usually solves this.", + RuntimeWarning) + # notice: the one below becomes 0 when one element in the s1 dataset is equal to one in the s2 dataset + # and k=1 (as the distance between those two would be 0, which gives infinity when dividing) + if np.any(nu == 0): + warnings.warn(f"The distance between an element of the second dataset and its {k}-th NN in the first " + f"dataset is 0; this causes divergences in the code, and it is usually due to equal " + f"elements" + f" in the two datasets. Increasing the value of k usually solves this.", RuntimeWarning) + second_estimator = np.sum((rho / nu) ** d_2) + second_estimator = 2 - 2 * np.sqrt((m - 1) / n) * second_estimator + + # average the two estimators: + final_estimator = 0.5 * (first_estimator + second_estimator) + + return final_estimator + + def _estimate_always_positive(self): + return True diff --git a/abcpy/inferences.py b/abcpy/inferences.py index 137302c8..922b1ca0 100644 --- a/abcpy/inferences.py +++ b/abcpy/inferences.py @@ -1,19 +1,24 @@ import copy import logging import time +import warnings from abc import abstractproperty import numpy as np from scipy import optimize +from tqdm import tqdm from abcpy.acceptedparametersmanager import * +from abcpy.backends import BackendDummy +from abcpy.distances import Divergence from abcpy.graphtools import GraphTools from abcpy.jointapprox_lhd import SumCombination from abcpy.jointdistances import LinearCombination from abcpy.output import Journal -from abcpy.perturbationkernel import DefaultKernel +from abcpy.perturbationkernel import DefaultKernel, JointPerturbationKernel from abcpy.probabilisticmodels import * from abcpy.utils import cached +from abcpy.transformers import BoundedVarTransformer, DummyTransformer class InferenceMethod(GraphTools, metaclass=ABCMeta): @@ -77,7 +82,7 @@ def kernel(self): """To be overwritten by any sub-class: an attribute specifying the transition or perturbation kernel.""" raise NotImplementedError - def perturb(self, column_index, epochs=10, rng=np.random.RandomState()): + def perturb(self, column_index, epochs=100, rng=np.random.RandomState(), accepted_parameters_manager=None): """ Perturbs all free parameters, given the current weights. Commonly used during inference. @@ -88,6 +93,9 @@ def perturb(self, column_index, epochs=10, rng=np.random.RandomState()): The index of the column in the accepted_parameters_bds that should be used for perturbation epochs: integer The number of times perturbation should happen before the algorithm is terminated + accepted_parameters_manager: AcceptedParametersManager + The AcceptedParametersManager to use; if not provided, use the one stored in + self.accepted_parameters_manager Returns ------- @@ -96,10 +104,13 @@ def perturb(self, column_index, epochs=10, rng=np.random.RandomState()): """ current_epoch = 0 + if accepted_parameters_manager is None: + accepted_parameters_manager = self.accepted_parameters_manager + while current_epoch < epochs: # Get new parameters of the graph - new_parameters = self.kernel.update(self.accepted_parameters_manager, column_index, rng=rng) + new_parameters = self.kernel.update(accepted_parameters_manager, column_index, rng=rng) self._reset_flags() @@ -149,10 +160,11 @@ class RejectionABC(InferenceMethod): Parameters ---------- - model: list + root_models: list A list of the Probabilistic models corresponding to the observed datasets - distance: abcpy.distances.Distance - Distance object defining the distance measure to compare simulated and observed data sets. + distances: list of abcpy.distances.Distance + List of Distance objects defining the distance measure to compare simulated and observed data sets; one for + each model. backend: abcpy.backends.Backend Backend object defining the backend to be used. seed: integer, optionaldistance @@ -184,44 +196,113 @@ def __init__(self, root_models, distances, backend, seed=None): # counts the number of simulate calls self.simulation_counter = 0 - def sample(self, observations, n_samples, n_samples_per_param, epsilon, full_output=0): + def sample(self, observations, n_samples=100, n_samples_per_param=1, epsilon=None, simulation_budget=None, + quantile=None, full_output=0, path_to_save_journal=None): """ Samples from the posterior distribution of the model parameter given the observed data observations. + You can either specify the required number of posterior samples `n_samples`, or the + `simulation_budget`. In the former case, the threshold value `epsilon` is required, so that the algorithm + will produce `n_samples` posterior samples for which the ABC distance was smaller than `epsilon`. In the latter + case, you can specify either `epsilon` or `quantile`; in this case, the number of simulations specified in the + simulation budget will be run, and only the parameter values for which the ABC distance was smaller than + `epsilon` (or alternatively the ones for which the ABC distance is in the smaller specified quantile) will + be returned. + Parameters ---------- observations: list - A list, containing lists describing the observed data sets + A list, containing lists describing the observed data sets; one for each model. n_samples: integer Number of samples to generate n_samples_per_param: integer Number of data points in each simulated data set. epsilon: float Value of threshold + simulation_budget : integer, optional + Simulation budget to be considered (ie number of parameter values for which the ABC distance is computed). + Alternative to `n_samples`, which needs to be set explicitly to None to use this. Defaults to None. + quantile : float, optional + If `simulation_budget` is used, only the samples which achieve performance less than the specified quantile + of the ABC distances, will be retained in the set of posterior samples. This is alternative to epsilon. + Defaults to None. full_output: integer, optional - If full_output==1, intermediate results are included in output journal. - The default value is 0, meaning the intermediate results are not saved. - + It is actually unused in RejectionABC but left here for general compatibility with the other inference + classes. + path_to_save_journal: str, optional + If provided, save the journal after inference at the provided path. + Returns ------- abcpy.output.Journal a journal containing simulation results, metadata and optionally intermediate results. """ + if (n_samples is None) == (simulation_budget is None): + raise RuntimeError("One and only one of `n_samples` and `simulation_budget` needs to be specified.") + if n_samples is not None and quantile is not None: + raise RuntimeError("`quantile` can be specified only when the simulation budget is fixed with " + "`simulation_budget`.") + if n_samples is not None and epsilon is None: + raise RuntimeError("`epsilon` needs to be specified when the `n_samples` is given ") + if simulation_budget is not None and ((quantile is None) == (epsilon is None)): + raise RuntimeError("One and only one of `quantile` and `epsilon` needs to be specified when " + "`simulation_budget` is used.") + + self.fixed_budget = True if simulation_budget is not None else False + self.accepted_parameters_manager.broadcast(self.backend, observations) - self.n_samples = n_samples self.n_samples_per_param = n_samples_per_param - self.epsilon = epsilon + # instantiate journal (common to both approaches); this will be overwritten if self.fixed_budget is True journal = Journal(full_output) - journal.configuration["n_samples"] = self.n_samples + journal.configuration["n_samples"] = n_samples journal.configuration["n_samples_per_param"] = self.n_samples_per_param - journal.configuration["epsilon"] = self.epsilon + journal.configuration["epsilon"] = epsilon + journal.configuration["type_model"] = [type(model).__name__ for model in self.model] + journal.configuration["type_dist_func"] = [type(distance).__name__ for distance in self.distance.distances] + journal.configuration["full_output"] = full_output - accepted_parameters = None + if self.fixed_budget: + # sample simulation_budget number of samples with super high epsilon + n_samples = simulation_budget + self.epsilon = 1.7976931348623157e+308 # max possible numpy value + + # call sample routine + journal = self._sample_n_samples_epsilon(n_samples, journal) + + # then replace that journal with a new one selecting the correct samples only, with either epsilon or + # quantile (only one of them will not be None) + journal = self._journal_cleanup(journal, quantile, epsilon) + else: + self.epsilon = epsilon + + # call sample routine + journal = self._sample_n_samples_epsilon(n_samples, journal) + + if path_to_save_journal is not None: + # save journal there + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else \ + path_to_save_journal + '.jnl' + journal.save(path_to_save_journal) + return journal + + def _sample_n_samples_epsilon(self, n_samples, journal): + """Obtains `n_samples` posterior samples with threshold `self.epsilon`, and stores them into the journal. + Parameters + ---------- + n_samples: integer + Number of samples to generate + journal: abcpy.output.Journal + Journal file where to store results + Returns + ------- + abcpy.output.Journal + a journal containing simulation results, metadata and optionally intermediate results. + """ # main Rejection ABC algorithm seed_arr = self.rng.randint(1, n_samples * n_samples, size=n_samples, dtype=np.int32) rng_arr = np.array([np.random.RandomState(seed) for seed in seed_arr]) @@ -248,6 +329,68 @@ def sample(self, observations, n_samples, n_samples_per_param, epsilon, full_out return journal + def _journal_cleanup(self, journal, quantile=None, threshold=None): + """This function takes a Journal file (typically produced by an Rejection ABC run with very large epsilon value) + and keeps only the samples which achieve performance less than either some quantile of the ABC distances, + or either some specified threshold. It is a very simple way to obtain a Rejection ABC which + works on a percentile of the obtained distances. + + It creates a new Journal file storing the results. + + Parameters + ---------- + journal: abcpy.output.Journal + Journal file where to store results + quantile : float, optional + If `simulation_budget` is used, only the samples which achieve performance less than the specified quantile + of the ABC distances, will be retained in the set of posterior samples. This is alternative to epsilon. + Defaults to None. + threshold: float + Value of threshold + + Returns + ------- + abcpy.output.Journal + a new journal containing simulation results, metadata and optionally intermediate results. + """ + + if quantile is not None: + distance_cutoff = np.quantile(journal.distances[-1], quantile) + else: + distance_cutoff = threshold + picked_simulations = journal.distances[-1] < distance_cutoff + new_distances = journal.distances[-1][picked_simulations] + n_reduced_samples = np.sum(picked_simulations) + if n_reduced_samples == 0: + raise RuntimeError( + "The specified value of threshold is too low, no simulations from the ones generated with the fixed " + "simulation budget are accepted." + ) + new_journal = Journal(journal._type) + new_journal.configuration["n_samples"] = n_reduced_samples + new_journal.configuration["n_samples_per_param"] = journal.configuration[ + "n_samples_per_param" + ] + new_journal.configuration["epsilon"] = distance_cutoff + + new_accepted_parameters = [] + param_names = journal.get_parameters().keys() + new_names_and_parameters = {name: [] for name in param_names} + for i in np.where(picked_simulations)[0]: + if picked_simulations[i]: + new_accepted_parameters.append(journal.get_accepted_parameters()[i]) + for name in param_names: + new_names_and_parameters[name].append(journal.get_parameters()[name][i]) + + new_journal.add_accepted_parameters(new_accepted_parameters) + new_journal.add_weights(np.ones((n_reduced_samples, 1))) + new_journal.add_ESS_estimate(np.ones((n_reduced_samples, 1))) + new_journal.add_distances(new_distances) + new_journal.add_user_parameters(new_names_and_parameters) + new_journal.number_of_simulations.append(journal.number_of_simulations[-1]) + + return new_journal + def _sample_parameter(self, rng, npc=None): """ Samples a single model parameter and simulates from it until @@ -305,14 +448,16 @@ class PMCABC(BaseDiscrepancy, InferenceMethod): Parameters ---------- - model : list + root_models: list A list of the Probabilistic models corresponding to the observed datasets - distance : abcpy.distances.Distance - Distance object defining the distance measure to compare simulated and observed data sets. - kernel : abcpy.distributions.Distribution - Distribution object defining the perturbation kernel needed for the sampling. + distances: list of abcpy.distances.Distance + List of Distance objects defining the distance measure to compare simulated and observed data sets; one for + each model. backend : abcpy.backends.Backend Backend object defining the backend to be used. + kernel : abcpy.perturbationkernel.PerturbationKernel, optional + PerturbationKernel object defining the perturbation kernel needed for the sampling. If not provided, the + DefaultKernel is used. seed : integer, optional Optional initial seed for the random number generator. The default value is generated randomly. """ @@ -350,16 +495,16 @@ def __init__(self, root_models, distances, backend, kernel=None, seed=None): self.simulation_counter = 0 def sample(self, observations, steps, epsilon_init, n_samples=10000, n_samples_per_param=1, epsilon_percentile=10, - covFactor=2, full_output=0, journal_file=None): + covFactor=2, full_output=0, journal_file=None, path_to_save_journal=None): """Samples from the posterior distribution of the model parameter given the observed data observations. Parameters ---------- observations : list - A list, containing lists describing the observed data sets + A list, containing lists describing the observed data sets; one for each model. steps : integer - Number of iterations in the sequential algoritm ("generations") + Number of iterations in the sequential algorithm ("generations") epsilon_init : numpy.ndarray An array of proposed values of epsilon to be used at each steps. Can be supplied A single value to be used as the threshold in Step 1 or a `steps`-dimensional array of values to be @@ -378,6 +523,11 @@ def sample(self, observations, steps, epsilon_init, n_samples=10000, n_samples_p journal_file: str, optional Filename of a journal file to read an already saved journal file, from which the first iteration will start. The default value is None. + path_to_save_journal: str, optional + If provided, save the journal at the provided path. The journal is saved (and overwritten) after each step + of the sequential inference routine, so that partial results are stored to the disk in case the + inference routine does not end correctly; recall that you need to set full_output=1 to obtain the + full partial results. Returns ------- @@ -388,14 +538,21 @@ def sample(self, observations, steps, epsilon_init, n_samples=10000, n_samples_p self.n_samples = n_samples self.n_samples_per_param = n_samples_per_param + if path_to_save_journal is not None: + # save journal there + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else path_to_save_journal + '.jnl' + if journal_file is None: journal = Journal(full_output) journal.configuration["type_model"] = [type(model).__name__ for model in self.model] - journal.configuration["type_dist_func"] = type(self.distance).__name__ + journal.configuration["type_dist_func"] = [type(distance).__name__ for distance in self.distance.distances] + journal.configuration["steps"] = steps + journal.configuration["epsilon_init"] = epsilon_init journal.configuration["n_samples"] = self.n_samples journal.configuration["n_samples_per_param"] = self.n_samples_per_param - journal.configuration["steps"] = steps journal.configuration["epsilon_percentile"] = epsilon_percentile + journal.configuration["covFactor"] = covFactor + journal.configuration["full_output"] = full_output else: journal = Journal.fromFile(journal_file) @@ -408,7 +565,7 @@ def sample(self, observations, steps, epsilon_init, n_samples=10000, n_samples_p epsilon_arr = epsilon_init elif len(epsilon_init) == 1: epsilon_arr = [None] * steps - epsilon_arr[0] = epsilon_init + epsilon_arr[0] = epsilon_init[0] else: raise ValueError("The length of epsilon_init can only be equal to 1 or steps.") @@ -420,7 +577,7 @@ def sample(self, observations, steps, epsilon_init, n_samples=10000, n_samples_p accepted_parameters = journal.get_accepted_parameters(-1) accepted_weights = journal.get_weights(-1) - if hasattr(journal, "distances"): + if hasattr(journal, "distances"): # if restarting from a journal, use the previous distances to check determine a new epsilon # (it if is larger than the epsilon_arr[0] provided here) self.logger.info("Calculating acceptances threshold from provided journal file") @@ -436,7 +593,7 @@ def sample(self, observations, steps, epsilon_init, n_samples=10000, n_samples_p self.accepted_parameters_manager.update_kernel_values(self.backend, kernel_parameters=kernel_parameters) # 3: calculate covariance - self.logger.info("Calculateing covariance matrix") + self.logger.info("Calculating covariance matrix") new_cov_mats = self.kernel.calculate_cov(self.accepted_parameters_manager) # Since each entry of new_cov_mats is a numpy array, we can multiply like this # accepted_cov_mats = [covFactor * new_cov_mat for new_cov_mat in new_cov_mats] @@ -523,11 +680,14 @@ def sample(self, observations, steps, epsilon_init, n_samples=10000, n_samples_p journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) - # Add epsilon_arr to the journal - if journal_file is not None and "epsilon_arr" in journal.configuration.keys(): - journal.configuration["epsilon_arr"] += epsilon_arr - else: - journal.configuration["epsilon_arr"] = epsilon_arr + # Add epsilon_arr to the journal + if journal_file is not None and "epsilon_arr" in journal.configuration.keys(): + journal.configuration["epsilon_arr"] += epsilon_arr + else: + journal.configuration["epsilon_arr"] = epsilon_arr + + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) return journal @@ -644,9 +804,9 @@ def _compute_accepted_cov_mats(self, covFactor, new_cov_mats): for new_cov_mat in new_cov_mats: if not (new_cov_mat.size == 1): accepted_cov_mats.append( - covFactor * new_cov_mat + 0.0001 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) + covFactor * new_cov_mat + 1e-20 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) else: - accepted_cov_mats.append((covFactor * new_cov_mat + 0.0001 * new_cov_mat).reshape(1, 1)) + accepted_cov_mats.append((covFactor * new_cov_mat + 1e-20 * new_cov_mat).reshape(1, 1)) return accepted_cov_mats @@ -655,7 +815,7 @@ class PMC(BaseLikelihood, InferenceMethod): Population Monte Carlo based inference scheme of Cappé et. al. [1]. This algorithm assumes a likelihood function is available and can be evaluated - at any parameter value given the oberved dataset. In absence of the + at any parameter value given the observed dataset. In absence of the likelihood function or when it can't be evaluated with a rational computational expenses, we use the approximated likelihood functions in abcpy.approx_lhd module, for which the argument of the consistency of the @@ -669,14 +829,15 @@ class PMC(BaseLikelihood, InferenceMethod): Parameters ---------- - model : list + root_models : list A list of the Probabilistic models corresponding to the observed datasets - loglikfun : abcpy.approx_lhd.Approx_likelihood - Approx_loglikelihood object defining the approximated loglikelihood to be used. - kernel : abcpy.distributions.Distribution - Distribution object defining the perturbation kernel needed for the sampling. + loglikfuns : list of abcpy.approx_lhd.Approx_likelihood + List of Approx_loglikelihood object defining the approximated loglikelihood to be used; one for each model. backend : abcpy.backends.Backend Backend object defining the backend to be used. + kernel : abcpy.perturbationkernel.PerturbationKernel, optional + PerturbationKernel object defining the perturbation kernel needed for the sampling. If not provided, the + DefaultKernel is used. seed : integer, optional Optional initial seed for the random number generator. The default value is generated randomly. @@ -717,23 +878,23 @@ def __init__(self, root_models, loglikfuns, backend, kernel=None, seed=None): self.simulation_counter = 0 def sample(self, observations, steps, n_samples=10000, n_samples_per_param=100, covFactors=None, iniPoints=None, - full_output=0, journal_file=None): + full_output=0, journal_file=None, path_to_save_journal=None): """Samples from the posterior distribution of the model parameter given the observed data observations. Parameters ---------- observations : list - A list, containing lists describing the observed data sets + A list, containing lists describing the observed data sets; one for each model. steps : integer - number of iterations in the sequential algoritm ("generations") + number of iterations in the sequential algorithm ("generations") n_samples : integer, optional number of samples to generate. The default value is 10000. n_samples_per_param : integer, optional number of data points in each simulated data set. The default value is 100. covFactors : list of float, optional scaling parameter of the covariance matrix. The default is a p dimensional array of 1 when p is the - dimension of the parameter. + dimension of the parameter. One for each perturbation kernel. inipoints : numpy.ndarray, optional parameter vaulues from where the sampling starts. By default sampled from the prior. full_output: integer, optional @@ -742,7 +903,11 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=100, journal_file: str, optional Filename of a journal file to read an already saved journal file, from which the first iteration will start. The default value is None. - + path_to_save_journal: str, optional + If provided, save the journal at the provided path. The journal is saved (and overwritten) after each step + of the sequential inference routine, so that partial results are stored to the disk in case the + inference routine does not end correctly; recall that you need to set full_output=1 to obtain the + full partial results. Returns ------- @@ -755,16 +920,19 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=100, self.n_samples = n_samples self.n_samples_per_param = n_samples_per_param + if path_to_save_journal is not None: + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else path_to_save_journal + '.jnl' + if journal_file is None: journal = Journal(full_output) journal.configuration["type_model"] = [type(model).__name__ for model in self.model] - journal.configuration["type_lhd_func"] = type(self.likfun).__name__ + journal.configuration["type_lhd_func"] = [type(likfun).__name__ for likfun in self.likfun.approx_lhds] + journal.configuration["steps"] = steps journal.configuration["n_samples"] = self.n_samples journal.configuration["n_samples_per_param"] = self.n_samples_per_param - journal.configuration["steps"] = steps journal.configuration["covFactor"] = covFactors journal.configuration["iniPoints"] = iniPoints - + journal.configuration["full_output"] = full_output else: journal = Journal.fromFile(journal_file) @@ -941,6 +1109,9 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=100, journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) + return journal # define helper functions for map step @@ -1038,9 +1209,9 @@ def _compute_accepted_cov_mats(self, covFactors, new_cov_mats): for covFactor, new_cov_mat in zip(covFactors, new_cov_mats): if not (new_cov_mat.size == 1): accepted_cov_mats.append( - covFactor * new_cov_mat + 0.0001 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) + covFactor * new_cov_mat + 1e-20 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) else: - accepted_cov_mats.append((covFactor * new_cov_mat + 0.0001 * new_cov_mat).reshape(1, 1)) + accepted_cov_mats.append((covFactor * new_cov_mat + 1e-20 * new_cov_mat).reshape(1, 1)) return accepted_cov_mats @@ -1054,14 +1225,16 @@ class SABC(BaseDiscrepancy, InferenceMethod): Parameters ---------- - model : list + root_models: list A list of the Probabilistic models corresponding to the observed datasets - distance : abcpy.distances.Distance - Distance object defining the distance measure used to compare simulated and observed data sets. - kernel : abcpy.distributions.Distribution - Distribution object defining the perturbation kernel needed for the sampling. + distances: list of abcpy.distances.Distance + List of Distance objects defining the distance measure to compare simulated and observed data sets; one for + each model. backend : abcpy.backends.Backend Backend object defining the backend to be used. + kernel : abcpy.perturbationkernel.PerturbationKernel, optional + PerturbationKernel object defining the perturbation kernel needed for the sampling. If not provided, the + DefaultKernel is used. seed : integer, optional Optional initial seed for the random number generator. The default value is generated randomly. """ @@ -1085,6 +1258,12 @@ def __init__(self, root_models, distances, backend, kernel=None, seed=None): # We define the joint Linear combination distance using all the distances for each individual models self.distance = LinearCombination(root_models, distances) + # check if the distance estimators are always positive + if np.any([isinstance(distance, Divergence) and not distance._estimate_always_positive() + for distance in distances]): + raise RuntimeError("SABC does not work with estimates of divergences which may be negative. Use another " + "inference algorithm or a different estimator.") + if kernel is None: mapping, garbage_index = self._get_mapping() @@ -1107,16 +1286,17 @@ def __init__(self, root_models, distances, backend, kernel=None, seed=None): self.simulation_counter = 0 def sample(self, observations, steps, epsilon, n_samples=10000, n_samples_per_param=1, beta=2, delta=0.2, - v=0.3, ar_cutoff=0.1, resample=None, n_update=None, full_output=0, journal_file=None): + v=0.3, ar_cutoff=0.1, resample=None, n_update=None, full_output=0, journal_file=None, + path_to_save_journal=None): """Samples from the posterior distribution of the model parameter given the observed data observations. Parameters ---------- observations : list - A list, containing lists describing the observed data sets + A list, containing lists describing the observed data sets; one for each model. steps : integer - Number of maximum iterations in the sequential algoritm ("generations") + Number of maximum iterations in the sequential algorithm ("generations") epsilon : numpy.float A proposed value of threshold to start with. n_samples : integer, optional @@ -1143,6 +1323,11 @@ def sample(self, observations, steps, epsilon, n_samples=10000, n_samples_per_pa journal_file: str, optional Filename of a journal file to read an already saved journal file, from which the first iteration will start. The default value is None. + path_to_save_journal: str, optional + If provided, save the journal at the provided path. The journal is saved (and overwritten) after each step + of the sequential inference routine, so that partial results are stored to the disk in case the + inference routine does not end correctly; recall that you need to set full_output=1 to obtain the + full partial results. Returns ------- @@ -1156,11 +1341,17 @@ def sample(self, observations, steps, epsilon, n_samples=10000, n_samples_per_pa self.n_samples = n_samples self.n_samples_per_param = n_samples_per_param + if path_to_save_journal is not None: + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else path_to_save_journal + '.jnl' + if journal_file is None: journal = Journal(full_output) journal.configuration["type_model"] = [type(model).__name__ for model in self.model] - journal.configuration["type_dist_func"] = type(self.distance).__name__ - journal.configuration["type_kernel_func"] = type(self.kernel) + journal.configuration["type_dist_func"] = [type(distance).__name__ for distance in self.distance.distances] + journal.configuration["type_kernel_func"] = [type(kernel).__name__ for kernel in self.kernel.kernels] if \ + isinstance(self.kernel, JointPerturbationKernel) else type(self.kernel) + journal.configuration["steps"] = steps + journal.configuration["epsilon"] = self.epsilon journal.configuration["n_samples"] = self.n_samples journal.configuration["n_samples_per_param"] = self.n_samples_per_param journal.configuration["beta"] = beta @@ -1373,21 +1564,24 @@ def sample(self, observations, steps, epsilon, n_samples=10000, n_samples_per_pa journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) - # Add epsilon_arr, number of final steps and final output to the journal - # print("INFO: Saving final configuration to output journal.") - if (full_output == 0) or (full_output == 1 and broken_preemptively and aStep <= steps - 1): - journal.add_accepted_parameters(copy.deepcopy(accepted_parameters)) - journal.add_weights(copy.deepcopy(accepted_weights)) - journal.add_ESS_estimate(accepted_weights) - journal.add_distances(copy.deepcopy(distances)) - self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, - accepted_weights=accepted_weights) - names_and_parameters = self._get_names_and_parameters() - journal.add_user_parameters(names_and_parameters) - journal.number_of_simulations.append(self.simulation_counter) + # Add epsilon_arr, number of final steps and final output to the journal + # print("INFO: Saving final configuration to output journal.") + if (full_output == 0) or (full_output == 1 and broken_preemptively and aStep <= steps - 1): + journal.add_accepted_parameters(copy.deepcopy(accepted_parameters)) + journal.add_weights(copy.deepcopy(accepted_weights)) + journal.add_ESS_estimate(accepted_weights) + journal.add_distances(copy.deepcopy(distances)) + self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, + accepted_weights=accepted_weights) + names_and_parameters = self._get_names_and_parameters() + journal.add_user_parameters(names_and_parameters) + journal.number_of_simulations.append(self.simulation_counter) - journal.configuration["steps"] = aStep + 1 - journal.configuration["epsilon"] = epsilon + journal.configuration["steps"] = aStep + 1 + journal.configuration["epsilon"] = epsilon + + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) return journal @@ -1558,9 +1752,9 @@ def _compute_accepted_cov_mats(self, beta, new_cov_mats): for new_cov_mat in new_cov_mats: if not (new_cov_mat.size == 1): accepted_cov_mats.append( - beta * new_cov_mat + 0.0001 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) + beta * new_cov_mat + 1e-20 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) else: - accepted_cov_mats.append((beta * new_cov_mat + 0.0001 * new_cov_mat).reshape(1, 1)) + accepted_cov_mats.append((beta * new_cov_mat + 1e-20 * new_cov_mat).reshape(1, 1)) return accepted_cov_mats @@ -1574,12 +1768,14 @@ class ABCsubsim(BaseDiscrepancy, InferenceMethod): ---------- model : list A list of the Probabilistic models corresponding to the observed datasets - distance : abcpy.distances.Distance - Distance object defining the distance used to compare the simulated and observed data sets. - kernel : abcpy.distributions.Distribution - Distribution object defining the perturbation kernel needed for the sampling. + distances: list of abcpy.distances.Distance + List of Distance objects defining the distance measure to compare simulated and observed data sets; one for + each model. backend : abcpy.backends.Backend Backend object defining the backend to be used. + kernel : abcpy.perturbationkernel.PerturbationKernel, optional + PerturbationKernel object defining the perturbation kernel needed for the sampling. If not provided, the + DefaultKernel is used. seed : integer, optional Optional initial seed for the random number generator. The default value is generated randomly. """ @@ -1623,16 +1819,16 @@ def __init__(self, root_models, distances, backend, kernel=None, seed=None): self.simulation_counter = 0 def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, chain_length=10, ap_change_cutoff=10, - full_output=0, journal_file=None): + full_output=0, journal_file=None, path_to_save_journal=None): """Samples from the posterior distribution of the model parameter given the observed data observations. Parameters ---------- observations : list - A list, containing lists describing the observed data sets + A list, containing lists describing the observed data sets; one for each model. steps : integer - Number of iterations in the sequential algoritm ("generations") + Number of iterations in the sequential algorithm ("generations") n_samples : integer, optional Number of samples to generate. The default value is 10000. n_samples_per_param : integer, optional @@ -1648,6 +1844,11 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ch journal_file: str, optional Filename of a journal file to read an already saved journal file, from which the first iteration will start. The default value is None. + path_to_save_journal: str, optional + If provided, save the journal at the provided path. The journal is saved (and overwritten) after each step + of the sequential inference routine, so that partial results are stored to the disk in case the + inference routine does not end correctly; recall that you need to set full_output=1 to obtain the + full partial results. Returns ------- @@ -1660,12 +1861,16 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ch self.chain_length = chain_length self.n_samples = n_samples self.n_samples_per_param = n_samples_per_param + if path_to_save_journal is not None: + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else path_to_save_journal + '.jnl' if journal_file is None: journal = Journal(full_output) journal.configuration["type_model"] = [type(model).__name__ for model in self.model] - journal.configuration["type_dist_func"] = type(self.distance).__name__ - journal.configuration["type_kernel_func"] = type(self.kernel) + journal.configuration["type_dist_func"] = [type(distance).__name__ for distance in self.distance.distances] + journal.configuration["type_kernel_func"] = [type(kernel).__name__ for kernel in self.kernel.kernels] if \ + isinstance(self.kernel, JointPerturbationKernel) else type(self.kernel) + journal.configuration["steps"] = steps journal.configuration["n_samples"] = self.n_samples journal.configuration["n_samples_per_param"] = self.n_samples_per_param journal.configuration["chain_length"] = self.chain_length @@ -1800,6 +2005,8 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ch self.logger.info(msg) if anneal_parameter_change_percentage < ap_change_cutoff: break + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) # Add anneal_parameter, number of final steps and final output to the journal # print("INFO: Saving final configuration to output journal.") @@ -1818,6 +2025,9 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ch journal.configuration["steps"] = aStep + 1 journal.configuration["anneal_parameter"] = anneal_parameter + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) + return journal # define helper functions for map step @@ -1971,12 +2181,14 @@ class RSMCABC(BaseDiscrepancy, InferenceMethod): ---------- model : list A list of the Probabilistic models corresponding to the observed datasets - distance : abcpy.distances.Distance - Distance object defining the distance measure used to compare simulated and observed data sets. - kernel : abcpy.distributions.Distribution - Distribution object defining the perturbation kernel needed for the sampling. + distances: list of abcpy.distances.Distance + List of Distance objects defining the distance measure to compare simulated and observed data sets; one for + each model. backend : abcpy.backends.Backend Backend object defining the backend to be used. + kernel : abcpy.perturbationkernel.PerturbationKernel, optional + PerturbationKernel object defining the perturbation kernel needed for the sampling. If not provided, the + DefaultKernel is used. seed : integer, optional Optional initial seed for the random number generator. The default value is generated randomly. """ @@ -2025,7 +2237,8 @@ def __init__(self, root_models, distances, backend, kernel=None, seed=None): self.simulation_counter = 0 def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, alpha=0.1, epsilon_init=100, - epsilon_final=0.1, const=0.01, covFactor=2.0, full_output=0, journal_file=None): + epsilon_final=0.1, const=0.01, covFactor=2.0, full_output=0, journal_file=None, + path_to_save_journal=None): """ Samples from the posterior distribution of the model parameter given the observed data observations. @@ -2033,9 +2246,9 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, al Parameters ---------- observations : list - A list, containing lists describing the observed data sets + A list, containing lists describing the observed data sets; one for each model. steps : integer - Number of iterations in the sequential algoritm ("generations") + Number of iterations in the sequential algorithm ("generations") n_samples : integer, optional Number of samples to generate. The default value is 10000. n_samples_per_param : integer, optional @@ -2056,6 +2269,11 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, al journal_file: str, optional Filename of a journal file to read an already saved journal file, from which the first iteration will start. The default value is None. + path_to_save_journal: str, optional + If provided, save the journal at the provided path. The journal is saved (and overwritten) after each step + of the sequential inference routine, so that partial results are stored to the disk in case the + inference routine does not end correctly; recall that you need to set full_output=1 to obtain the + full partial results. Returns ------- @@ -2069,14 +2287,24 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, al self.alpha = alpha self.n_samples = n_samples self.n_samples_per_param = n_samples_per_param + if path_to_save_journal is not None: + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else path_to_save_journal + '.jnl' if journal_file is None: journal = Journal(full_output) journal.configuration["type_model"] = [type(model).__name__ for model in self.model] - journal.configuration["type_dist_func"] = type(self.distance).__name__ + journal.configuration["type_dist_func"] = [type(distance).__name__ for distance in self.distance.distances] + journal.configuration["type_kernel_func"] = [type(kernel).__name__ for kernel in self.kernel.kernels] if \ + isinstance(self.kernel, JointPerturbationKernel) else type(self.kernel) + journal.configuration["steps"] = steps journal.configuration["n_samples"] = self.n_samples journal.configuration["n_samples_per_param"] = self.n_samples_per_param - journal.configuration["steps"] = steps + journal.configuration["alpha"] = alpha + journal.configuration["epsilon_init"] = epsilon_init + journal.configuration["epsilon_final"] = epsilon_final + journal.configuration["const"] = const + journal.configuration["covFactor"] = covFactor + journal.configuration["full_output"] = full_output else: journal = Journal.fromFile(journal_file) @@ -2217,8 +2445,11 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, al self.accepted_parameters_manager.update_broadcast(self.backend, accepted_weights=accepted_weights, accepted_parameters=accepted_parameters) - # Add epsilon_arr to the journal - journal.configuration["epsilon_arr"] = epsilon + # Add epsilon_arr to the journal + journal.configuration["epsilon_arr"] = epsilon + + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) return journal @@ -2313,9 +2544,9 @@ def _compute_accepted_cov_mats(self, covFactor, new_cov_mats): for new_cov_mat in new_cov_mats: if not (new_cov_mat.size == 1): accepted_cov_mats.append( - covFactor * new_cov_mat + 0.0001 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) + covFactor * new_cov_mat + 1e-20 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) else: - accepted_cov_mats.append((covFactor * new_cov_mat + 0.0001 * new_cov_mat).reshape(1, 1)) + accepted_cov_mats.append((covFactor * new_cov_mat + 1e-20 * new_cov_mat).reshape(1, 1)) return accepted_cov_mats @@ -2330,12 +2561,14 @@ class APMCABC(BaseDiscrepancy, InferenceMethod): ---------- model : list A list of the Probabilistic models corresponding to the observed datasets - distance : abcpy.distances.Distance - Distance object defining the distance measure used to compare simulated and observed data sets. - kernel : abcpy.distributions.Distribution - Distribution object defining the perturbation kernel needed for the sampling. + distances: list of abcpy.distances.Distance + List of Distance objects defining the distance measure to compare simulated and observed data sets; one for + each model. backend : abcpy.backends.Backend Backend object defining the backend to be used. + kernel : abcpy.perturbationkernel.PerturbationKernel, optional + PerturbationKernel object defining the perturbation kernel needed for the sampling. If not provided, the + DefaultKernel is used. seed : integer, optional Optional initial seed for the random number generator. The default value is generated randomly. """ @@ -2383,16 +2616,16 @@ def __init__(self, root_models, distances, backend, kernel=None, seed=None): self.simulation_counter = 0 def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, alpha=0.1, acceptance_cutoff=0.03, - covFactor=2.0, full_output=0, journal_file=None): + covFactor=2.0, full_output=0, journal_file=None, path_to_save_journal=None): """Samples from the posterior distribution of the model parameter given the observed data observations. Parameters ---------- observations : list - A list, containing lists describing the observed data sets + A list, containing lists describing the observed data sets; one for each model. steps : integer - Number of iterations in the sequential algoritm ("generations") + Number of iterations in the sequential algorithm ("generations") n_samples : integer, optional Number of samples to generate. The default value is 10000. n_samples_per_param : integer, optional @@ -2409,6 +2642,11 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, al journal_file: str, optional Filename of a journal file to read an already saved journal file, from which the first iteration will start. The default value is None. + path_to_save_journal: str, optional + If provided, save the journal at the provided path. The journal is saved (and overwritten) after each step + of the sequential inference routine, so that partial results are stored to the disk in case the + inference routine does not end correctly; recall that you need to set full_output=1 to obtain the + full partial results. Returns ------- @@ -2421,14 +2659,22 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, al self.alpha = alpha self.n_samples = n_samples self.n_samples_per_param = n_samples_per_param + if path_to_save_journal is not None: + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else path_to_save_journal + '.jnl' if journal_file is None: journal = Journal(full_output) journal.configuration["type_model"] = [type(model).__name__ for model in self.model] - journal.configuration["type_dist_func"] = type(self.distance).__name__ + journal.configuration["type_dist_func"] = [type(distance).__name__ for distance in self.distance.distances] + journal.configuration["type_kernel_func"] = [type(kernel).__name__ for kernel in self.kernel.kernels] if \ + isinstance(self.kernel, JointPerturbationKernel) else type(self.kernel) + journal.configuration["steps"] = steps journal.configuration["n_samples"] = self.n_samples journal.configuration["n_samples_per_param"] = self.n_samples_per_param - journal.configuration["steps"] = steps + journal.configuration["alpha"] = self.alpha + journal.configuration["acceptance_cutoff"] = acceptance_cutoff + journal.configuration["covFactor"] = covFactor + journal.configuration["full_output"] = full_output else: journal = Journal.fromFile(journal_file) @@ -2557,8 +2803,11 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, al if prob_acceptance < acceptance_cutoff: break - # Add epsilon_arr to the journal - journal.configuration["epsilon_arr"] = epsilon + # Add epsilon_arr to the journal + journal.configuration["epsilon_arr"] = epsilon + + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) return journal @@ -2647,32 +2896,46 @@ def _compute_accepted_cov_mats(self, covFactor, new_cov_mats): for new_cov_mat in new_cov_mats: if not (new_cov_mat.size == 1): accepted_cov_mats.append( - covFactor * new_cov_mat + 0.0001 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) + covFactor * new_cov_mat + 1e-20 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) else: - accepted_cov_mats.append((covFactor * new_cov_mat + 0.0001 * new_cov_mat).reshape(1, 1)) + accepted_cov_mats.append((covFactor * new_cov_mat + 1e-20 * new_cov_mat).reshape(1, 1)) return accepted_cov_mats class SMCABC(BaseDiscrepancy, InferenceMethod): - """This class implements Sequential Monte Carlo Approximate Bayesian computation of - Del Moral et al. [1]. + """This class implements two versions of Sequential Monte Carlo Approximate Bayesian computation, following either + the original in Del Moral et al. [1] or the newer version in Bernton et al. [3]. The first one is commonly used + when standard statistics based ABC is done (for instance with Euclidean distance), while the second one is instead + used when divergence-ABC is done (for instance, ABC with the Wasserstein distance, with the KL divergence and so + on). The first implementation in fact does not work in that case, as it assumes the distance can be computed with + respect to one single observation. Additionally, our implementation of the algorithm by Bernton et al. + does not work with the standard MCMC kernel, but requires using the r-hit kernel [2], which is arguably more + efficient (even if it leads to an algorithm for which the number of simulations is not known a priori). [1] P. Del Moral, A. Doucet, A. Jasra, An adaptive sequential Monte Carlo method for approximate Bayesian computation. Statistics and Computing, 22(5):1009–1020, 2012. - [2] Lee, Anthony. "n the choice of MCMC kernels for approximate Bayesian computation with SMC samplers. + [2] Lee, Anthony. "On the choice of MCMC kernels for approximate Bayesian computation with SMC samplers. Proceedings of the 2012 Winter Simulation Conference (WSC). IEEE, 2012. + [3] Bernton, E., Jacob, P. E., Gerber, M., & Robert, C. P. (2019). Approximate Bayesian computation with the + Wasserstein distance. Journal of the Royal Statistical Society Series B, 81(2), 235-269. + Parameters ---------- model : list A list of the Probabilistic models corresponding to the observed datasets - distance : abcpy.distances.Distance - Distance object defining the distance measure used to compare simulated and observed data sets. - kernel : abcpy.distributions.Distribution - Distribution object defining the perturbation kernel needed for the sampling. + distances: list of abcpy.distances.Distance + List of Distance objects defining the distance measure to compare simulated and observed data sets; one for + each model. backend : abcpy.backends.Backend Backend object defining the backend to be used. + kernel : abcpy.perturbationkernel.PerturbationKernel, optional + PerturbationKernel object defining the perturbation kernel needed for the sampling. If not provided, the + DefaultKernel is used. + version : string, optional + Denotes which version to use, either the one by Del Moral et al. [1] ("DelMoral") or the newer version in + Bernton et al. [3] ("Bernton"). By default, uses "DelMoral". seed : integer, optional Optional initial seed for the random number generator. The default value is generated randomly. """ @@ -2691,7 +2954,7 @@ class SMCABC(BaseDiscrepancy, InferenceMethod): backend = None - def __init__(self, root_models, distances, backend, kernel=None, seed=None): + def __init__(self, root_models, distances, backend, kernel=None, version="DelMoral", seed=None): self.model = root_models # We define the joint Linear combination distance using all the distances for each individual models self.distance = LinearCombination(root_models, distances) @@ -2706,6 +2969,22 @@ def __init__(self, root_models, distances, backend, kernel=None, seed=None): self.kernel = kernel self.backend = backend + if version not in ["DelMoral", "Bernton"]: + raise RuntimeError('The only implemented SMCABC methods are the one by Del Moral et al. ("DelMoral") or' + ' the one by Bernton et al. [3] ("Bernton"). Please set version' + ' to one of these two.') + else: + self.bernton = version == "Bernton" + self.version = version + + # check now if we are using divergences and in that case return error if using DelMoral version as that does + # not work: + if not self.bernton and np.any([isinstance(distance, Divergence) for distance in distances]): + raise RuntimeError("You requested to use the SMCABC algorithm by Del Moral et al. " + "together with divergences " + "between empirical measures. That algorithm however only works for standard statistics " + "based ABC.") + self.logger = logging.getLogger(__name__) self.epsilon = None @@ -2718,17 +2997,18 @@ def __init__(self, root_models, distances, backend, kernel=None, seed=None): self.simulation_counter = 0 - def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, epsilon_final=0.1, alpha=0.95, - covFactor=2, resample=None, full_output=0, which_mcmc_kernel=0, journal_file=None): + def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, epsilon_final=0.1, alpha=None, + covFactor=2, resample=None, full_output=0, which_mcmc_kernel=None, r=None, journal_file=None, + path_to_save_journal=None): """Samples from the posterior distribution of the model parameter given the observed data observations. Parameters ---------- observations : list - A list, containing lists describing the observed data sets + A list, containing lists describing the observed data sets; one for each model. steps : integer - Number of iterations in the sequential algoritm ("generations") + Number of iterations in the sequential algorithm ("generations") n_samples : integer, optional Number of samples to generate. The default value is 10000. n_samples_per_param : integer, optional @@ -2737,43 +3017,74 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ep The final threshold value of epsilon to be reached; if at some iteration you reach a lower epsilon than epsilon_final, the algorithm will stop and not proceed with further iterations. The default value is 0.1. alpha : float, optional - A parameter taking values between [0,1], determinining the rate of change of the threshold epsilon. The - default value is 0.95. + A parameter taking values between [0,1], determining the rate of change of the threshold epsilon. If + the algorithm by Bernton et al. is used, epsilon is chosen in order to guarantee a + proportion of unique particles equal to alpha + after resampling. If instead the algorithm by Del Moral et al. is used, + epsilon is chosen such that the ESS (Effective Sample Size) with + the new threshold value is alpha times the ESS with the old threshold value. The default value is None, + in which case 0.95 is used for Del Moral et al. algorithm, or 0.5 for Bernton et al. algorithm. covFactor : float, optional scaling parameter of the covariance matrix. The default value is 2. resample : float, optional It defines the resample step: introduce a resample step, after the particles have been - perturbed and the new weights have been computed, if the effective sample size is smaller than resample. If - not provided, resample is set to 0.5 * n_samples. + perturbed and the new weights have been computed, if the effective sample size is smaller than resample. + Notice that the algorithm by Bernton et al. always uses resample (as the weight values in that setup can + only be equal to 0 or 1), so that this parameter is ignored in that case. + If not provided, resample is set to 0.5 * n_samples. full_output: integer, optional If full_output==1, intermediate results are included in output journal. The default value is 0, meaning the intermediate results are not saved. which_mcmc_kernel: integer, optional - Specifies which MCMC kernel to be used: '0' kernel suggested in [1], any other value will use r-hit kernel - suggested by Anthony Lee [2]. The default value is 0. + Specifies which MCMC kernel to be used: '0' is the standard MCMC kernel used in the algorithm by + Del Moral et al [1], '1' uses the first version of r-hit kernel suggested by Anthony Lee (Alg. 5 in [2]), + while '2' uses the second version of r-hit kernel suggested by Anthony Lee (Alg. 6 in [2]). + The default value is 2 if the used algorithm is the one by Bernton et al, and is 0 for the algorithm by + Del Moral et al. + r: integer, optional: + Specifies the value of 'r' (the number of wanted hits) in the r-hits kernels. It is therefore ignored if + 'which_mcmc_kernel==0'. If no value is provided, the first version of r-hit kernel uses r=3, while the + second uses r=2. The default value is None. journal_file: str, optional Filename of a journal file to read an already saved journal file, from which the first iteration will start. The default value is None. + path_to_save_journal: str, optional + If provided, save the journal at the provided path. The journal is saved (and overwritten) after each step + of the sequential inference routine, so that partial results are stored to the disk in case the + inference routine does not end correctly; recall that you need to set full_output=1 to obtain the + full partial results. Returns ------- abcpy.output.Journal A journal containing simulation results, metadata and optionally intermediate results. """ - self.sample_from_prior(rng=self.rng) self.accepted_parameters_manager.broadcast(self.backend, observations) self.n_samples = n_samples self.n_samples_per_param = n_samples_per_param + if path_to_save_journal is not None: + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else path_to_save_journal + '.jnl' if journal_file is None: journal = Journal(full_output) journal.configuration["type_model"] = [type(model).__name__ for model in self.model] - journal.configuration["type_dist_func"] = type(self.distance).__name__ + journal.configuration["type_dist_func"] = [type(distance).__name__ for distance in self.distance.distances] + journal.configuration["type_kernel_func"] = [type(kernel).__name__ for kernel in self.kernel.kernels] if \ + isinstance(self.kernel, JointPerturbationKernel) else type(self.kernel) + journal.configuration["steps"] = steps journal.configuration["n_samples"] = self.n_samples journal.configuration["n_samples_per_param"] = self.n_samples_per_param - journal.configuration["steps"] = steps - # maybe add which kernel I am using? + journal.configuration["epsilon_final"] = epsilon_final + journal.configuration["alpha"] = alpha + journal.configuration["covFactor"] = covFactor + journal.configuration["resample"] = resample + journal.configuration["which_mcmc_kernel"] = which_mcmc_kernel + journal.configuration["r"] = r + journal.configuration["full_output"] = full_output + journal.configuration["version"] = self.version + self.sample_from_prior(rng=self.rng) # initialize only if you are not restarting from a journal, in order + # to ensure reproducibility else: journal = Journal.fromFile(journal_file) @@ -2786,8 +3097,29 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ep if resample is None: resample = n_samples * 0.5 - # Define epsilon_init - epsilon = [10000] + if alpha is None: + alpha = 0.5 if self.bernton else 0.95 + + # Define maximum value of epsilon + if not np.isinf(self.distance.dist_max()): + epsilon = [self.distance.dist_max()] + else: + epsilon = [1e5] + + if which_mcmc_kernel is None: + which_mcmc_kernel = 2 if self.bernton else 0 + + if which_mcmc_kernel not in [0, 1, 2]: + raise NotImplementedError("'which_mcmc_kernel' was given wrong value. It specifies which MCMC kernel to be" + " used: '0' kernel suggested in [1], '1' uses the first version of r-hit" + "kernel suggested by Anthony Lee (Alg. 5 in [2]), while '2' uses the second " + "version of r-hit kernel" + "suggested by Anthony Lee (Alg. 6 in [2]). The default value is 0.") + + if self.bernton and which_mcmc_kernel == 0: + raise RuntimeError("The algorithm by Bernton et al. does not work with the standard MCMC kernel.") + + self.r = r # main SMC ABC algorithm for aStep in range(0, steps): @@ -2797,6 +3129,9 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ep accepted_parameters = journal.get_accepted_parameters(-1) accepted_weights = journal.get_weights(-1) accepted_y_sim = journal.opt_values[-1] + distances = journal.get_distances(-1) + + epsilon = journal.configuration["epsilon_arr"] self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, accepted_weights=accepted_weights) @@ -2820,39 +3155,57 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ep # 0: Compute the Epsilon if accepted_y_sim != None: - self.logger.info("Compute epsilon, might take a while") - # Compute epsilon for next step - fun = lambda epsilon_var: self._compute_epsilon(epsilon_var, \ - epsilon, observations, accepted_y_sim, accepted_weights, - n_samples, n_samples_per_param, alpha) - epsilon_new = self._bisection(fun, epsilon_final, epsilon[-1], 0.001) + self.logger.info( + "Compute epsilon, might take a while; previous epsilon value: {:.4f}".format(epsilon[-1])) + if self.bernton: + # for the Bernton algorithm, the distances have already been computed before during the acceptance + # step. This however holds only when the r-hit kernels are used + + # first compute the distances for the current set of parameters and observations + # (notice that may have already been done somewhere before!): + # current_distance_matrix = self._compute_distance_matrix_divergence(observations, accepted_y_sim, + # n_samples) + # assert np.allclose(current_distance_matrix, distances) + current_distance_matrix = distances + + # Compute epsilon for next step + fun = self._def_compute_epsilon_divergence_unique_particles(n_samples, + current_distance_matrix, alpha) + epsilon_new = self._bisection(fun, epsilon_final, epsilon[-1], 0.001) + else: + # first compute the distances for the current set of parameters and observations + # (notice that may have already been done somewhere before!): + current_distance_matrix = self._compute_distance_matrix(observations, accepted_y_sim, n_samples, + n_samples_per_param) + fun = self._def_compute_epsilon(epsilon, accepted_weights, n_samples, current_distance_matrix, + alpha) + epsilon_new = self._bisection(fun, epsilon_final, epsilon[-1], 0.001) if epsilon_new < epsilon_final: epsilon_new = epsilon_final epsilon.append(epsilon_new) # 1: calculate weights for new parameters self.logger.info("Calculating weights") - if accepted_y_sim != None: - new_weights = np.zeros(shape=n_samples, ) - for ind1 in range(n_samples): - numerator = 0.0 - denominator = 0.0 - for ind2 in range(n_samples_per_param): - numerator += (self.distance.distance(observations, [[accepted_y_sim[ind1][0][ind2]]]) < epsilon[ - -1]) - denominator += ( - self.distance.distance(observations, [[accepted_y_sim[ind1][0][ind2]]]) < epsilon[-2]) - if denominator != 0.0: - new_weights[ind1] = accepted_weights[ind1] * (numerator / denominator) - else: - new_weights[ind1] = 0 + if accepted_y_sim is not None: + if self.bernton: + new_weights = (current_distance_matrix < epsilon[-1]) * 1 + else: + numerators = np.sum(current_distance_matrix < epsilon[-1], axis=1) + denominators = np.sum(current_distance_matrix < epsilon[-2], axis=1) + + non_zero_denominator = denominators != 0 + new_weights = np.zeros(shape=n_samples) + + new_weights[non_zero_denominator] = accepted_weights.flatten()[non_zero_denominator] * ( + numerators[non_zero_denominator] / denominators[non_zero_denominator]) new_weights = new_weights / sum(new_weights) else: new_weights = np.ones(shape=n_samples, ) * (1.0 / n_samples) - - # 2: Resample - if accepted_y_sim is not None and pow(sum(pow(new_weights, 2)), -1) < resample: + # 2: Resample; we resample always when using the Bernton et al. algorithm, as in that case weights + # can only be proportional to 1 or 0; if we use the Del Moral version, instead, the + # weights can have fractional values -> use the # resample threshold + if accepted_y_sim is not None and (self.bernton or pow(sum(pow(new_weights, 2)), -1) < resample): self.logger.info("Resampling") # Weighted resampling: index_resampled = self.rng.choice(n_samples, n_samples, replace=True, p=new_weights) @@ -2897,11 +3250,13 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ep self._update_broadcasts(accepted_y_sim) # calculate resample parameters - self.logger.info("Drawing perturbed sampless") + self.logger.info("Drawing perturbed samples") if which_mcmc_kernel == 0: params_and_ysim_pds = self.backend.map(self._accept_parameter, rng_and_index_pds) - else: + elif which_mcmc_kernel == 1: params_and_ysim_pds = self.backend.map(self._accept_parameter_r_hit_kernel, rng_and_index_pds) + elif which_mcmc_kernel == 2: + params_and_ysim_pds = self.backend.map(self._accept_parameter_r_hit_kernel_version_2, rng_and_index_pds) params_and_ysim = self.backend.collect(params_and_ysim_pds) new_parameters, new_y_sim, distances, counter = [list(t) for t in zip(*params_and_ysim)] distances = np.array(distances) @@ -2926,57 +3281,140 @@ def sample(self, observations, steps, n_samples=10000, n_samples_per_param=1, ep journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) - # Add epsilon_arr to the journal - journal.configuration["epsilon_arr"] = epsilon + # Add epsilon_arr to the journal + journal.configuration["epsilon_arr"] = epsilon + + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) return journal - def _compute_epsilon(self, epsilon_new, epsilon, observations, accepted_y_sim, accepted_weights, n_samples, - n_samples_per_param, alpha): + def _compute_distance_matrix(self, observations, accepted_y_sim, n_samples, n_samples_per_param): + distance_matrix = np.zeros((n_samples, n_samples_per_param)) + + for ind1 in range(n_samples): + for ind2 in range(n_samples_per_param): + distance_matrix[ind1, ind2] = self.distance.distance(observations, [[accepted_y_sim[ind1][0][ind2]]]) + self.logger.debug( + 'Computed distance inside for weights:' + str( + distance_matrix[ind1, ind2])) + return distance_matrix + + def _compute_distance_matrix_divergence(self, observations, accepted_y_sim, n_samples): + distance_matrix = np.zeros(n_samples) + for ind1 in range(n_samples): + distance_matrix[ind1] = self.distance.distance(observations, [accepted_y_sim[ind1][0]]) + self.logger.debug('Computed distance matrix for weights:' + str(distance_matrix[ind1])) + return distance_matrix + + @staticmethod + def _def_compute_epsilon(epsilon, accepted_weights, n_samples, distance_matrix, alpha): """ + Returns a function of 'epsilon_new' that is used in the bisection routine. + The distances are computed just once in the definition of the function; this therefore avoids computing them + repeatedly during the _bisection routine for the same input values, which is inefficient. + Parameters ---------- - epsilon_new: float - New value for epsilon. epsilon: float Current threshold. - observations: numpy.ndarray - Observed data. - accepted_y_sim: numpy.ndarray - Accepted simulated data. accepted_weights: numpy.ndarray Accepted weights. n_samples: integer Number of samples to generate. - n_samples_per_param: integer - Number of data points in each simulated data set. + distance_matrix: np.ndarray: + the distance matrix between observation and data used to compute the weights alpha: float Returns ------- - float - Newly computed value for threshold. + callable + The function used in the bisection routine """ RHS = alpha * pow(sum(pow(accepted_weights, 2)), -1) - LHS = np.zeros(shape=n_samples, ) - for ind1 in range(n_samples): - numerator = 0.0 - denominator = 0.0 - for ind2 in range(n_samples_per_param): - numerator += (self.distance.distance(observations, [[accepted_y_sim[ind1][0][ind2]]]) < epsilon_new) - denominator += (self.distance.distance(observations, [[accepted_y_sim[ind1][0][ind2]]]) < epsilon[-1]) - if denominator == 0: - LHS[ind1] = 0 + denominators = np.sum(distance_matrix < epsilon[-1], axis=1) + non_zero_denominator = denominators != 0 + + def _compute_epsilon(epsilon_new): + """ + Parameters + ---------- + epsilon_new: float + New value for epsilon. + + Returns + ------- + float + Newly computed value for threshold. + """ + # old version (not optimized): + # for ind1 in range(n_samples): + # numerator = 0.0 + # denominator = 0.0 + # for ind2 in range(n_samples_per_param): + # numerator += (distance_matrix[ind1, ind2] < epsilon_new) + # denominator += (distance_matrix[ind1, ind2] < epsilon[-1]) + # if denominator == 0: + # LHS[ind1] = 0 + # else: + # LHS[ind1] = accepted_weights[ind1] * (numerator / denominator) + + numerators = np.sum(distance_matrix < epsilon_new, axis=1) + + LHS = np.zeros(shape=n_samples) + + LHS[non_zero_denominator] = accepted_weights.flatten()[non_zero_denominator] * ( + numerators[non_zero_denominator] / denominators[non_zero_denominator]) + if sum(LHS) == 0: + result = RHS else: - LHS[ind1] = accepted_weights[ind1] * (numerator / denominator) - if sum(LHS) == 0: - result = RHS - else: - LHS = LHS / sum(LHS) - LHS = pow(sum(pow(LHS, 2)), -1) - result = RHS - LHS - return result + LHS = LHS / sum(LHS) # normalize weights. + LHS = pow(sum(pow(LHS, 2)), -1) + result = RHS - LHS + + return result + + return _compute_epsilon + + def _def_compute_epsilon_divergence_unique_particles(self, n_samples, distance_matrix, alpha): + """ + Parameters + ---------- + n_samples: integer + Number of samples to generate. + alpha: float + + Returns + ------- + callable + The function used in the bisection routine + """ + + def _compute_epsilon_divergence_unique_particles(epsilon_new): + """ + Parameters + ---------- + epsilon_new: float + New value for epsilon. + Returns + ------- + float + proportion of unique particles after resampling + """ + new_weights = (distance_matrix < epsilon_new) * 1 + self.logger.debug('New weights:' + str(new_weights)) + if sum(new_weights) != 0: + new_weights = new_weights / sum(new_weights) + rng = np.random.RandomState(1) # this fixes the randomness across iterations; it makes sense therefore + # Here we want a proportion of unique particles equal to alpha after resampling + result = (len( + np.unique(rng.choice(n_samples, n_samples, replace=True, p=new_weights))) / n_samples) - alpha + else: + result = - alpha + return result + + return _compute_epsilon_divergence_unique_particles def _bisection(self, func, low, high, tol): # cache computed values, as we call func below @@ -3054,6 +3492,8 @@ def _accept_parameter(self, rng_and_index, npc=None): for ind in range(self.n_samples_per_param): numerator += (self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), [[y_sim[0][ind]]]) < self.epsilon[-1]) + # we have most likely already computed this distance before, but hard to keep track. Moreover this + # is parallelized -> should not impact too much on computing time denominator += (self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), [[y_sim_old[0][ind]]]) < self.epsilon[-1]) if denominator == 0: @@ -3084,7 +3524,7 @@ def _accept_parameter(self, rng_and_index, npc=None): def _accept_parameter_r_hit_kernel(self, rng_and_index, npc=None): """ This implements algorithm 5 in Lee (2012) [2] which is used as an MCMC kernel in SMCABC. This implementation - uses r=3. + uses r=3 as default value. Parameters ---------- @@ -3105,7 +3545,7 @@ def _accept_parameter_r_hit_kernel(self, rng_and_index, npc=None): rng.seed(rng.randint(np.iinfo(np.uint32).max, dtype=np.uint32)) # Set value of r for r-hit kernel - r = 3 + r = 3 if self.r is None else self.r mapping_for_kernels, garbage_index = self.accepted_parameters_manager.get_mapping( self.accepted_parameters_manager.model) @@ -3114,35 +3554,41 @@ def _accept_parameter_r_hit_kernel(self, rng_and_index, npc=None): if self.accepted_parameters_manager.accepted_parameters_bds is None: self.sample_from_prior(rng=rng) y_sim = self.simulate(self.n_samples_per_param, rng=rng, npc=npc) + # the following is was probably already computed before, but hard to keep track: + distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) counter += 1 else: if self.accepted_parameters_manager.accepted_weights_bds.value()[index] > 0: theta = self.accepted_parameters_manager.accepted_parameters_bds.value()[index] - # Sample from theta until we get 'r-1' y_sim inside the epsilon ball + # Sample from theta until we get 'r-1' y_sim inside the epsilon ball (line 4 in Alg 5 in [3]) self.set_parameters(theta) - accept_old_arr, y_sim_old_arr, N_old = [], [], 0 - while sum(accept_old_arr) < r - 1: + accept_old_arr, N_old = [], 0 + # y_sim_old_arr = [] this is actually not used. + while len(accept_old_arr) < r - 1: y_sim = self.simulate(self.n_samples_per_param, rng=rng, npc=npc) - y_sim_old_arr.append(y_sim) + # y_sim_old_arr.append(y_sim) if self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) < self.epsilon[-1]: accept_old_arr.append(N_old) N_old += 1 counter += 1 - # Perturb and sample from the perturbed theta until we get 'r' y_sim inside the epsilon ball + # Perturb and sample from the perturbed theta until we get 'r' y_sim inside the epsilon ball + # (line 2 in Alg 5 in [3]) while True: perturbation_output = self.perturb(index, rng=rng) if perturbation_output[0] and self.pdf_of_prior(self.model, perturbation_output[1]) != 0: break - accept_new_arr, y_sim_new_arr, N = [], [], 0 - while sum(accept_new_arr) < r: + accept_new_arr, y_sim_new_arr, distance_new_arr, N = [], [], [], 0 + while len(accept_new_arr) < r: y_sim = self.simulate(self.n_samples_per_param, rng=rng, npc=npc) y_sim_new_arr.append(y_sim) - if self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), - y_sim) < self.epsilon[-1]: + distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), + y_sim) + if distance < self.epsilon[-1]: accept_new_arr.append(N) + distance_new_arr.append(distance) counter += 1 N += 1 @@ -3159,16 +3605,156 @@ def _accept_parameter_r_hit_kernel(self, rng_and_index, npc=None): if rng.binomial(1, acceptance_prob) == 1: self.set_parameters(perturbation_output[1]) - # Randomly sample index J - J = rng.choice(accept_new_arr).astype(int) + # Randomly sample index J between the first r-1 hits + J = rng.choice(accept_new_arr[:-1]).astype(int) y_sim = y_sim_new_arr[J] + distance = distance_new_arr[J] else: self.set_parameters(theta) y_sim = self.accepted_y_sim_bds.value()[index] + # the following is was probably already computed before, but hard to keep track: + distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) else: self.set_parameters(self.accepted_parameters_manager.accepted_parameters_bds.value()[index]) y_sim = self.accepted_y_sim_bds.value()[index] - distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) + # the following distance was probably already computed before? + distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) + return self.get_parameters(), y_sim, distance, counter + + def _accept_parameter_r_hit_kernel_version_2(self, rng_and_index, npc=None): + """ + This implements algorithm 6 in Lee (2012) [2] which is used as an MCMC kernel in SMCABC. This implementation + uses r=2 as default value. + + Parameters + ---------- + rng_and_index: numpy.ndarray + 2 dimensional array. The first entry is a random number generator. + The second entry defines the index in the data set. + + Returns + ------- + Tuple + The first entry of the tuple is the accepted parameters. The second entry is the simulated data set. + The third one is the distance between the simulated data set and the observation, while the fourth one is + the number of simulations needed to obtain the accepted parameter. + """ + + rng = rng_and_index[0] + index = rng_and_index[1] + rng.seed(rng.randint(np.iinfo(np.uint32).max, dtype=np.uint32)) + + # Set value of r for r-hit kernel + r = 2 if self.r is None else self.r + mapping_for_kernels, garbage_index = self.accepted_parameters_manager.get_mapping( + self.accepted_parameters_manager.model) + + counter = 0 + # print("on seed " + str(seed) + " distance: " + str(distance) + " epsilon: " + str(self.epsilon)) + if self.accepted_parameters_manager.accepted_parameters_bds is None: + self.sample_from_prior(rng=rng) + y_sim = self.simulate(self.n_samples_per_param, rng=rng, npc=npc) + distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) + # the following is was probably already computed before, but hard to keep track: + counter += 1 + else: + if self.accepted_parameters_manager.accepted_weights_bds.value()[index] > 0: + theta = self.accepted_parameters_manager.accepted_parameters_bds.value()[index] + + # Generate different perturbed values from theta and sample from model until we get 'r' y_sim + # inside the epsilon ball (line 1 in Alg 6 in [3]) + accept_prime_arr, z_sim_arr, theta_prime_arr, distance_arr, N_prime = [], [], [], [], 0 + while len(accept_prime_arr) < r: + # first perturb: + while True: + # this perturbs using the current theta value + perturbation_output = self.perturb(index, rng=rng) + if perturbation_output[0] and self.pdf_of_prior(self.model, perturbation_output[1]) != 0: + break + theta_prime_arr.append(perturbation_output[1]) + # now simulate + z_sim = self.simulate(self.n_samples_per_param, rng=rng, npc=npc) + z_sim_arr.append(z_sim) # could store here only the accepted ones... Can improve + distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), z_sim) + if distance < self.epsilon[-1]: + accept_prime_arr.append(N_prime) + distance_arr.append(distance) + N_prime += 1 + counter += 1 + # select the index among the first r-1 hits (line 2 in Alg 6) + L = rng.choice(accept_prime_arr[:-1]).astype(int) + theta_prime_L = theta_prime_arr[L] + + # create a new AcceptedParametersManager storing the theta_prime_L in order to draw perturbations from + # it: + inner_accepted_parameters_manager = AcceptedParametersManager(self.model) + # define a dummy backend: + backend_inner = BackendDummy() + # need to pass the covariance matrix (take it from the overall AcceptedParametersManager) and the + # theta_prime_L: + inner_accepted_parameters_manager.update_broadcast(backend_inner, accepted_parameters=[theta_prime_L], + accepted_cov_mats=self.accepted_parameters_manager.accepted_cov_mats_bds.value()) + + # in kernel_parameters you need to store theta_prime_L + # inner_accepted_parameters_manager.update_kernel_values(backend_inner, + # kernel_parameters=[[theta_prime_L]]) + + # update the kernel parameters (which is used to perturb - this is more general than the one above which + # works with one kernel only) + kernel_parameters = [] + for kernel in self.kernel.kernels: + kernel_parameters.append( + inner_accepted_parameters_manager.get_accepted_parameters_bds_values(kernel.models)) + + inner_accepted_parameters_manager.update_kernel_values(backend_inner, + kernel_parameters=kernel_parameters) + + # in kernel_parameters you need to store theta_prime_L + # inner_accepted_parameters_manager.update_kernel_values(backend_inner, + # kernel_parameters=[[theta_prime_L]]) + + # Generate different perturbed values from the parameter value selected above and sample from model + # until we get 'r-1' y_sim inside the epsilon ball (line 3 in Alg 6 in [3]) + accept_arr, N = [], 0 + while len(accept_arr) < r - 1: + while True: + perturbation_output = self.perturb(0, rng=rng, + accepted_parameters_manager=inner_accepted_parameters_manager) + if perturbation_output[0] and self.pdf_of_prior(self.model, perturbation_output[1]) != 0: + break + x_sim = self.simulate(self.n_samples_per_param, rng=rng, npc=npc) # + # y_sim_new_arr.append(y_sim) + if self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), + x_sim) < self.epsilon[-1]: + accept_arr.append(N) + counter += 1 + N += 1 + + # Calculate acceptance probability (line 4 in Alg 6) + ratio_prior_prob = self.pdf_of_prior(self.model, theta_prime_L) / self.pdf_of_prior(self.model, + theta) + kernel_numerator = self.kernel.pdf(mapping_for_kernels, self.accepted_parameters_manager, + theta_prime_L, theta) + kernel_denominator = self.kernel.pdf(mapping_for_kernels, self.accepted_parameters_manager, theta, + theta_prime_L) + ratio_likelihood_prob = kernel_numerator / kernel_denominator + + acceptance_prob = min(1, (N / (N_prime - 1)) * ratio_prior_prob * ratio_likelihood_prob) + + if rng.binomial(1, acceptance_prob) == 1: + self.set_parameters(theta_prime_L) + y_sim = z_sim_arr[L] + distance = distance_arr[L] + else: + self.set_parameters(theta) + y_sim = self.accepted_y_sim_bds.value()[index] + # the following is was probably already computed before, but hard to keep track: + distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) + else: + self.set_parameters(self.accepted_parameters_manager.accepted_parameters_bds.value()[index]) + y_sim = self.accepted_y_sim_bds.value()[index] + # the following is was probably already computed before, but hard to keep track: + distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) return self.get_parameters(), y_sim, distance, counter def _compute_accepted_cov_mats(self, covFactor, new_cov_mats): @@ -3192,7 +3778,538 @@ def _compute_accepted_cov_mats(self, covFactor, new_cov_mats): for new_cov_mat in new_cov_mats: if not (new_cov_mat.size == 1): accepted_cov_mats.append( - covFactor * new_cov_mat + 0.0001 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) + covFactor * new_cov_mat + 1e-20 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) + else: + accepted_cov_mats.append((covFactor * new_cov_mat + 1e-20 * new_cov_mat).reshape(1, 1)) + return accepted_cov_mats + + +class MCMCMetropoliHastings(BaseLikelihood, InferenceMethod): + """ + Simple Metropolis-Hastings MCMC working with the approximate likelihood functions Approx_likelihood, with + multivariate normal proposals. + + Parameters + ---------- + root_models : list + A list of the Probabilistic models corresponding to the observed datasets + loglikfuns : list of abcpy.approx_lhd.Approx_likelihood + List of Approx_loglikelihood object defining the approximated loglikelihood to be used; one for each model. + backend : abcpy.backends.Backend + Backend object defining the backend to be used. + kernel : abcpy.perturbationkernel.PerturbationKernel, optional + PerturbationKernel object defining the perturbation kernel needed for the sampling. If not provided, the + DefaultKernel is used. + seed : integer, optional + Optional initial seed for the random number generator. The default value is generated randomly. + + """ + + model = None + likfun = None + kernel = None + rng = None + + n_samples = None + n_samples_per_param = None + + backend = None + + def __init__(self, root_models, loglikfuns, backend, kernel=None, seed=None): + self.model = root_models + # We define the joint Sum of Loglikelihood functions using all the loglikelihoods for each individual models + self.likfun = SumCombination(root_models, loglikfuns) + + mapping, garbage_index = self._get_mapping() + models = [] + self.parameter_names_with_index = {} + for mdl, mdl_index in mapping: + models.append(mdl) + self.parameter_names_with_index[mdl.name] = mdl_index # dict storing param names with index + + self.parameter_names = [model.name for model in models] # store parameter names + + if kernel is None: + kernel = DefaultKernel(models) + + self.kernel = kernel + self.backend = backend + self.rng = np.random.RandomState(seed) + self.logger = logging.getLogger(__name__) + + # these are usually big tables, so we broadcast them to have them once + # per executor instead of once per task + self.accepted_parameters_manager = AcceptedParametersManager(self.model) + # this is used to handle the data for adapting the covariance: + self.accepted_parameters_manager_adaptive_cov = AcceptedParametersManager(self.model) + + self.simulation_counter = 0 + + def sample(self, observations, n_samples, n_samples_per_param=100, burnin=1000, cov_matrices=None, iniPoint=None, + adapt_proposal_cov_interval=None, covFactor=None, bounds=None, speedup_dummy=True, use_tqdm=True, + journal_file=None, path_to_save_journal=None): + """Samples from the posterior distribution of the model parameter given the observed + data observations. The MCMC is run for burnin + n_samples steps, and n_samples_per_param are used at each step + to estimate the approximate loglikelihood. The burnin steps are then discarded from the chain stored in the + journal file. + + During burnin, the covariance matrix is adapted from the steps generated up to that point, in a way similar to + what suggested in [1], after each adapt_proposal_cov_interval steps. Differently from the original algorithm in + [1], here the proposal covariance matrix is fixed after the end of the burnin steps. + + In case the original parameter space is bounded (for instance with uniform prior on an interval), the MCMC can + be optionally run on a transformed space. Therefore, the covariance matrix describes proposals on the + transformed space; the acceptance rate then takes into account the Jacobian of the transformation. In order to + use MCMC with transformed space, you need to specify lower and upper bounds in the corresponding parameters (see + details in the description of `bounds`). + + The returned journal file contains also information on acceptance rates (in the configuration dictionary). + + [1] Haario, H., Saksman, E., & Tamminen, J. (2001). An adaptive Metropolis algorithm. Bernoulli, 7(2), 223-242. + + Parameters + ---------- + observations : list + A list, containing lists describing the observed data sets; one for each model. + n_samples : integer, optional + number of samples to generate. The default value is 10000. + n_samples_per_param : integer, optional + number of data points in each simulated data set. The default value is 100. + burnin : integer, optional + Number of burnin steps to discard. Defaults to 1000. + cov_matrices : list of matrices, optional + list of initial covariance matrices for the proposals. If not provided, identity matrices are used. If the + sample routine is restarting from a journal file and cov_matrices is not provided, cov_matrices is set to + the value used for sampling after burnin in the previous journal file (ie what is stored in + `journal.configuration["actual_cov_matrices"]`). + iniPoint : numpy.ndarray, optional + parameter value from where the sampling starts. By default sampled from the prior. Not used if journal_file + is passed. + adapt_proposal_cov_interval : integer, optional + the proposal covariance matrix is adapted each adapt_cov_matrix steps during burnin, by using the chain up + to that point. If None, no adaptation is done. Default value is None. Use with care as, if the likelihood + estimate is very noisy, the adaptation may work pooly (see `covFactor` parameter). + covFactor : float, optional + the factor by which to scale the empirical covariance matrix in order to obtain the covariance matrix for + the proposal, whenever that is updated during the burnin steps. If not provided, we use the default value + 2.4 ** 2 / dim_theta suggested in [1]. + Notice that this value was shown to be optimal (at least in some + limit sense) in the case in which the likelihood is fully known. In the present case, in which the + likelihood is estimated from data, that value may turn out to be too large; specifically, if + the likelihood estimate is very noisy, that choice may lead to a very bad adaptation which may give rise + to an MCMC which does not explore the space well (for instance, the obtained covariance matrix may turn out + to be too small). If that happens, we suggest to set covFactor to a smaller value than the default one, in + which case the acceptance rate of the chain will likely be smaller but the exploration will be better. + Alternatively, it is possible to reduce the noise in the likelihood estimate by increasing + `n_samples_per_param`. + bounds : dictionary, optional + dictionary containing the lower and upper bound for the transformation to be applied to the parameters. The + key of each entry is the name of the parameter as defined in the model, while the value if a tuple (or list) + with `(lower_bound, upper_bound)` content. If the parameter is bounded on one side only, the other bound + should be set to 'None'. If a parameter is not in this dictionary, no transformation is applied to it. + If a parameter is bounded on two sides, the used transformation is based on the logit. If conversely it is + lower bounded, we apply instead a log transformation. Notice that we do not implement yet the transformation + for upper bounded variables. If no value is provided, the default value is None, which means no + transformation at all is applied. + speedup_dummy: boolean, optional. + If set to True, the map function is not used to parallelize simulations (for the new parameter value) when + the backend is Dummy. This can improve performance as it can exploit potential vectorization in the model. + However, this breaks reproducibility when using, for instance, BackendMPI with respect to BackendDummy, due + to the different way the random seeds are used when speedup_dummy is set to True. Please set this to False + if you are interested in preserving reproducibility across MPI and Dummy backend. Defaults to True. + use_tqdm : boolean, optional + Whether using tqdm or not to display progress. Defaults to True. + journal_file: str, optional + Filename of a journal file to read an already saved journal file, from which the first iteration will start. + That's the only information used (it does not use the previous covariance matrix). + The default value is None. + path_to_save_journal: str, optional + If provided, save the journal at the provided path at the end of the inference routine. + + Returns + ------- + abcpy.output.Journal + A journal containing simulation results, metadata and optionally intermediate results. + """ + + self.observations = observations + self.n_samples = n_samples + self.n_samples_per_param = n_samples_per_param + self.speedup_dummy = speedup_dummy + # we use this in all places which require a backend but which are not parallelized in MCMC: + self.dummy_backend = BackendDummy() + dim = len(self.parameter_names) + + if path_to_save_journal is not None: + path_to_save_journal = path_to_save_journal if '.jnl' in path_to_save_journal else path_to_save_journal + '.jnl' + + if bounds is None: + # no transformation is performed + self.transformer = DummyTransformer() + else: + if not isinstance(bounds, dict): + raise TypeError("Argument `bounds` need to be a dictionary") + bounds_keys = bounds.keys() + for key in bounds_keys: + if key not in self.parameter_names: + raise KeyError("The keys in argument `bounds` need to correspond to the parameter names used " + "in defining the model") + if not hasattr(bounds[key], "__len__") or len(bounds[key]) != 2: + raise RuntimeError("Each entry in `bounds` need to be a tuple with 2 value, representing the lower " + "and upper bound of the corresponding parameter. If the parameter is bounded on " + "one side only, the other bound should be set to 'None'.") + + # create lower_bounds and upper_bounds_vector: + lower_bound_transformer = np.array([None] * dim) + upper_bound_transformer = np.array([None] * dim) + + for key in bounds_keys: + lower_bound_transformer[self.parameter_names_with_index[key]] = bounds[key][0] + upper_bound_transformer[self.parameter_names_with_index[key]] = bounds[key][1] + + # initialize transformer: + self.transformer = BoundedVarTransformer(np.array(lower_bound_transformer), + np.array(upper_bound_transformer)) + + accepted_parameters = [] + accepted_parameters_burnin = [] + if journal_file is None: + journal = Journal(0) + journal.configuration["type_model"] = [type(model).__name__ for model in self.model] + journal.configuration["type_lhd_func"] = [type(likfun).__name__ for likfun in self.likfun.approx_lhds] + journal.configuration["type_kernel_func"] = [type(kernel).__name__ for kernel in self.kernel.kernels] if \ + isinstance(self.kernel, JointPerturbationKernel) else type(self.kernel) + journal.configuration["n_samples"] = self.n_samples + journal.configuration["n_samples_per_param"] = self.n_samples_per_param + journal.configuration["burnin"] = burnin + journal.configuration["cov_matrices"] = cov_matrices + journal.configuration["iniPoint"] = iniPoint + journal.configuration["adapt_proposal_cov_interval"] = adapt_proposal_cov_interval + journal.configuration["covFactor"] = covFactor + journal.configuration["bounds"] = bounds + journal.configuration["speedup_dummy"] = speedup_dummy + journal.configuration["use_tqdm"] = use_tqdm + journal.configuration["acceptance_rates"] = [] + # Initialize chain: when not supplied, randomly draw it from prior distribution + # It is an MCMC chain: weights are always 1; forget about them + # accepted_parameter will keep track of the chain position + if iniPoint is None: + self.sample_from_prior(rng=self.rng) + accepted_parameter = self.get_parameters() + else: + accepted_parameter = iniPoint + if burnin == 0: + accepted_parameters_burnin.append(accepted_parameter) + self.logger.info("Calculate approximate loglikelihood") + approx_log_likelihood_accepted_parameter = self._simulate_and_compute_log_lik(accepted_parameter) + # update the number of simulations (this tracks the number of parameters for which simulations are done; + # the actual number of simulations is this times n_samples_per_param)) + self.simulation_counter += 1 + self.acceptance_rate = 0 + else: + # check the following: + self.logger.info("Restarting from previous journal") + journal = Journal.fromFile(journal_file) + # this is used to compute the overall acceptance rate: + self.acceptance_rate = journal.configuration["acceptance_rates"][-1] * journal.configuration["n_samples"] + accepted_parameter = journal.get_accepted_parameters(-1)[-1] # go on from last MCMC step + journal.configuration["n_samples"] += self.n_samples # add the total number of samples + journal.configuration["burnin"] = burnin + if journal.configuration["n_samples_per_param"] != self.n_samples_per_param: + warnings.warn("You specified a different n_samples_per_param from the one used in the passed " + "journal_file; the algorithm will still work fine.") + journal.configuration["n_samples_per_param"] = self.n_samples_per_param + journal.configuration["cov_matrices"] = cov_matrices + journal.configuration["bounds"] = bounds # overwrite + if cov_matrices is None: # use the previously stored one unless the user defines it + cov_matrices = journal.configuration["actual_cov_matrices"] + journal.configuration["speedup_dummy"] = speedup_dummy + approx_log_likelihood_accepted_parameter = journal.final_step_loglik + self.simulation_counter = journal.number_of_simulations[-1] # update the number of simulations + + if covFactor is None: + covFactor = 2.4 ** 2 / dim + + accepted_parameter_prior_pdf = self.pdf_of_prior(self.model, accepted_parameter) + + # set the accepted parameter in the kernel (in order to correctly generate next proposal) + # do this on transformed parameter + accepted_parameter_transformed = self.transformer.transform(accepted_parameter) + self._update_kernel_parameters(accepted_parameter_transformed) + + # compute jacobian + log_det_jac_accepted_param = self.transformer.jac_log_det(accepted_parameter) + + # 3: calculate covariance + self.logger.info("Set kernel covariance matrix ") + if cov_matrices is None: + # need to set that to some value (we use identity matrices). Be careful, there need to be one + # covariance matrix for each kernel; not sure that this works in case of multivariate parameters. + + # the kernel parameters are only used to get the exact shape of cov_matrices + cov_matrices = [np.eye(len(self.accepted_parameters_manager.kernel_parameters_bds.value()[0][kernel_index])) + for kernel_index in range(len(self.kernel.kernels))] + + self.accepted_parameters_manager.update_broadcast(self.dummy_backend, accepted_cov_mats=cov_matrices) + + # main MCMC algorithm + self.logger.info("Starting MCMC") + for aStep in tqdm(range(burnin + n_samples), disable=not use_tqdm): + + self.logger.debug("Step {} of MCMC algorithm".format(aStep)) + + # 1: Resample parameters + self.logger.debug("Generate proposal") + + # perturb element 0 of accepted_parameters_manager.kernel_parameters_bds: + # new_parameter = self.perturb(0, rng=self.rng)[1] # do not use this as it leads to some weird error. + # rather do: + new_parameters_transformed = self.kernel.update(self.accepted_parameters_manager, 0, rng=self.rng) + + self._reset_flags() # not sure whether this is needed, leave it anyway + + # Order the parameters provided by the kernel in depth-first search order + new_parameter_transformed = self.get_correct_ordering(new_parameters_transformed) + + # transform back + new_parameter = self.transformer.inverse_transform(new_parameter_transformed) + + # for now we are only using a simple MVN proposal. For bounded parameter values, this is not great; we + # could also implement a proposal on transformed space, which would be better. + new_parameter_prior_pdf = self.pdf_of_prior(self.model, new_parameter) + if new_parameter_prior_pdf == 0: + self.logger.debug("Proposal parameter at step {} is out of prior region.".format(aStep)) + if aStep >= burnin: + accepted_parameters.append(accepted_parameter) + else: + accepted_parameters_burnin.append(accepted_parameter) + continue + + # 2: calculate approximate likelihood for new parameter. If the backend is MPI, we distribute simulations + # and then compute the approx likelihood locally + self.logger.debug("Calculate approximate loglikelihood") + approx_log_likelihood_new_parameter = self._simulate_and_compute_log_lik(new_parameter) + self.simulation_counter += 1 # update the number of simulations + + log_det_jac_new_param = self.transformer.jac_log_det(new_parameter) + # log_det_jac_accepted_param = self.transformer.jac_log_det(accepted_parameter) + log_jac_term = log_det_jac_accepted_param - log_det_jac_new_param + + # compute acceptance rate: + alpha = np.exp( + log_jac_term + approx_log_likelihood_new_parameter - approx_log_likelihood_accepted_parameter) * ( + new_parameter_prior_pdf) / (accepted_parameter_prior_pdf) # assumes symmetric kernel + + # Metropolis-Hastings step: + if self.rng.uniform() < alpha: + # update param value and approx likelihood + accepted_parameter_transformed = new_parameter_transformed + accepted_parameter = new_parameter + approx_log_likelihood_accepted_parameter = approx_log_likelihood_new_parameter + accepted_parameter_prior_pdf = new_parameter_prior_pdf + log_det_jac_accepted_param = log_det_jac_new_param + # set the accepted parameter in the kernel (in order to correctly generate next proposal) + self._update_kernel_parameters(accepted_parameter_transformed) + if aStep >= burnin: + self.acceptance_rate += 1 + + # save to the trace: + if aStep >= burnin: + accepted_parameters.append(accepted_parameter) + else: + accepted_parameters_burnin.append(accepted_parameter) + + # adapt covariance of proposal: + if adapt_proposal_cov_interval is not None and (aStep + 1) % adapt_proposal_cov_interval == 0: + # store the accepted_parameters for adapting the covariance in the kernel. + # I use this piece of code as it formats the data in the right way + # for the sake of using them to compute the kernel cov: + self.accepted_parameters_manager_adaptive_cov.update_broadcast( + self.dummy_backend, accepted_parameters=accepted_parameters_burnin) + kernel_parameters = [] + for kernel in self.kernel.kernels: + kernel_parameters.append( + self.accepted_parameters_manager_adaptive_cov.get_accepted_parameters_bds_values( + kernel.models)) + self.accepted_parameters_manager_adaptive_cov.update_kernel_values( + self.dummy_backend, kernel_parameters=kernel_parameters) + + self.logger.info("Updating covariance matrix") + cov_matrices = self.kernel.calculate_cov(self.accepted_parameters_manager_adaptive_cov) + # this scales with the cov_Factor: + cov_matrices = self._compute_accepted_cov_mats(covFactor, cov_matrices) + # store it in the main AcceptedParametersManager in order to perturb data with it in the following: + self.accepted_parameters_manager.update_broadcast(self.dummy_backend, + accepted_cov_mats=cov_matrices) + + self.acceptance_rate /= journal.configuration["n_samples"] + self.logger.info("Saving results to output journal") + self.accepted_parameters_manager.update_broadcast(self.dummy_backend, accepted_parameters=accepted_parameters) + names_and_parameters = self._get_names_and_parameters() + if journal_file is not None: # concatenate chains + journal.add_accepted_parameters(journal.get_accepted_parameters() + copy.deepcopy(accepted_parameters)) + names_and_parameters = [(names_and_parameters[i][0], + journal.get_parameters()[names_and_parameters[i][0]] + names_and_parameters[i][1]) + for i in range(len(names_and_parameters))] + journal.add_user_parameters(names_and_parameters) + else: + journal.add_accepted_parameters(copy.deepcopy(accepted_parameters)) + journal.add_user_parameters(names_and_parameters) + journal.number_of_simulations.append(self.simulation_counter) + journal.configuration["acceptance_rates"].append(self.acceptance_rate) + journal.add_weights(np.ones((journal.configuration['n_samples'], 1))) + # store the final loglik to be able to restart the journal correctly + journal.final_step_loglik = approx_log_likelihood_accepted_parameter + # store the final actual cov_matrices, in order to use this when restarting from journal + journal.configuration["actual_cov_matrices"] = cov_matrices + + if path_to_save_journal is not None: # save journal + journal.save(path_to_save_journal) + + return journal + + def _sample_parameter(self, rng, npc=None): + """ + Generate a simulation from the model with the current value of accepted_parameter + + Parameters + ---------- + rng: random number generator + The random number generator to be used. + Returns + ------- + np.array + accepted parameter + """ + + # get the new parameter value + theta = self.new_parameter_bds.value() + # Simulate the fake data from the model given the parameter value theta + self.logger.debug("Simulate model for parameter " + str(theta)) + acc = self.set_parameters(theta) + if acc is False: + self.logger.debug("Parameter " + str(theta) + "has not been accepted") + y_sim = self.simulate(1, rng=rng, npc=npc) + + return y_sim + + def _approx_log_lik_calc(self, y_sim, npc=None): + """ + Compute likelihood for new parameters using approximate likelihood function + + Parameters + ---------- + y_sim: list + A list containing self.n_samples_per_param simulations for the new parameter value + Returns + ------- + float + The approximated likelihood function + """ + self.logger.debug("Extracting observation.") + obs = self.observations + + self.logger.debug("Computing likelihood...") + loglhd = self.likfun.loglikelihood(obs, y_sim) + + self.logger.debug("LogLikelihood is :" + str(loglhd)) + + return loglhd + + def _simulate_and_compute_log_lik(self, new_parameter): + """Helper function which simulates data from `new_parameter` and computes the approximate loglikelihood. + In case the backend is not BackendDummy (ie parallelization is available) this parallelizes the different + simulations (which are all for the same parameter value). + + Notice that, according to the used model, spreading the simulations in different tasks can be more inefficient + than using one single call, according to the level of vectorization that the model uses and the overhead + associated. For this reason, we do not split the simulations in different tasks when the backend is + BackendDummy. + + Parameters + ---------- + new_parameter + Parameter value from which to generate data with which to compute the approximate loglikelihood. + + Returns + ------- + float + The approximated likelihood function + """ + if isinstance(self.backend, BackendDummy) and self.speedup_dummy: + # do all the required simulations here without parallellizing; however this gives different result + # from the other option due to the way random seeds are handled. + self.logger.debug('simulations') + theta = new_parameter + # Simulate the fake data from the model given the parameter value theta + self.logger.debug("Simulate model for parameter " + str(theta)) + acc = self.set_parameters(theta) + if acc is False: + self.logger.debug("Parameter " + str(theta) + "has not been accepted") + simulations_from_new_parameter = self.simulate(n_samples_per_param=self.n_samples_per_param, rng=self.rng) + else: + self.logger.debug('parallelize simulations for fixed parameter value') + seed_arr = self.rng.randint(0, np.iinfo(np.uint32).max, size=self.n_samples_per_param, dtype=np.uint32) + rng_arr = np.array([np.random.RandomState(seed) for seed in seed_arr]) + rng_pds = self.backend.parallelize(rng_arr) + + # need first to broadcast the new_parameter value: + self.new_parameter_bds = self.backend.broadcast(new_parameter) + + # map step: + simulations_from_new_parameter_pds = self.backend.map(self._sample_parameter, rng_pds) + self.logger.debug("collect simulations from pds") + simulations_from_new_parameter = self.backend.collect(simulations_from_new_parameter_pds) + # now need to reshape that correctly. The first index has to be the model, then we need to have + # n_samples_per_param and then the size of the simulation + simulations_from_new_parameter = [ + [simulations_from_new_parameter[sample_index][model_index][0] for sample_index in + range(self.n_samples_per_param)] for model_index in range(len(self.model))] + approx_log_likelihood_new_parameter = self._approx_log_lik_calc(simulations_from_new_parameter) + + return approx_log_likelihood_new_parameter + + def _update_kernel_parameters(self, accepted_parameter): + """This stores the last accepted parameter in the kernel so that it will be used to generate the new proposal + with self.perturb. + + Parameters + ---------- + accepted_parameter + Parameter value from which you want to generate proposal at next iteration of MCMC. + """ + # I use this piece of code as it formats the data in the right way for the sake of using it in the kernel + self.accepted_parameters_manager.update_broadcast(self.dummy_backend, accepted_parameters=[accepted_parameter]) + + kernel_parameters = [] + for kernel in self.kernel.kernels: + kernel_parameters.append( + self.accepted_parameters_manager.get_accepted_parameters_bds_values(kernel.models)) + self.accepted_parameters_manager.update_kernel_values(self.dummy_backend, kernel_parameters=kernel_parameters) + + @staticmethod + def _compute_accepted_cov_mats(covFactor, new_cov_mats): + """ + Update the covariance matrices computed from data by multiplying them with covFactor and adding a small term in + the diagonal for numerical stability. + + Parameters + ---------- + covFactor : float + factor to correct the covariance matrices + new_cov_mats : list + list of covariance matrices computed from data + Returns + ------- + list + List of new accepted covariance matrices + """ + # accepted_cov_mats = [covFactor * cov_mat for cov_mat in accepted_cov_mats] + accepted_cov_mats = [] + for new_cov_mat in new_cov_mats: + if not (new_cov_mat.size == 1): + accepted_cov_mats.append( + covFactor * new_cov_mat + 1e-20 * np.trace(new_cov_mat) * np.eye(new_cov_mat.shape[0])) else: - accepted_cov_mats.append((covFactor * new_cov_mat + 0.0001 * new_cov_mat).reshape(1, 1)) + accepted_cov_mats.append((covFactor * new_cov_mat + 1e-20 * new_cov_mat).reshape(1, 1)) return accepted_cov_mats diff --git a/abcpy/output.py b/abcpy/output.py index 887a6880..13e05b71 100644 --- a/abcpy/output.py +++ b/abcpy/output.py @@ -386,7 +386,7 @@ def posterior_histogram(self, iteration=None, n_bins=10): H, edges = np.histogramdd(np.hstack(joined_params), bins=n_bins, weights=weights.reshape(len(weights), )) return [H, edges] - # TODO this does not work for multidimensional parameters + # TODO this does not work for multidimensional parameters and discrete distributions def plot_posterior_distr(self, parameters_to_show=None, ranges_parameters=None, iteration=None, show_samples=None, single_marginals_only=False, double_marginals_only=False, write_posterior_mean=True, show_posterior_mean=True, true_parameter_values=None, contour_levels=14, figsize=None, @@ -873,7 +873,58 @@ def Wass_convergence_plot(self, num_iter_max=1e8, **kwargs): ax.set_xlabel("Iteration") ax.set_ylabel("Wasserstein distance") ax.legend() - # put horizontal line at the largest value ESS can get: ax.set_xticks(np.arange(len(self.weights) - 1) + 1) return fig, ax, wass_dist_lists + + def traceplot(self, parameters_to_show=None, iteration=None, **kwargs): + """ + Produces a traceplot for the MCMC inference scheme. This only works for journal files which were created by the + `MCMCMetropolisHastings` inference scheme. + + + Parameters + ---------- + parameters_to_show : list, optional + a list of the parameters for which you want to plot the traceplot. For each parameter, you need + to provide the name string as it was defined in the model. For instance, + `jrnl.traceplot(parameters_to_show=["mu"])` will only plot the traceplot for the + parameter named "mu" in the list of parameters. + If `None`, then all parameters will be displayed. + iteration : int, optional + specify the iteration for which to plot the posterior distribution, in the case of a sequential algorithm. + If `None`, then the last iteration will be used. + kwargs + Additional arguments passed to matplotlib.pyplot.plot + + Returns + ------- + tuple + a tuple containing the matplotlib "fig, axes" objects defining the plot. Can be useful for further + modifications. + """ + + if not "acceptance_rates" in self.configuration.keys(): + raise RuntimeError("The traceplot can be produced only for journal files which were created by the" + " MCMCMetropolisHastings inference scheme") + + if parameters_to_show is None: + parameters_to_show = list(self.names_and_parameters[0].keys()) + else: + for param_name in parameters_to_show: + if param_name not in self.names_and_parameters[0].keys(): + raise KeyError("Parameter '{}' is not a parameter of the model.".format(param_name)) + + param_dict = self.get_parameters(iteration) + n_plots = len(parameters_to_show) + + fig, ax = plt.subplots(1, n_plots, figsize=(4 * n_plots, 4)) + fig.suptitle("Traceplot") + if n_plots == 1: + ax = [ax] + for i, name in enumerate(parameters_to_show): + ax[i].plot(np.squeeze(param_dict[name])) + ax[i].set_title(name) + ax[i].set_xlabel("MCMC step") + + return fig, ax diff --git a/abcpy/perturbationkernel.py b/abcpy/perturbationkernel.py index 796d077d..6c03cb56 100644 --- a/abcpy/perturbationkernel.py +++ b/abcpy/perturbationkernel.py @@ -63,7 +63,7 @@ def update(self, accepted_parameters_manager, row_index, rng): raise NotImplementedError - def pdf(self, accepted_parameters_manager, kernel_index, row_index, x): + def pdf(self, accepted_parameters_manager, kernel_index, mean, x): """ Calculates the pdf of the kernel at point x. @@ -73,8 +73,8 @@ def pdf(self, accepted_parameters_manager, kernel_index, row_index, x): The accepted parameters manager that manages all bds objects. kernel_index: integer The index of the kernel in the list of kernels of the joint perturbation kernel. - row_index: integer - The index of the accepted parameters bds for which the pdf should be evaluated. + mean: np array, np.float or np.integer + The reference point of the kernel x: list or float The point at which the pdf should be evaluated. @@ -86,7 +86,7 @@ def pdf(self, accepted_parameters_manager, kernel_index, row_index, x): """ if isinstance(self, DiscreteKernel): - return self.pmf(accepted_parameters_manager, kernel_index, row_index, x) + return self.pmf(accepted_parameters_manager, kernel_index, mean, x) else: raise NotImplementedError @@ -95,7 +95,7 @@ class ContinuousKernel(metaclass=ABCMeta): """This abstract base class represents all perturbation kernels acting on continuous parameters.""" @abstractmethod - def pdf(self, accepted_parameters_manager, kernel_index, index, x): + def pdf(self, accepted_parameters_manager, kernel_index, mean, x): raise NotImplementedError @@ -103,7 +103,7 @@ class DiscreteKernel(metaclass=ABCMeta): """This abstract base class represents all perturbation kernels acting on discrete parameters.""" @abstractmethod - def pmf(self, accepted_parameters_manager, kernel_index, index, x): + def pmf(self, accepted_parameters_manager, kernel_index, mean, x): raise NotImplementedError @@ -213,8 +213,8 @@ def pdf(self, mapping, accepted_parameters_manager, mean, x): index in the accepted_parameters_bds list corresponding to an output of this model. accepted_parameters_manager: abcpy.AcceptedParametersManager object The AcceptedParametersManager to be used. - index: integer - The row to be considered in the accepted_parameters_bds matrix. + mean: np array, np.float or np.integer + The reference point of the kernel x: The point at which the pdf should be evaluated. Returns @@ -243,7 +243,6 @@ class MultivariateNormalKernel(PerturbationKernel, ContinuousKernel): def __init__(self, models): self.models = models - def calculate_cov(self, accepted_parameters_manager, kernel_index): """ Calculates the covariance matrix relevant to this kernel. @@ -338,8 +337,8 @@ def pdf(self, accepted_parameters_manager, kernel_index, mean, x): The AcceptedParametersManager to be used. kernel_index: integer The index of the kernel in the list of kernels in the joint kernel. - index: integer - The row to be considered in the accepted_parameters_bds matrix. + mean: np array, np.float or np.integer + The reference point of the kernel x: The point at which the pdf should be evaluated. Returns @@ -483,8 +482,8 @@ def pdf(self, accepted_parameters_manager, kernel_index, mean, x): The AcceptedParametersManager to be used. kernel_index: integer The index of the kernel in the list of kernels in the joint kernel. - index: integer - The row to be considered in the accepted_parameters_bds matrix. + mean: np array, np.float or np.integer + The reference point of the kernel x: The point at which the pdf should be evaluated. Returns @@ -526,7 +525,7 @@ def pdf(self, accepted_parameters_manager, kernel_index, mean, x): class RandomWalkKernel(PerturbationKernel, DiscreteKernel): - def __init__(self, models): + def __init__(self, models, jump=1): """ This class defines a kernel perturbing discrete parameters using a naive random walk. @@ -537,15 +536,18 @@ def __init__(self, models): """ self.models = models + self.jump = jump def update(self, accepted_parameters_manager, kernel_index, row_index, rng=np.random.RandomState()): """ - Updates the parameter values contained in the accepted_paramters_manager using a random walk. + Updates the parameter values contained in the accepted_paramters_manager using a multivariate normal distribution. Parameters ---------- accepted_parameters_manager: abcpy.AcceptedParametersManager object Defines the AcceptedParametersManager to be used. + kernel_index: integer + The index of the kernel in the list of kernels in the joint kernel. row_index: integer The index of the row that should be considered from the accepted_parameters_bds matrix. rng: random number generator @@ -557,16 +559,37 @@ def update(self, accepted_parameters_manager, kernel_index, row_index, rng=np.ra The perturbed parameter values. """ + # Get all parameters relevant to this kernel + discrete_model_values = accepted_parameters_manager.kernel_parameters_bds.value()[kernel_index][row_index] - # Get parameter values relevant to this kernel - discrete_model_values = accepted_parameters_manager.kernel_parameters_bds.value()[kernel_index] - - perturbed_discrete_values = [] - discrete_model_values = np.array(discrete_model_values)[row_index] + if isinstance(discrete_model_values[0], + (np.float, np.float32, np.float64, np.int, np.int32, np.int64)): + # Perturb + discrete_model_values = np.array(discrete_model_values) + perturbed_discrete_values = [] + # Implement a random walk for the discrete parameter values + for discrete_value in discrete_model_values: + perturbed_discrete_values.append( + rng.randint(discrete_value - self.jump, discrete_value + self.jump + 1)) + perturbed_discrete_values = np.array(perturbed_discrete_values) + else: + # Learn the structure + struct = [[] for i in range(len(discrete_model_values))] + for i in range(len(discrete_model_values)): + struct[i] = discrete_model_values[i].shape[0] + struct = np.array(struct).cumsum() + discrete_model_values = np.concatenate(discrete_model_values) - # Implement a random walk for the discrete parameter values - for discrete_value in discrete_model_values: - perturbed_discrete_values.append(np.array([rng.randint(discrete_value - 1, discrete_value + 2)])) + # Perturb + discrete_model_values = np.array(discrete_model_values) + perturbed_discrete_values = [] + # Implement a random walk for the discrete parameter values + for discrete_value in discrete_model_values: + perturbed_discrete_values.append( + rng.randint(discrete_value - self.jump, discrete_value + self.jump + 1)) + perturbed_discrete_values = np.array(perturbed_discrete_values) + # Perturbed values anc split according to the structure + perturbed_discrete_values = np.split(perturbed_discrete_values, struct)[:-1] return perturbed_discrete_values @@ -590,8 +613,8 @@ def pmf(self, accepted_parameters_manager, kernel_index, mean, x): The AcceptedParametersManager to be used. kernel_index: integer The index of the kernel in the list of kernels of the joint kernel. - index: integer - The row to be considered in the accepted_parameters_bds matrix. + mean: integer + The reference point of the kernel x: The point at which the pdf should be evaluated. Returns @@ -599,8 +622,152 @@ def pmf(self, accepted_parameters_manager, kernel_index, mean, x): float The pmf evaluated at point x. """ + if np.abs(mean[0] - x[0]) > self.jump: + return 0 + else: + return 1. / (2 * self.jump + 1) + + +class NetworkRandomWalkKernel(PerturbationKernel, DiscreteKernel): + def __init__(self, models, network, name_weight): + """ + This class defines a kernel perturbing discrete parameters on a provided network with moves proportional to an attribute of the edegs of the network. + + Parameters + ---------- + models: list + List of abcpy.ProbabilisticModel objects + network: A network + Networkx object + name_weight: string + name of the attribute of the network to be used as probability + + """ + + self.models = models + self.network = network + self.name_weight = name_weight + + def update(self, accepted_parameters_manager, kernel_index, row_index, rng=np.random.RandomState()): + """ + Updates the parameter values contained in the accepted_paramters_manager. + + Parameters + ---------- + accepted_parameters_manager: abcpy.AcceptedParametersManager object + Defines the AcceptedParametersManager to be used. + kernel_index: integer + The index of the kernel in the list of kernels in the joint kernel. + row_index: integer + The index of the row that should be considered from the accepted_parameters_bds matrix. + rng: random number generator + The random number generator to be used. + + Returns + ------- + np.ndarray + The perturbed parameter values. + + """ + # Get all parameters relevant to this kernel + discrete_model_values = accepted_parameters_manager.kernel_parameters_bds.value()[kernel_index][row_index] + + if isinstance(discrete_model_values[0], + (np.float, np.float32, np.float64, np.int, np.int32, np.int64)): + # Perturb + discrete_model_values = np.array(discrete_model_values) + perturbed_discrete_values = [] + # Implement a random walk for the discrete parameter values + for discrete_value in discrete_model_values: + nodes_proposed = list(self.network.neighbors(discrete_value)) + weight = np.zeros(shape=(len(nodes_proposed),)) + for ind in range(len(nodes_proposed)): + weight[ind] = self.network[1][nodes_proposed[ind]][self.name_weight] + weight = weight / sum(weight) + perturbed_discrete_values.append(np.random.choice(nodes_proposed, 1, p=weight)[0]) + perturbed_discrete_values = np.array(perturbed_discrete_values) + else: + # Learn the structure + struct = [[] for i in range(len(discrete_model_values))] + for i in range(len(discrete_model_values)): + struct[i] = discrete_model_values[i].shape[0] + struct = np.array(struct).cumsum() + discrete_model_values = np.concatenate(discrete_model_values) + + # Perturb + discrete_model_values = np.array(discrete_model_values) + perturbed_discrete_values = [] + # Implement a random walk for the discrete parameter values + for discrete_value in discrete_model_values: + nodes_proposed = list(self.network.neighbors(discrete_value)) + weight = np.zeros(shape=(len(nodes_proposed),)) + for ind in range(len(nodes_proposed)): + weight[ind] = self.network[discrete_value][nodes_proposed[ind]][self.name_weight] + weight = weight / sum(weight) + perturbed_discrete_values.append(np.random.choice(nodes_proposed, 1, p=weight)[0]) + perturbed_discrete_values = np.array(perturbed_discrete_values) + # Perturbed values anc split according to the structure + perturbed_discrete_values = np.split(perturbed_discrete_values, struct)[:-1] + + return perturbed_discrete_values + + def calculate_cov(self, accepted_parameters_manager, kernel_index): + """ + Calculates the covariance matrix of this kernel. Since there is no covariance matrix associated with this + random walk, it returns an empty list. + """ - return 1. / 3 + return np.array([0]).reshape(-1, ) + + def pmf(self, accepted_parameters_manager, kernel_index, mean, x): + """Calculates the pdf of the kernel. + Commonly used to calculate weights. + + Parameters + ---------- + accepted_parameters_manager: abcpy.AcceptedParametersManager object + The AcceptedParametersManager to be used. + kernel_index: integer + The index of the kernel in the list of kernels in the joint kernel. + mean: np array, np.float or np.integer + The reference point of the kernel + x: The point at which the pdf should be evaluated. + + Returns + ------- + float + The pdf evaluated at point x. + """ + density = 1 + if isinstance(mean[0], + (np.float, np.float32, np.float64, np.int, np.int32, np.int64)): + mean = np.array(mean).astype(int) + for ind1 in range(len(mean)): + discrete_value = mean[ind1] + nodes_proposed = list(self.network.neighbors(discrete_value)) + weight = np.zeros(shape=(len(nodes_proposed),)) + for ind2 in range(len(nodes_proposed)): + weight[ind2] = self.network[discrete_value][nodes_proposed[ind2]][self.name_weight] + weight = weight / sum(weight) + if x[ind1] in nodes_proposed: + density = density * weight[np.where(nodes_proposed == x[ind1])[0]][0] + else: + density = density * 0 + return density + else: + mean = np.array(np.concatenate(mean)).astype(int) + for ind1 in range(len(mean)): + discrete_value = mean[ind1] + nodes_proposed = list(self.network.neighbors(discrete_value)) + weight = np.zeros(shape=(len(nodes_proposed),)) + for ind2 in range(len(nodes_proposed)): + weight[ind2] = self.network[discrete_value][nodes_proposed[ind2]][self.name_weight] + weight = weight / sum(weight) + if x[ind1] in nodes_proposed: + density = density * weight[np.where(nodes_proposed == x[ind1])[0]][0] + else: + density = density * 0 + return density class DefaultKernel(JointPerturbationKernel): diff --git a/abcpy/statisticslearning.py b/abcpy/statisticslearning.py index 36ef0914..f667253a 100644 --- a/abcpy/statisticslearning.py +++ b/abcpy/statisticslearning.py @@ -286,7 +286,7 @@ class StatisticsLearningNN(StatisticsLearning, GraphTools): def __init__(self, model, statistics_calc, backend, training_routine, distance_learning, embedding_net=None, n_samples=1000, n_samples_val=0, n_samples_per_param=1, parameters=None, simulations=None, parameters_val=None, simulations_val=None, seed=None, cuda=None, scale_samples=True, quantile=0.1, - **training_routine_kwargs): + use_tqdm=True, **training_routine_kwargs): """ Parameters ---------- @@ -355,6 +355,8 @@ def __init__(self, model, statistics_calc, backend, training_routine, distance_l Default value is True. quantile: float, optional quantile used to define the similarity set if distance_learning is True. Default to 0.1. + use_tqdm : boolean, optional + Whether using tqdm or not to display progress. Defaults to True. training_routine_kwargs: additional kwargs to be passed to the underlying training routine. """ @@ -439,12 +441,12 @@ def __init__(self, model, statistics_calc, backend, training_routine, distance_l if distance_learning: self.embedding_net = training_routine(simulations, similarity_set, embedding_net=self.embedding_net, - cuda=cuda, samples_val=simulations_val, + cuda=cuda, samples_val=simulations_val, use_tqdm=use_tqdm, similarity_set_val=similarity_set_val, **training_routine_kwargs) else: self.embedding_net = training_routine(simulations, target, embedding_net=self.embedding_net, cuda=cuda, samples_val=simulations_val, target_val=target_val, - **training_routine_kwargs) + use_tqdm=use_tqdm, **training_routine_kwargs) self.logger.info("Finished learning the transformation.") @@ -486,8 +488,8 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample n_samples_per_param=1, parameters=None, simulations=None, parameters_val=None, simulations_val=None, early_stopping=False, epochs_early_stopping_interval=1, start_epoch_early_stopping=10, seed=None, cuda=None, scale_samples=True, batch_size=16, n_epochs=200, load_all_data_GPU=False, - lr=1e-3, optimizer=None, - scheduler=None, start_epoch_training=0, optimizer_kwargs={}, scheduler_kwargs={}, loader_kwargs={}): + lr=1e-3, optimizer=None, scheduler=None, start_epoch_training=0, use_tqdm=True, + optimizer_kwargs={}, scheduler_kwargs={}, loader_kwargs={}): """ Parameters ---------- @@ -583,6 +585,8 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample the scheduler and the optimizer at each epoch. Default to 0. verbose: boolean, optional if True, prints more information from the training routine. Default to False. + use_tqdm : boolean, optional + Whether using tqdm or not to display progress. Defaults to True. optimizer_kwargs: Python dictionary, optional dictionary containing optional keyword arguments for the optimizer. scheduler_kwargs: Python dictionary, optional @@ -602,7 +606,7 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample seed=seed, cuda=cuda, scale_samples=scale_samples, batch_size=batch_size, n_epochs=n_epochs, load_all_data_GPU=load_all_data_GPU, lr=lr, optimizer=optimizer, scheduler=scheduler, - start_epoch_training=start_epoch_training, + start_epoch_training=start_epoch_training, use_tqdm=use_tqdm, optimizer_kwargs=optimizer_kwargs, scheduler_kwargs=scheduler_kwargs, loader_kwargs=loader_kwargs) @@ -625,7 +629,8 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample early_stopping=False, epochs_early_stopping_interval=1, start_epoch_early_stopping=10, seed=None, cuda=None, scale_samples=True, quantile=0.1, batch_size=16, n_epochs=200, load_all_data_GPU=False, margin=1., lr=None, optimizer=None, - scheduler=None, start_epoch_training=0, optimizer_kwargs={}, scheduler_kwargs={}, loader_kwargs={}): + scheduler=None, start_epoch_training=0, use_tqdm=True, optimizer_kwargs={}, scheduler_kwargs={}, + loader_kwargs={}): """ Parameters ---------- @@ -726,6 +731,8 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample the scheduler and the optimizer at each epoch. Default to 0. verbose: boolean, optional if True, prints more information from the training routine. Default to False. + use_tqdm : boolean, optional + Whether using tqdm or not to display progress. Defaults to True. optimizer_kwargs: Python dictionary, optional dictionary containing optional keyword arguments for the optimizer. scheduler_kwargs: Python dictionary, optional @@ -748,6 +755,7 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample quantile=quantile, batch_size=batch_size, n_epochs=n_epochs, load_all_data_GPU=load_all_data_GPU, margin=margin, lr=lr, optimizer=optimizer, scheduler=scheduler, + use_tqdm=use_tqdm, start_epoch_training=start_epoch_training, optimizer_kwargs=optimizer_kwargs, scheduler_kwargs=scheduler_kwargs, loader_kwargs=loader_kwargs) @@ -772,7 +780,7 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample early_stopping=False, epochs_early_stopping_interval=1, start_epoch_early_stopping=10, seed=None, cuda=None, scale_samples=True, quantile=0.1, batch_size=16, n_epochs=200, positive_weight=None, load_all_data_GPU=False, margin=1., lr=None, optimizer=None, scheduler=None, - start_epoch_training=0, optimizer_kwargs={}, scheduler_kwargs={}, loader_kwargs={}): + start_epoch_training=0, use_tqdm=True, optimizer_kwargs={}, scheduler_kwargs={}, loader_kwargs={}): """ Parameters ---------- @@ -851,6 +859,8 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample the scheduler and the optimizer at each epoch. Default to 0. verbose: boolean, optional if True, prints more information from the training routine. Default to False. + use_tqdm : boolean, optional + Whether using tqdm or not to display progress. Defaults to True. optimizer_kwargs: Python dictionary, optional dictionary containing optional keyword arguments for the optimizer. scheduler_kwargs: Python dictionary, optional @@ -876,6 +886,7 @@ def __init__(self, model, statistics_calc, backend, embedding_net=None, n_sample load_all_data_GPU=load_all_data_GPU, margin=margin, lr=lr, optimizer=optimizer, scheduler=scheduler, start_epoch_training=start_epoch_training, + use_tqdm=use_tqdm, optimizer_kwargs=optimizer_kwargs, scheduler_kwargs=scheduler_kwargs, loader_kwargs=loader_kwargs) diff --git a/abcpy/transformers.py b/abcpy/transformers.py new file mode 100644 index 00000000..95ba4675 --- /dev/null +++ b/abcpy/transformers.py @@ -0,0 +1,197 @@ +import numpy as np + + +# these transformers are used in the MCMC inference scheme, in order to run MCMC of an unbounded transformed space in +# case the original space is bounded. It therefore also implements the jacobian terms which appear in the acceptance +# rate. + +class BoundedVarTransformer: + """ + This scaler implements both lower bounded and two sided bounded transformations according to the provided bounds; + """ + + def __init__(self, lower_bound, upper_bound): + + # upper and lower bounds can be both scalar or array-like with size the size of the variable + self.lower_bound = lower_bound + self.upper_bound = upper_bound + + if not hasattr(lower_bound, "shape") or not hasattr(upper_bound, "shape"): + raise RuntimeError("Provided lower and upper bounds need to be arrays.") + elif hasattr(lower_bound, "shape") and hasattr(upper_bound, "shape") and lower_bound.shape != upper_bound.shape: + raise RuntimeError("Provided lower and upper bounds need to have same shape.") + + # note that == None checks if the array is None element wise. + self.unbounded_vars = np.logical_and(np.equal(lower_bound, None), np.equal(upper_bound, None)) + self.lower_bounded_vars = np.logical_and(np.not_equal(lower_bound, None), np.equal(upper_bound, None)) + self.upper_bounded_vars = np.logical_and(np.equal(lower_bound, None), np.not_equal(upper_bound, None)) + self.two_sided_bounded_vars = np.logical_and(np.not_equal(lower_bound, None), np.not_equal(upper_bound, None)) + if self.upper_bounded_vars.any(): + raise NotImplementedError("We do not yet implement the transformation for upper bounded random variables") + + self.lower_bound_lower_bounded = self.lower_bound[self.lower_bounded_vars].astype("float32") + self.lower_bound_two_sided = self.lower_bound[self.two_sided_bounded_vars].astype("float32") + self.upper_bound_two_sided = self.upper_bound[self.two_sided_bounded_vars].astype("float32") + + @staticmethod + def logit(x): + return np.log(x) - np.log(1 - x) + + def _check_data_in_bounds(self, x): + if np.any(x[self.lower_bounded_vars] <= self.lower_bound_lower_bounded): + raise RuntimeError("The provided data are out of the bounds.") + if (x[self.two_sided_bounded_vars] <= self.lower_bound[self.two_sided_bounded_vars]).any() or ( + x[self.two_sided_bounded_vars] >= self.upper_bound_two_sided).any(): + raise RuntimeError("The provided data is out of the bounds.") + + def _apply_nonlinear_transf(self, x): + # apply the different scalers to the different kind of variables: + x_transf = x.copy() + x_transf[self.lower_bounded_vars] = np.log(x[self.lower_bounded_vars] - self.lower_bound_lower_bounded) + x_transf[self.two_sided_bounded_vars] = self.logit( + (x[self.two_sided_bounded_vars] - self.lower_bound_two_sided) / ( + self.upper_bound_two_sided - self.lower_bound_two_sided)) + return x_transf + + @staticmethod + def _array_from_list(x): + return np.array(x).reshape(-1) + + @staticmethod + def _list_from_array(x_arr, x): + # transforms the array x to the list structure that contains x + x_new = [None] * len(x) + for i in range(len(x)): + if isinstance(x[i], np.ndarray): + x_new[i] = np.array(x_arr[i].reshape(x[i].shape)) + else: + x_new[i] = x_arr[i] + return x_new + + def transform(self, x): + """Scale features of x according to feature_range. + Parameters + ---------- + x : list of len n_parameters + Input data that will be transformed. + Returns + ------- + Xt : array-like of shape (n_samples, n_features) + Transformed data. + """ + # convert data to array: + x_arr = self._array_from_list(x) + + # need to check if we can apply the log first: + self._check_data_in_bounds(x_arr) + + # we transform the data with the log transformation: + x_arr = self._apply_nonlinear_transf(x_arr) + + # convert back to the list structure: + x = self._list_from_array(x_arr, x) + + return x + + def inverse_transform(self, x): + """Undo the scaling of x according to feature_range. + Parameters + ---------- + x : list of len n_parameters + Input data that will be transformed. It cannot be sparse. + Returns + ------- + Xt : array-like of shape (n_samples, n_features) + Transformed data. + """ + # now apply the inverse transform + x_arr = self._array_from_list(x) + + inv_x = x_arr.copy() + inv_x[self.two_sided_bounded_vars] = (self.upper_bound_two_sided - self.lower_bound_two_sided) * np.exp( + x_arr[self.two_sided_bounded_vars]) / (1 + np.exp( + x_arr[self.two_sided_bounded_vars])) + self.lower_bound_two_sided + inv_x[self.lower_bounded_vars] = np.exp(x_arr[self.lower_bounded_vars]) + self.lower_bound_lower_bounded + + # convert back to the list structure: + inv_x = self._list_from_array(inv_x, x) + + return inv_x + + def jac_log_det(self, x): + """Returns the log determinant of the Jacobian: log |J_t(x)|. + + Parameters + ---------- + x : list of len n_parameters + Input data, living in the original space (with lower bound constraints). + Returns + ------- + res : float + log determinant of the jacobian + """ + x = self._array_from_list(x) + self._check_data_in_bounds(x) + + results = np.zeros_like(x) + results[self.two_sided_bounded_vars] = np.log( + (self.upper_bound_two_sided - self.lower_bound_two_sided).astype("float64") / ( + (x[self.two_sided_bounded_vars] - self.lower_bound_two_sided) * ( + self.upper_bound_two_sided - x[self.two_sided_bounded_vars]))) + results[self.lower_bounded_vars] = - np.log(x[self.lower_bounded_vars] - self.lower_bound_lower_bounded) + + return np.sum(results) + + def jac_log_det_inverse_transform(self, x): + """Returns the log determinant of the Jacobian evaluated in the inverse transform: + log |J_t(t^{-1}(x))| = - log |J_{t^{-1}}(x)| + + Parameters + ---------- + x : list of len n_parameters + Input data, living in the transformed space (spanning the whole R^d). + Returns + ------- + res : float + log determinant of the jacobian evaluated in t^{-1}(x) + """ + x = self._array_from_list(x) + + results = np.zeros_like(x) + results[self.lower_bounded_vars] = - x[self.lower_bounded_vars] + # two sided: need some tricks to avoid numerical issues: + results[self.two_sided_bounded_vars] = - np.log( + self.upper_bound_two_sided - self.lower_bound_two_sided) + + indices = x[self.two_sided_bounded_vars] < 100 # for avoiding numerical overflow + res_b = np.copy(x)[self.two_sided_bounded_vars] + res_b[indices] = np.log(1 + np.exp(x[self.two_sided_bounded_vars][indices])) + results[self.two_sided_bounded_vars] += res_b + + indices = x[self.two_sided_bounded_vars] > - 100 # for avoiding numerical overflow + res_c = np.copy(- x)[self.two_sided_bounded_vars] + res_c[indices] = np.log(1 + np.exp(- x[self.two_sided_bounded_vars][indices])) + results[self.two_sided_bounded_vars] += res_c + + # res = res_b + res_c - res_a + + return np.sum(results) + + +class DummyTransformer: + """Dummy transformer which does nothing, and for which the jacobian is 1""" + + def __init__(self): + pass + + def transform(self, x): + return x + + def inverse_transform(self, x): + return x + + def jac_log_det(self, x): + return 0 + + def jac_log_det_inverse_transform(self, x): + return 0 diff --git a/doc/source/abcpy.rst b/doc/source/abcpy.rst index 4f332e70..eb4c6974 100644 --- a/doc/source/abcpy.rst +++ b/doc/source/abcpy.rst @@ -196,3 +196,12 @@ abcpy.statisticslearning module :special-members: __init__ :undoc-members: :show-inheritance: + +abcpy.transformers module +------------------------------- + +.. automodule:: abcpy.transformers + :members: + :special-members: __init__ + :undoc-members: + :show-inheritance: diff --git a/doc/source/class-diagram.png b/doc/source/class-diagram.png index 767fbff9..33517d0c 100644 Binary files a/doc/source/class-diagram.png and b/doc/source/class-diagram.png differ diff --git a/doc/source/class-diagram.svg b/doc/source/class-diagram.svg new file mode 100644 index 00000000..997e28ef --- /dev/null +++ b/doc/source/class-diagram.svg @@ -0,0 +1,7130 @@ + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ContrastiveDistLearn + SemiautomaticNN + TripletDistLearn + + + NeuralEmbedding + LinearTransformation + + + + + + + + + + + + + + + + StatisticsLearning + Hyperparameter + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + MCMCMetropolisHastings + + + + + SemiParametricBSL + + + + + + + + + LogNormal + + + + Exponential + + + + + + + + + + + + SlicedWasserstein + GammaDivergence + KLDivergence + MMD + EnergyDistance + SquaredHellinger + + + + + + + + + + + DiscreteUniform + + + + + diff --git a/doc/source/getting_started.rst b/doc/source/getting_started.rst index 08d14e30..392e5a1a 100644 --- a/doc/source/getting_started.rst +++ b/doc/source/getting_started.rst @@ -289,24 +289,31 @@ In ABCpy, we implement widely used and advanced variants of ABC inferential sche * Rejection ABC :py:class:`abcpy.inferences.RejectionABC`, * Population Monte Carlo ABC :py:class:`abcpy.inferences.PMCABC`, -* Sequential Monte Carlo ABC :py:class:`abcpy.inferences.SMCABC`, +* Two versions of Sequential Monte Carlo ABC :py:class:`abcpy.inferences.SMCABC`, * Replenishment sequential Monte Carlo ABC (RSMC-ABC) :py:class:`abcpy.inferences.RSMCABC`, * Adaptive population Monte Carlo ABC (APMC-ABC) :py:class:`abcpy.inferences.APMCABC`, * ABC with subset simulation (ABCsubsim) :py:class:`abcpy.inferences.ABCsubsim`, and * Simulated annealing ABC (SABC) :py:class:`abcpy.inferences.SABC`. -To perform ABC algorithms, we provide different standard distance functions between datasets, e.g., a discrepancy -measured by achievable classification accuracy between two datasets +To perform ABC algorithms, we provide different the standard Euclidean distance function +(:py:class:`abcpy.distances.Euclidean`) as well as discrepancy measures between empirical datasets +(eg, achievable classification accuracy between two datasets -* :py:class:`abcpy.distances.Euclidean`, -* :py:class:`abcpy.distances.LogReg`, -* :py:class:`abcpy.distances.PenLogReg`. +* classification accuracy (:py:class:`abcpy.distances.LogReg`, :py:class:`abcpy.distances.PenLogReg`) +* :py:class:`abcpy.distances.Wasserstein` +* :py:class:`abcpy.distances.SlicedWasserstein` +* :py:class:`abcpy.distances.GammaDivergence` +* :py:class:`abcpy.distances.KLDivergence` +* :py:class:`abcpy.distances.MMD` +* :py:class:`abcpy.distances.EnergyDistance` +* :py:class:`abcpy.distances.SquaredHellingerDistance`. -We also have implemented the population Monte Carlo :py:class:`abcpy.inferences.PMC` algorithm to infer parameters when +We also have implemented the standard Metropolis-Hastings MCMC :py:class:`abcpy.inferences.MCMCMetropoliHastings` and population Monte Carlo :py:class:`abcpy.inferences.PMC` algorithm to infer parameters when the likelihood or approximate likelihood function is available. For approximation of the likelihood function we provide -two methods: +the following methods methods: -* Synthetic likelihood approximation :py:class:`abcpy.approx_lhd.SynLikelihood`, and another method using +* Synthetic likelihood approximation :py:class:`abcpy.approx_lhd.SynLikelihood`, +* Semiparametric Synthetic likelihood :py:class:`abcpy.approx_lhd.SemiParametricSynLikelihood`, and another method using * penalized logistic regression :py:class:`abcpy.approx_lhd.PenLogReg`. Next we explain how we can use PMC algorithm using approximation of the @@ -336,6 +343,7 @@ Further possibilities of combination will be made available in later versions of The source code can be found in `examples/approx_lhd/pmc_hierarchical_models.py`. +For an example using MCMC instead of PMC, check out `examples/approx_lhd/mcmc_hierarchical_models.py`. Statistics Learning ~~~~~~~~~~~~~~~~~~~ diff --git a/doc/source/installation.rst b/doc/source/installation.rst index 9b2c0708..5e8c4e28 100644 --- a/doc/source/installation.rst +++ b/doc/source/installation.rst @@ -36,7 +36,7 @@ To create a package and install it, do make package pip3 install wheel - pip3 install build/dist/abcpy-0.6.1-py3-none-any.whl + pip3 install build/dist/abcpy-0.6.2-py3-none-any.whl ``wheel`` is required to install in this way. diff --git a/examples/approx_lhd/mcmc_hierarchical_models.py b/examples/approx_lhd/mcmc_hierarchical_models.py new file mode 100644 index 00000000..72c17841 --- /dev/null +++ b/examples/approx_lhd/mcmc_hierarchical_models.py @@ -0,0 +1,132 @@ +"""An example showing how to implement a bayesian network in ABCpy. We consider here a model of school grades which +depend on some variables.""" +import logging + + +def infer_parameters(n_samples=500, burnin=250, n_samples_per_param=40, logging_level=logging.WARN): + """Perform inference for this example. + + Parameters + ---------- + n_samples : integer, optional + Number of posterior samples to generate. The default value is 500. + burnin : integer, optional + Number of burnin steps for the MCMC. The default value is 250. + n_samples_per_param : integer, optional + Number of data points in each simulated data set. The default value is 10. + + Returns + ------- + abcpy.output.Journal + A journal containing simulation results, metadata and optionally intermediate results. + """ + logging.basicConfig(level=logging_level) + # Observed data corresponding to model_1 defined below + grades_obs = [3.872486707973337, 4.6735380808674405, 3.9703538990858376, 4.11021272048805, 4.211048655421368, + 4.154817956586653, 4.0046893064392695, 4.01891381384729, 4.123804757702919, 4.014941267301294, + 3.888174595940634, 4.185275142948246, 4.55148774469135, 3.8954427675259016, 4.229264035335705, + 3.839949451328312, 4.039402553532825, 4.128077814241238, 4.361488645531874, 4.086279074446419, + 4.370801602256129, 3.7431697332475466, 4.459454162392378, 3.8873973643008255, 4.302566721487124, + 4.05556051626865, 4.128817316703757, 3.8673704442215984, 4.2174459453805015, 4.202280254493361, + 4.072851400451234, 3.795173229398952, 4.310702877332585, 4.376886328810306, 4.183704734748868, + 4.332192463368128, 3.9071312388426587, 4.311681374107893, 3.55187913252144, 3.318878360783221, + 4.187850500877817, 4.207923106081567, 4.190462065625179, 4.2341474252986036, 4.110228694304768, + 4.1589891480847765, 4.0345604687633045, 4.090635481715123, 3.1384654393449294, 4.20375641386518, + 4.150452690356067, 4.015304457401275, 3.9635442007388195, 4.075915739179875, 3.5702080541929284, + 4.722333310410388, 3.9087618197155227, 4.3990088006390735, 3.968501165774181, 4.047603645360087, + 4.109184340976979, 4.132424805281853, 4.444358334346812, 4.097211737683927, 4.288553086265748, + 3.8668863066511303, 3.8837108501541007] + + # The prior information changing the class size and social background, depending on school location + from abcpy.continuousmodels import Uniform, Normal + school_location = Uniform([[0.2], [0.3]], name="school_location") + + # The average class size of a certain school + class_size = Normal([[school_location], [0.1]], name="class_size") + + # The social background of a student + background = Normal([[school_location], [0.1]], name="background") + + # The grade a student would receive without any bias + grade_without_additional_effects = Normal([[4.5], [0.25]], name="grade_without_additional_effects") + + # The grade a student of a certain school receives; this defined a new random variable by subtraction + final_grade = grade_without_additional_effects - class_size - background + + # Observed data corresponding to model_2 defined below + scholarship_obs = [2.7179657436207805, 2.124647285937229, 3.07193407853297, 2.335024761813643, 2.871893855192, + 3.4332002458233837, 3.649996835818173, 3.50292335102711, 2.815638168018455, 2.3581613289315992, + 2.2794821846395568, 2.8725835459926503, 3.5588573782815685, 2.26053126526137, 1.8998143530749971, + 2.101110815311782, 2.3482974964831573, 2.2707679029919206, 2.4624550491079225, 2.867017757972507, + 3.204249152084959, 2.4489542437714213, 1.875415915801106, 2.5604889644872433, 3.891985093269989, + 2.7233633223405205, 2.2861070389383533, 2.9758813233490082, 3.1183403287267755, + 2.911814060853062, 2.60896794303205, 3.5717098647480316, 3.3355752461779824, 1.99172284546858, + 2.339937680892163, 2.9835630207301636, 2.1684912355975774, 3.014847335983034, 2.7844122961916202, + 2.752119871525148, 2.1567428931391635, 2.5803629307680644, 2.7326646074552103, 2.559237193255186, + 3.13478196958166, 2.388760269933492, 3.2822443541491815, 2.0114405441787437, 3.0380056368041073, + 2.4889680313769724, 2.821660164621084, 3.343985964873723, 3.1866861970287808, 4.4535037154856045, + 3.0026333138006027, 2.0675706089352612, 2.3835301730913185, 2.584208398359566, 3.288077633446465, + 2.6955853384148183, 2.918315169739928, 3.2464814419322985, 2.1601516779909433, 3.231003347780546, + 1.0893224045062178, 0.8032302688764734, 2.868438615047827] + + # A quantity that determines whether a student will receive a scholarship + scholarship_without_additional_effects = Normal([[2], [0.5]], name="scholarship_without_additional_effects") + + # A quantity determining whether a student receives a scholarship, including his social background + final_scholarship = scholarship_without_additional_effects + 3 * background + + # Define a summary statistics for final grade and final scholarship + from abcpy.statistics import Identity + statistics_calculator_final_grade = Identity(degree=2, cross=False) + statistics_calculator_final_scholarship = Identity(degree=3, cross=False) + + # Define an approximate likelihood for final grade and final scholarship + from abcpy.approx_lhd import SynLikelihood + approx_lhd_final_grade = SynLikelihood(statistics_calculator_final_grade) + approx_lhd_final_scholarship = SynLikelihood(statistics_calculator_final_scholarship) + + # Define a backend + from abcpy.backends import BackendDummy as Backend + backend = Backend() + + # Define a perturbation kernel to explore parameter space + from abcpy.perturbationkernel import DefaultKernel + kernel = DefaultKernel([school_location, class_size, grade_without_additional_effects, + background, scholarship_without_additional_effects]) + + # Define sampler to use with the + from abcpy.inferences import MCMCMetropoliHastings + sampler = MCMCMetropoliHastings([final_grade, final_scholarship], + [approx_lhd_final_grade, approx_lhd_final_scholarship], backend, kernel, seed=1) + + # Sample + journal = sampler.sample([grades_obs, scholarship_obs], n_samples=n_samples, burnin=burnin, + n_samples_per_param=n_samples_per_param) + return journal + + +def analyse_journal(journal): + # output parameters and weights + print(journal.get_parameters()) + print(journal.weights) + + # do post analysis + print(journal.posterior_mean()) + print(journal.posterior_cov()) + + # print configuration + print(journal.configuration) + + # plot posterior + journal.plot_posterior_distr(path_to_save="posterior.png") + + # save and load journal + journal.save("experiments.jnl") + + from abcpy.output import Journal + new_journal = Journal.fromFile('experiments.jnl') + + +if __name__ == "__main__": + journal = infer_parameters(logging_level=logging.INFO) + analyse_journal(journal) diff --git a/setup.py b/setup.py index 21b36bcb..608a6104 100644 --- a/setup.py +++ b/setup.py @@ -63,6 +63,7 @@ 'Programming Language :: Python :: 3.6', 'Programming Language :: Python :: 3.7', 'Programming Language :: Python :: 3.8', + 'Programming Language :: Python :: 3.9', ], # What does your project relate to? diff --git a/tests/approx_lhd_tests.py b/tests/approx_lhd_tests.py index 2942f007..4b14d14d 100644 --- a/tests/approx_lhd_tests.py +++ b/tests/approx_lhd_tests.py @@ -2,7 +2,7 @@ import numpy as np -from abcpy.approx_lhd import PenLogReg, SynLikelihood +from abcpy.approx_lhd import PenLogReg, SynLikelihood, SemiParametricSynLikelihood from abcpy.continuousmodels import Normal from abcpy.continuousmodels import Uniform from abcpy.statistics import Identity @@ -56,13 +56,20 @@ def test_likelihood(self): self.assertAlmostEqual(comp_likelihood_biv, np.log(expected_likelihood_biv)) def test_likelihood_multiple_observations(self): - comp_likelihood = self.likfun.loglikelihood(self.y_obs, self.y_sim) - expected_likelihood = 9.77317308598673e-08 - self.assertAlmostEqual(comp_likelihood, np.log(expected_likelihood)) + comp_likelihood = self.likfun.likelihood(self.y_obs_double, self.y_sim) + expected_likelihood = 7.337876253225462e-10 + self.assertAlmostEqual(comp_likelihood, expected_likelihood) expected_likelihood_biv = 0.9999999999999979 - comp_likelihood_biv = self.likfun_bivariate.loglikelihood(self.y_obs_bivariate, self.y_sim_bivariate) - self.assertAlmostEqual(comp_likelihood_biv, np.log(expected_likelihood_biv)) + comp_likelihood_biv = self.likfun_bivariate.likelihood(self.y_obs_bivariate_double, self.y_sim_bivariate) + self.assertAlmostEqual(comp_likelihood_biv, expected_likelihood_biv) + + def test_loglikelihood_additive(self): + comp_loglikelihood_a = self.likfun.loglikelihood([self.y_obs_double[0]], self.y_sim) + comp_loglikelihood_b = self.likfun.loglikelihood([self.y_obs_double[1]], self.y_sim) + comp_loglikelihood_two = self.likfun.loglikelihood(self.y_obs_double, self.y_sim) + + self.assertAlmostEqual(comp_loglikelihood_two, comp_loglikelihood_a + comp_loglikelihood_b) class SynLikelihoodTests(unittest.TestCase): @@ -85,18 +92,72 @@ def test_likelihood(self): # create observed data y_obs = [1.8] # calculate the statistics of the observed data - comp_likelihood = self.likfun.loglikelihood(y_obs, self.y_sim) - expected_likelihood = 0.20963610211945238 + comp_loglikelihood = self.likfun.loglikelihood(y_obs, self.y_sim) + expected_loglikelihood = -0.6434435652263701 # This checks whether it computes a correct value and dimension is right - self.assertAlmostEqual(comp_likelihood, np.log(expected_likelihood)) + self.assertAlmostEqual(comp_loglikelihood, expected_loglikelihood) def test_likelihood_multiple_observations(self): y_obs = [1.8, 0.9] - comp_likelihood = self.likfun.loglikelihood(y_obs, self.y_sim) - print(comp_likelihood) - expected_likelihood = 0.04457899184856649 + comp_loglikelihood = self.likfun.loglikelihood(y_obs, self.y_sim) + expected_loglikelihood = -1.2726154993040115 # This checks whether it computes a correct value and dimension is right - self.assertAlmostEqual(comp_likelihood, np.log(expected_likelihood)) + self.assertAlmostEqual(comp_loglikelihood, expected_loglikelihood) + + def test_loglikelihood_additive(self): + y_obs = [1.8, 0.9] + comp_loglikelihood_a = self.likfun.loglikelihood([y_obs[0]], self.y_sim) + comp_loglikelihood_b = self.likfun.loglikelihood([y_obs[1]], self.y_sim) + comp_loglikelihood_two = self.likfun.loglikelihood(y_obs, self.y_sim) + + self.assertAlmostEqual(comp_loglikelihood_two, comp_loglikelihood_a + comp_loglikelihood_b) + + +class SemiParametricSynLikelihoodTests(unittest.TestCase): + def setUp(self): + self.mu = Uniform([[-5.0], [5.0]], name='mu') + self.sigma = Uniform([[5.0], [10.0]], name='sigma') + self.model = Normal([self.mu, self.sigma]) + self.stat_calc_1 = Identity(degree=1, cross=False) + self.likfun_1 = SemiParametricSynLikelihood(self.stat_calc_1) + self.stat_calc = Identity(degree=2, cross=False) + self.likfun = SemiParametricSynLikelihood(self.stat_calc) + # create fake simulated data + self.mu._fixed_values = [1.1] + self.sigma._fixed_values = [1.0] + self.y_sim = self.model.forward_simulate(self.model.get_input_values(), 100, rng=np.random.RandomState(1)) + + def test_likelihood(self): + # Checks whether wrong input type produces error message + self.assertRaises(TypeError, self.likfun.loglikelihood, 3.4, [2, 1]) + self.assertRaises(TypeError, self.likfun.loglikelihood, [2, 4], 3.4) + + # create observed data + y_obs = [1.8] + + # check whether it raises correct error with input of wrong size + self.assertRaises(RuntimeError, self.likfun_1.loglikelihood, y_obs, self.y_sim) + + # calculate the statistics of the observed data + comp_loglikelihood = self.likfun.loglikelihood(y_obs, self.y_sim) + expected_loglikelihood = -2.3069321875272815 + # This checks whether it computes a correct value and dimension is right + self.assertAlmostEqual(comp_loglikelihood, expected_loglikelihood) + + def test_likelihood_multiple_observations(self): + y_obs = [1.8, 0.9] + comp_loglikelihood = self.likfun.loglikelihood(y_obs, self.y_sim) + expected_loglikelihood = -3.7537571275591683 + # This checks whether it computes a correct value and dimension is right + self.assertAlmostEqual(comp_loglikelihood, expected_loglikelihood) + + def test_loglikelihood_additive(self): + y_obs = [1.8, 0.9] + comp_loglikelihood_a = self.likfun.loglikelihood([y_obs[0]], self.y_sim) + comp_loglikelihood_b = self.likfun.loglikelihood([y_obs[1]], self.y_sim) + comp_loglikelihood_two = self.likfun.loglikelihood(y_obs, self.y_sim) + + self.assertAlmostEqual(comp_loglikelihood_two, comp_loglikelihood_a + comp_loglikelihood_b) if __name__ == '__main__': diff --git a/tests/continuousmodels_tests.py b/tests/continuousmodels_tests.py index ac35f9a7..c3e910f8 100644 --- a/tests/continuousmodels_tests.py +++ b/tests/continuousmodels_tests.py @@ -31,6 +31,16 @@ class MultiStudentTAPITests(AbstractAPIImplementationTests, unittest.TestCase): model_inputs = [[[1, 0], [[1, 0], [0, 1]], 3]] +class LogNormalTAPITests(AbstractAPIImplementationTests, unittest.TestCase): + model_types = [LogNormal] + model_inputs = [[0, 1]] + + +class ExponentialTAPITests(AbstractAPIImplementationTests, unittest.TestCase): + model_types = [Exponential] + model_inputs = [[0.4]] + + class CheckParametersAtInitializationTests(unittest.TestCase): """Tests that no probabilistic model with invalid parameters can be initialized.""" @@ -72,6 +82,14 @@ def test_MultiStudentT(self): with self.assertRaises(ValueError): MultiStudentT([[1, 0], [[1, 0], [0, 1]], -1]) + def test_LogNormal(self): + with self.assertRaises(ValueError): + LogNormal([1, -1]) + + def test_Exponential(self): + with self.assertRaises(ValueError): + Exponential([[1], [-1]]) + class DimensionTests(unittest.TestCase): """Tests whether the dimensions of all continuous models are defined in the correct way.""" @@ -96,6 +114,14 @@ def test_MultiStudentT(self): M = MultiStudentT([[1, 0], [[0.1, 0], [0, 0.1]], 1]) self.assertTrue(M.get_output_dimension() == 2) + def test_LogNormal(self): + LN = LogNormal([3, 1]) + self.assertTrue(LN.get_output_dimension() == 1) + + def test_LogNormal(self): + EXP = Exponential([3]) + self.assertTrue(EXP.get_output_dimension() == 1) + class SampleFromDistributionTests(unittest.TestCase): """Tests the return value of forward_simulate for all continuous distributions.""" @@ -130,6 +156,18 @@ def test_Uniform(self): self.assertTrue(isinstance(samples, list)) self.assertTrue(len(samples) == 3) + def test_LogNormal(self): + LN = LogNormal([3, 1]) + samples = LN.forward_simulate(LN.get_input_values(), 3) + self.assertTrue(isinstance(samples, list)) + self.assertTrue(len(samples) == 3) + + def test_LogNormal(self): + EXP = Exponential([3]) + samples = EXP.forward_simulate(EXP.get_input_values(), 3) + self.assertTrue(isinstance(samples, list)) + self.assertTrue(len(samples) == 3) + class CheckParametersBeforeSamplingTests(unittest.TestCase): """Tests whether False will be returned if the input parameters of _check_parameters_before_sampling are not accepted.""" @@ -162,6 +200,15 @@ def test_MultiStudentT(self): self.assertFalse(M._check_input([[1, 0], [[1, 0], [0, 1]], -1])) + def test_LogNormal(self): + LN = LogNormal([3, 1]) + self.assertFalse(LN._check_input([3, -1])) + + def test_Exponential(self): + EXP = Exponential([3]) + self.assertFalse(EXP._check_input([-3])) + self.assertFalse(EXP._check_input([-3, 1])) + if __name__ == '__main__': unittest.main() diff --git a/tests/distances_tests.py b/tests/distances_tests.py index 349d25a6..c75cf7ce 100644 --- a/tests/distances_tests.py +++ b/tests/distances_tests.py @@ -1,8 +1,9 @@ -import unittest - import numpy as np +import unittest +import warnings -from abcpy.distances import Euclidean, PenLogReg, LogReg, Wasserstein +from abcpy.distances import Euclidean, PenLogReg, LogReg, Wasserstein, SlicedWasserstein, GammaDivergence, MMD, \ + EnergyDistance, KLDivergence, SquaredHellingerDistance from abcpy.statistics import Identity @@ -111,7 +112,6 @@ def test_distance(self): self.assertRaises(TypeError, self.distancefunc.distance, 3.4, d2) self.assertRaises(TypeError, self.distancefunc.distance, d1, 3.4) - # completely separable datasets should have a distance of 1.0 self.assertEqual(self.distancefunc.distance(d1, d2), 28.623685155319652) # equal data sets should have a distance of approximately 0.0; it won't be exactly 0 due to numerical rounding @@ -121,5 +121,255 @@ def test_dist_max(self): self.assertTrue(self.distancefunc.dist_max() == np.inf) +class SlicedWassersteinTests(unittest.TestCase): + def setUp(self): + self.stat_calc = Identity(degree=2, cross=False) + self.rng = np.random.RandomState(1) + self.distancefunc = SlicedWasserstein(self.stat_calc, rng=self.rng) + + def test_distance(self): + d1 = 0.5 * self.rng.randn(100, 2) - 10 + d2 = 0.5 * self.rng.randn(100, 2) + 10 + + d1 = d1.tolist() + d2 = d2.tolist() + + # Checks whether wrong input type produces error message + self.assertRaises(TypeError, self.distancefunc.distance, 3.4, d2) + self.assertRaises(TypeError, self.distancefunc.distance, d1, 3.4) + + self.assertAlmostEqual(self.distancefunc.distance(d1, d2), 12.604402810464576) + + # equal data sets should have a distance of approximately 0.0; it won't be exactly 0 due to numerical rounding + self.assertAlmostEqual(self.distancefunc.distance(d1, d1), 0.0, delta=1e-5) + + def test_dist_max(self): + self.assertTrue(self.distancefunc.dist_max() == np.inf) + + +class GammaDivergenceTests(unittest.TestCase): + def setUp(self): + self.stat_calc = Identity(degree=2, cross=False) + self.rng = np.random.RandomState(1) + # notice: the GammaDivergence estimator becomes infinity when one element in the s1 dataset is equal to one in + # the s2 dataset and k=1 (as the distance between those two would be 0, which gives infinity when dividing) + self.distancefunc_k1 = GammaDivergence(self.stat_calc, gam=0.1, k=1) + self.distancefunc_k2 = GammaDivergence(self.stat_calc, gam=0.1, k=2) + + def test_distance(self): + d1 = 0.5 * self.rng.randn(100, 2) - 10 + d2 = 0.5 * self.rng.randn(100, 2) + 10 + d3 = 0.5 * self.rng.randn(100, 2) + 10 + d3[0] = d1[0] # set one single element to be equal + d4 = 0.5 * self.rng.randn(100, 2) + 10 + d4[0] = d4[1] = d1[0] # set two elements to be equal + d5 = 0.5 * self.rng.randn(100, 2) + 10 + d5[0] = d5[1] = d5[2] = d1[0] # set three elements to be equal + + d1 = d1.tolist() + d2 = d2.tolist() + d3 = d3.tolist() + d4 = d4.tolist() + d5 = d5.tolist() + + # Checks whether wrong input type produces error message + self.assertRaises(TypeError, self.distancefunc_k1.distance, 3.4, d2) + self.assertRaises(TypeError, self.distancefunc_k1.distance, d1, 3.4) + + self.assertRaises(ValueError, self.distancefunc_k2.distance, [d1[i] for i in range(2)], + [d2[i] for i in range(2)]) + + # print(self.distancefunc.distance(d1, d2)) + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d2), 11.781045117610248) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d2), 9.898652219975974) + + # check identical dataset + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d1), - np.inf) + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d3), - np.inf) + self.assertAlmostEqual(self.distancefunc_k1.distance(d2, d4), + np.inf) + self.assertAlmostEqual(self.distancefunc_k2.distance(d2, d5), + np.inf) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d1), -1.8398859674447485) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d3), 9.904082614090905) + + def test_dist_max(self): + self.assertTrue(self.distancefunc_k1.dist_max() == np.inf) + + +class KLDivergenceTests(unittest.TestCase): + def setUp(self): + self.stat_calc = Identity(degree=2, cross=False) + self.rng = np.random.RandomState(1) + # notice: the KLDivergenceTests estimator becomes infinity when one element in the s1 dataset is equal to one + # in the s2 dataset and k=1 (as the distance between those two would be 0, which gives infinity when dividing) + self.distancefunc_k1 = KLDivergence(self.stat_calc, k=1) + self.distancefunc_k2 = KLDivergence(self.stat_calc, k=2) + + def test_distance(self): + d1 = 0.5 * self.rng.randn(100, 2) - 10 + d2 = 0.5 * self.rng.randn(100, 2) + 10 + d3 = 0.5 * self.rng.randn(100, 2) + 10 + d3[0] = d1[0] # set one single element to be equal + + d1 = d1.tolist() + d2 = d2.tolist() + d3 = d3.tolist() + + # Checks whether wrong input type produces error message + self.assertRaises(TypeError, self.distancefunc_k1.distance, 3.4, d2) + self.assertRaises(TypeError, self.distancefunc_k1.distance, d1, 3.4) + + self.assertRaises(ValueError, self.distancefunc_k2.distance, [d1[i] for i in range(2)], + [d2[i] for i in range(2)]) + + # print(self.distancefunc.distance(d1, d2)) + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d2), 11.268656821978732) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d2), 9.632589300292643) + + # check identical dataset + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d1), - np.inf) + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d3), - np.inf) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d1), -1.6431603119265736) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d3), 9.629190258363401) + + def test_dist_max(self): + self.assertTrue(self.distancefunc_k1.dist_max() == np.inf) + + +class MMDTests(unittest.TestCase): + def setUp(self): + self.stat_calc = Identity(degree=2, cross=False) + self.rng = np.random.RandomState(1) + self.distancefunc = MMD(self.stat_calc) + self.distancefunc_biased = MMD(self.stat_calc, biased_estimator=True) + + def test_initialization(self): + self.assertRaises(RuntimeError, MMD, self.stat_calc, 5) + self.assertRaises(NotImplementedError, MMD, self.stat_calc, "ciao") + + def test_distance(self): + d1 = 0.5 * self.rng.randn(100, 2) - 10 + d2 = 0.5 * self.rng.randn(100, 2) + 10 + + d1 = d1.tolist() + d2 = d2.tolist() + + # Checks whether wrong input type produces error message + self.assertRaises(TypeError, self.distancefunc.distance, 3.4, d2) + self.assertRaises(TypeError, self.distancefunc.distance, d1, 3.4) + + # print(self.distancefunc.distance(d1, d2)) + self.assertAlmostEqual(self.distancefunc.distance(d1, d2), 0.01078569482366823) + self.assertAlmostEqual(self.distancefunc_biased.distance(d1, d2), 0.030677837875431546) + + # check identical dataset + self.assertAlmostEqual(self.distancefunc.distance(d1, d1), - 0.019872124487359623) + self.assertAlmostEqual(self.distancefunc_biased.distance(d1, d1), 0.0) + + # check if it is symmetric: + self.assertAlmostEqual(self.distancefunc.distance(d1, d2), self.distancefunc.distance(d2, d1)) + self.assertAlmostEqual(self.distancefunc_biased.distance(d1, d2), self.distancefunc_biased.distance(d2, d1)) + + def test_dist_max(self): + self.assertTrue(self.distancefunc.dist_max() == np.inf) + + +class EnergyDistanceTests(unittest.TestCase): + def setUp(self): + self.stat_calc = Identity(degree=2, cross=False) + self.rng = np.random.RandomState(1) + self.distancefunc = EnergyDistance(self.stat_calc) + self.distancefunc_unbiased = EnergyDistance(self.stat_calc, biased_estimator=False) + + def test_initialization(self): + self.assertRaises(RuntimeError, MMD, self.stat_calc, 5) + self.assertRaises(NotImplementedError, MMD, self.stat_calc, "ciao") + + def test_distance(self): + d1 = 0.5 * self.rng.randn(100, 2) - 10 + d2 = 0.5 * self.rng.randn(100, 2) + 10 + + d1 = d1.tolist() + d2 = d2.tolist() + + # Checks whether wrong input type produces error message + self.assertRaises(TypeError, self.distancefunc.distance, 3.4, d2) + self.assertRaises(TypeError, self.distancefunc.distance, d1, 3.4) + + # print(self.distancefunc.distance(d1, d2)) + self.assertAlmostEqual(self.distancefunc.distance(d1, d2), 33.76256466781305) + self.assertAlmostEqual(self.distancefunc_unbiased.distance(d1, d2), 33.41913579594933) + + # check identical dataset + self.assertAlmostEqual(self.distancefunc.distance(d1, d1), 0.0) + self.assertAlmostEqual(self.distancefunc_unbiased.distance(d1, d1), -0.31969528899748667) + + # check if it is symmetric: + self.assertAlmostEqual(self.distancefunc.distance(d1, d2), self.distancefunc.distance(d2, d1)) + self.assertAlmostEqual(self.distancefunc_unbiased.distance(d1, d2), + self.distancefunc_unbiased.distance(d2, d1), ) + + def test_dist_max(self): + self.assertTrue(self.distancefunc.dist_max() == np.inf) + + +class SquaredHellingerDistanceTests(unittest.TestCase): + def setUp(self): + self.stat_calc = Identity(degree=1, cross=False) + self.rng = np.random.RandomState(1) + # notice: the SquaredHellingerDistance estimator becomes infinity when one element in the s1 dataset is equal to one in + # the s2 dataset and k=1 (as the distance between those two would be 0, which gives infinity when dividing) + self.distancefunc_k1 = SquaredHellingerDistance(self.stat_calc, k=1) + self.distancefunc_k2 = SquaredHellingerDistance(self.stat_calc, k=2) + + # def test_different_distances(self): + # d1 = 0.5 * self.rng.randn(100, 2) - 10 + # d2 = 0.5 * self.rng.randn(100, 2) - 10 + # d1 = d1.tolist() + # for i in range(-20, 40): + # d2_list = (d2 + i).tolist() + # print(i, self.distancefunc_k1.distance(d1, d2_list)) + # + # # the estimator seems to be OK, even if it reaches super small values when the two datasets are from same + # # parameters + + def test_distance(self): + d1 = 0.5 * self.rng.randn(100, 2) - 10 + d2 = 0.5 * self.rng.randn(100, 2) + 10 + d3 = 0.5 * self.rng.randn(100, 2) + 10 + d3[0] = d1[0] # set one single element to be equal + d1 = d1.tolist() + d2 = d2.tolist() + d3 = d3.tolist() + + # Checks whether wrong input type produces error message + self.assertRaises(TypeError, self.distancefunc_k1.distance, 3.4, d2) + self.assertRaises(TypeError, self.distancefunc_k1.distance, d1, 3.4) + + self.assertRaises(ValueError, self.distancefunc_k2.distance, [d1[i] for i in range(2)], + [d2[i] for i in range(2)]) + + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d2), 1.123969688201248) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d2), 0.7593257032591699) + + # check infinities when there are same elements in the two datasets + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d1), - np.inf) + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d3), - np.inf) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d1), -324.97681257084014) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d3), -274.4911947796949) + + # check if it is symmetric: + self.assertAlmostEqual(self.distancefunc_k1.distance(d1, d2), self.distancefunc_k1.distance(d2, d1)) + self.assertAlmostEqual(self.distancefunc_k2.distance(d1, d2), self.distancefunc_k2.distance(d2, d1)) + + def test_dist_max(self): + self.assertTrue(self.distancefunc_k1.dist_max() == 2) + + if __name__ == '__main__': unittest.main() diff --git a/tests/inferences_tests.py b/tests/inferences_tests.py index ffd1ce3e..da824ae3 100644 --- a/tests/inferences_tests.py +++ b/tests/inferences_tests.py @@ -6,13 +6,14 @@ from abcpy.backends import BackendDummy from abcpy.continuousmodels import Normal from abcpy.continuousmodels import Uniform -from abcpy.distances import Euclidean -from abcpy.inferences import RejectionABC, PMC, PMCABC, SABC, ABCsubsim, SMCABC, APMCABC, RSMCABC +from abcpy.distances import Euclidean, MMD +from abcpy.inferences import RejectionABC, PMC, PMCABC, SABC, ABCsubsim, SMCABC, APMCABC, RSMCABC, \ + MCMCMetropoliHastings from abcpy.statistics import Identity class RejectionABCTest(unittest.TestCase): - def test_sample(self): + def setUp(self): # setup backend dummy = BackendDummy() @@ -23,33 +24,85 @@ def test_sample(self): self.model = Normal([mu, sigma]) # define sufficient statistics for the model - stat_calc = Identity(degree=2, cross=0) + stat_calc = Identity(degree=2, cross=False) # define a distance function dist_calc = Euclidean(stat_calc) # create fake observed data - y_obs = [np.array(9.8)] + self.y_obs = [np.array(9.8)] + # for correct seeding define 2 samplers + self.sampler = RejectionABC([self.model], [dist_calc], dummy, seed=1) + self.sampler2 = RejectionABC([self.model], [dist_calc], dummy, seed=1) + + def test_sample_n_samples(self): # use the rejection sampling scheme - sampler = RejectionABC([self.model], [dist_calc], dummy, seed=1) - journal = sampler.sample([y_obs], 10, 1, 10) + journal = self.sampler.sample([self.y_obs], 10, 1, 10, path_to_save_journal="tmp.jnl") mu_sample = np.array(journal.get_parameters()['mu']) sigma_sample = np.array(journal.get_parameters()['sigma']) # test shape of samples mu_shape, sigma_shape = (len(mu_sample), mu_sample[0].shape[1]), \ - (len(sigma_sample), - sigma_sample[0].shape[1]) + (len(sigma_sample), sigma_sample[0].shape[1]) self.assertEqual(mu_shape, (10, 1)) self.assertEqual(sigma_shape, (10, 1)) # Compute posterior mean # self.assertAlmostEqual(np.average(np.asarray(samples[:,0])),1.22301,10e-2) - self.assertLess(np.average(mu_sample) - 1.22301, 1e-2) - self.assertLess(np.average(sigma_sample) - 6.992218, 10e-2) + self.assertAlmostEqual(np.average(mu_sample), 1.223012836345375) + self.assertAlmostEqual(np.average(sigma_sample), 6.992218962395242) - self.assertFalse(journal.number_of_simulations == 0) + self.assertFalse(journal.number_of_simulations[0] == 0) + self.assertEqual(journal.configuration["epsilon"], 10) + + def test_sample_simulation_budget(self): + # use the rejection sampling scheme with epsilon first + journal = self.sampler2.sample([self.y_obs], n_samples=None, simulation_budget=100, epsilon=20) + mu_sample = np.array(journal.get_parameters()['mu']) + sigma_sample = np.array(journal.get_parameters()['sigma']) + + mu_shape, sigma_shape = (len(mu_sample), mu_sample[0].shape[1]), \ + (len(sigma_sample), sigma_sample[0].shape[1]) + self.assertEqual(mu_shape, (3, 1)) + self.assertEqual(sigma_shape, (3, 1)) + + # Compute posterior mean + # self.assertAlmostEqual(np.average(np.asarray(samples[:,0])),1.22301,10e-2) + self.assertAlmostEqual(np.average(mu_sample), 0.8175361535037666) + self.assertAlmostEqual(np.average(sigma_sample), 8.155647092489977) + + self.assertEqual(journal.number_of_simulations[0], 100) + self.assertEqual(journal.configuration["epsilon"], 20) + + # use the rejection sampling scheme with the quantile + journal = self.sampler2.sample([self.y_obs], n_samples=None, simulation_budget=100, quantile=0.1) + mu_sample = np.array(journal.get_parameters()['mu']) + sigma_sample = np.array(journal.get_parameters()['sigma']) + + mu_shape, sigma_shape = (len(mu_sample), mu_sample[0].shape[1]), \ + (len(sigma_sample), sigma_sample[0].shape[1]) + self.assertEqual(mu_shape, (10, 1)) + self.assertEqual(sigma_shape, (10, 1)) + + # Compute posterior mean + # self.assertAlmostEqual(np.average(np.asarray(samples[:,0])),1.22301,10e-2) + self.assertAlmostEqual(np.average(mu_sample), 0.10394992719538543) + self.assertAlmostEqual(np.average(sigma_sample), 6.746940834914168) + + self.assertEqual(journal.number_of_simulations[0], 200) + + def test_errors(self): + with self.assertRaises(RuntimeError): + self.sampler.sample([self.y_obs], n_samples=10, simulation_budget=10) + with self.assertRaises(RuntimeError): + self.sampler.sample([self.y_obs], n_samples=10, quantile=0.1) + with self.assertRaises(RuntimeError): + self.sampler.sample([self.y_obs], n_samples=10, epsilon=None) + with self.assertRaises(RuntimeError): + self.sampler.sample([self.y_obs], n_samples=None, simulation_budget=100, quantile=0.1, epsilon=1) + with self.assertRaises(RuntimeError): + self.sampler.sample([self.y_obs], n_samples=None, simulation_budget=100, quantile=None, epsilon=None) class PMCTests(unittest.TestCase): @@ -65,7 +118,7 @@ def test_sample(self): self.model = Normal([mu, sigma]) # define sufficient statistics for the model - stat_calc = Identity(degree=2, cross=0) + stat_calc = Identity(degree=2, cross=False) # create fake observed data # y_obs = self.model.forward_simulate(1, np.random.RandomState(1))[0].tolist() @@ -77,7 +130,7 @@ def test_sample(self): T, n_sample, n_samples_per_param = 1, 10, 100 sampler = PMC([self.model], [likfun], backend, seed=1) journal = sampler.sample([y_obs], T, n_sample, n_samples_per_param, covFactors=np.array([.1, .1]), - iniPoints=None) + iniPoints=None, path_to_save_journal="tmp.jnl") mu_post_sample, sigma_post_sample, post_weights = np.array(journal.get_parameters()['mu']), np.array( journal.get_parameters()['sigma']), np.array(journal.get_weights()) @@ -91,8 +144,8 @@ def test_sample(self): self.assertEqual(mu_sample_shape, (10, 1)) self.assertEqual(sigma_sample_shape, (10, 1)) self.assertEqual(weights_sample_shape, (10, 1)) - self.assertLess(abs(mu_post_mean - (-3.373004641385251)), 1e-3) - self.assertLess(abs(sigma_post_mean - 6.519325027532673), 1e-3) + self.assertAlmostEqual(mu_post_mean, -3.3711206204663773, delta=1e-3) + self.assertAlmostEqual(sigma_post_mean, 6.519325027532673, delta=1e-3) self.assertFalse(journal.number_of_simulations == 0) @@ -114,12 +167,166 @@ def test_sample(self): self.assertEqual(mu_sample_shape, (10, 1)) self.assertEqual(sigma_sample_shape, (10, 1)) self.assertEqual(weights_sample_shape, (10, 1)) - self.assertLess(abs(mu_post_mean - (-3.2517600952705257)), 1e-3) - self.assertLess(abs(sigma_post_mean - 6.9214661382633365), 1e-3) - + self.assertAlmostEqual(mu_post_mean, -3.2517600952705257, delta=1e-3) + self.assertAlmostEqual(sigma_post_mean, 6.9214661382633365, delta=1e-3) self.assertFalse(journal.number_of_simulations == 0) +class MCMCMetropoliHastingsTests(unittest.TestCase): + # test if MCMCMetropoliHastings works + + def setUp(self): + # setup backend + self.backend = BackendDummy() + + # define a uniform prior distribution + mu = Uniform([[-5.0], [5.0]], name='mu') + sigma = Uniform([[0.0], [10.0]], name='sigma') + # define a Gaussian model + self.model = Normal([mu, sigma]) + self.model2 = Normal([mu, sigma]) + + # define sufficient statistics for the model + stat_calc = Identity(degree=2, cross=False) + + # create fake observed data + # y_obs = self.model.forward_simulate(1, np.random.RandomState(1))[0].tolist() + self.y_obs = [np.array(9.8)] + self.y_obs2 = [np.array(3.4)] + + # Define the likelihood function + self.likfun = SynLikelihood(stat_calc) + self.likfun2 = SynLikelihood(stat_calc) + + self.bounds = {"mu": [-5, 5], "sigma": (0, 10)} + + def test_sample(self): + n_sample, n_samples_per_param = 50, 20 + + sampler = MCMCMetropoliHastings([self.model], [self.likfun], self.backend, seed=1) + journal = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, + iniPoint=None, burnin=10, adapt_proposal_cov_interval=5, use_tqdm=False, + path_to_save_journal="tmp.jnl") + # without speedup_dummy + sampler = MCMCMetropoliHastings([self.model], [self.likfun], self.backend, seed=1) # to reset seed + journal_repr = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, use_tqdm=False, + iniPoint=None, speedup_dummy=False, burnin=10, adapt_proposal_cov_interval=5) + # Compute posterior mean + mu_post_mean_1, sigma_post_mean_1 = journal.posterior_mean()['mu'], journal.posterior_mean()['sigma'] + mu_post_mean_2, sigma_post_mean_2 = journal_repr.posterior_mean()['mu'], journal_repr.posterior_mean()['sigma'] + + self.assertAlmostEqual(mu_post_mean_1, 3.735372456561986) + self.assertAlmostEqual(mu_post_mean_2, -0.6946660151693353) + self.assertAlmostEqual(sigma_post_mean_1, 5.751158868437219) + self.assertAlmostEqual(sigma_post_mean_2, 8.103358539327967) + + def test_sample_with_transformer(self): + n_sample, n_samples_per_param = 50, 20 + + sampler = MCMCMetropoliHastings([self.model], [self.likfun], self.backend, seed=1) + journal = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, + iniPoint=None, burnin=10, adapt_proposal_cov_interval=5, use_tqdm=False, + bounds=self.bounds) + # Compute posterior mean + mu_post_mean, sigma_post_mean = journal.posterior_mean()['mu'], journal.posterior_mean()['sigma'] + + self.assertAlmostEqual(mu_post_mean, 1.3797371606192235) + self.assertAlmostEqual(sigma_post_mean, 8.097776586316062) + + # test raises correct errors: + with self.assertRaises(TypeError): + journal = sampler.sample([self.y_obs], n_sample, bounds=[0, 1]) + with self.assertRaises(KeyError): + journal = sampler.sample([self.y_obs], n_sample, bounds={"hello": (0, 1)}) + with self.assertRaises(RuntimeError): + journal = sampler.sample([self.y_obs], n_sample, bounds={"mu": (0)}) + with self.assertRaises(RuntimeError): + journal = sampler.sample([self.y_obs], n_sample, bounds={"mu": (0, 1, 2)}) + + def test_sample_two_models(self): + n_sample, n_samples_per_param = 50, 20 + + sampler = MCMCMetropoliHastings([self.model, self.model2], [self.likfun, self.likfun2], self.backend, + seed=1) + journal = sampler.sample([self.y_obs, self.y_obs2], n_sample, n_samples_per_param, cov_matrices=None, + iniPoint=None, burnin=10, adapt_proposal_cov_interval=5, use_tqdm=False) + # without speedup_dummy + sampler = MCMCMetropoliHastings([self.model, self.model2], [self.likfun, self.likfun2], self.backend, + seed=1) # to reset seed + journal_repr = sampler.sample([self.y_obs, self.y_obs2], n_sample, n_samples_per_param, cov_matrices=None, + iniPoint=None, speedup_dummy=False, burnin=10, adapt_proposal_cov_interval=5, + use_tqdm=False) + # Compute posterior mean + mu_post_mean_1, sigma_post_mean_1 = journal.posterior_mean()['mu'], journal.posterior_mean()['sigma'] + mu_post_mean_2, sigma_post_mean_2 = journal_repr.posterior_mean()['mu'], journal_repr.posterior_mean()['sigma'] + + self.assertAlmostEqual(mu_post_mean_1, 0.1920594166217264) + self.assertAlmostEqual(mu_post_mean_2, -1.0095854412936525) + self.assertAlmostEqual(sigma_post_mean_1, 9.143353645946233) + self.assertAlmostEqual(sigma_post_mean_2, 7.539268611159257) + + def test_restart_from_journal(self): + for speedup_dummy in [True, False]: + # do at once: + n_sample, n_samples_per_param = 40, 20 + sampler = MCMCMetropoliHastings([self.model], [self.likfun], self.backend, seed=1) + journal_at_once = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, + iniPoint=None, speedup_dummy=speedup_dummy, burnin=20, + adapt_proposal_cov_interval=10, use_tqdm=False) + # do separate: + n_sample, n_samples_per_param = 20, 20 + sampler = MCMCMetropoliHastings([self.model], [self.likfun], self.backend, seed=1) + journal = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, iniPoint=None, + speedup_dummy=speedup_dummy, burnin=20, adapt_proposal_cov_interval=10, + use_tqdm=False, path_to_save_journal="tmp.jnl") + + journal_separate = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, + iniPoint=None, journal_file="tmp.jnl", burnin=0, + speedup_dummy=speedup_dummy, use_tqdm=False) # restart from this journal + + self.assertEqual(journal_separate.configuration['n_samples'], journal_at_once.configuration['n_samples']) + self.assertEqual(journal_separate.number_of_simulations[-1], journal_at_once.number_of_simulations[-1]) + self.assertEqual(journal_separate.configuration["acceptance_rates"][-1], + journal_at_once.configuration["acceptance_rates"][-1]) + self.assertEqual(len(journal_separate.get_parameters()), len(journal_at_once.get_parameters())) + self.assertEqual(len(journal_separate.get_parameters()['mu']), len(journal_at_once.get_parameters()['mu'])) + self.assertEqual(len(journal_separate.get_accepted_parameters()), + len(journal_at_once.get_accepted_parameters())) + self.assertEqual(journal_separate.get_weights().shape, journal_at_once.get_weights().shape) + self.assertEqual(journal_separate.posterior_mean()['mu'], journal_at_once.posterior_mean()['mu']) + + def test_restart_from_journal_with_transformer(self): + for speedup_dummy in [True, False]: + # do at once: + n_sample, n_samples_per_param = 40, 20 + sampler = MCMCMetropoliHastings([self.model], [self.likfun], self.backend, seed=1) + journal_at_once = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, + iniPoint=None, speedup_dummy=speedup_dummy, burnin=20, + adapt_proposal_cov_interval=10, use_tqdm=False, bounds=self.bounds) + # do separate: + n_sample, n_samples_per_param = 20, 20 + sampler = MCMCMetropoliHastings([self.model], [self.likfun], self.backend, seed=1) + journal = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, iniPoint=None, + speedup_dummy=speedup_dummy, burnin=20, adapt_proposal_cov_interval=10, + use_tqdm=False, bounds=self.bounds, path_to_save_journal="tmp.jnl") + + journal_separate = sampler.sample([self.y_obs], n_sample, n_samples_per_param, cov_matrices=None, + iniPoint=None, journal_file="tmp.jnl", burnin=0, + speedup_dummy=speedup_dummy, use_tqdm=False, + bounds=self.bounds) # restart from this journal + + self.assertEqual(journal_separate.configuration['n_samples'], journal_at_once.configuration['n_samples']) + self.assertEqual(journal_separate.number_of_simulations[-1], journal_at_once.number_of_simulations[-1]) + self.assertEqual(journal_separate.configuration["acceptance_rates"][-1], + journal_at_once.configuration["acceptance_rates"][-1]) + self.assertEqual(len(journal_separate.get_parameters()), len(journal_at_once.get_parameters())) + self.assertEqual(len(journal_separate.get_parameters()['mu']), len(journal_at_once.get_parameters()['mu'])) + self.assertEqual(len(journal_separate.get_accepted_parameters()), + len(journal_at_once.get_accepted_parameters())) + self.assertEqual(journal_separate.get_weights().shape, journal_at_once.get_weights().shape) + self.assertAlmostEqual(journal_separate.posterior_mean()['mu'], journal_at_once.posterior_mean()['mu']) + + class PMCABCTests(unittest.TestCase): def setUp(self): # find spark and initialize it @@ -133,7 +340,7 @@ def setUp(self): self.model = Normal([mu, sigma]) # define a distance function - stat_calc = Identity(degree=2, cross=0) + stat_calc = Identity(degree=2, cross=False) self.dist_calc = Euclidean(stat_calc) # create fake observed data @@ -211,14 +418,39 @@ def test_sample(self): self.assertFalse(journal.number_of_simulations == 0) + # use the PMCABC scheme for T = 2 providing only first value for eps_arr + T, n_sample, n_simulate, eps_arr, eps_percentile = 2, 10, 1, [10], 10 + sampler = PMCABC([self.model], [self.dist_calc], self.backend, seed=1) + sampler.sample_from_prior(rng=np.random.RandomState(1)) + journal = sampler.sample([self.observation], T, eps_arr, n_sample, n_simulate, eps_percentile) + mu_post_sample, sigma_post_sample, post_weights = np.array(journal.get_parameters()['mu']), np.array( + journal.get_parameters()['sigma']), np.array(journal.get_weights()) + + # Compute posterior mean + mu_post_mean, sigma_post_mean = journal.posterior_mean()['mu'], journal.posterior_mean()['sigma'] + + # test shape of sample + mu_sample_shape, sigma_sample_shape, weights_sample_shape = (len(mu_post_sample), mu_post_sample[0].shape[1]), \ + (len(sigma_post_sample), + sigma_post_sample[0].shape[1]), post_weights.shape + + self.assertEqual(mu_sample_shape, (10, 1)) + self.assertEqual(sigma_sample_shape, (10, 1)) + self.assertEqual(weights_sample_shape, (10, 1)) + self.assertLess(mu_post_mean - 0.9356, 10e-2) + self.assertLess(sigma_post_mean - 7.819, 10e-2) + + self.assertFalse(journal.number_of_simulations == 0) + def test_restart_from_journal(self): # test with value of eps_arr_2 > percentile of distances n_sample, n_simulate, eps_arr, eps_percentile = 10, 1, [10, 5], 10 # 2 steps with intermediate journal: sampler = PMCABC([self.model], [self.dist_calc], self.backend, seed=1) sampler.sample_from_prior(rng=np.random.RandomState(1)) - journal_intermediate = sampler.sample([self.observation], 1, [eps_arr[0]], n_sample, n_simulate, eps_percentile) - journal_intermediate.save("tmp.jnl") + journal_intermediate = sampler.sample([self.observation], 1, [eps_arr[0]], n_sample, n_simulate, eps_percentile, + path_to_save_journal="tmp.jnl") + journal_final_1 = sampler.sample([self.observation], 1, [eps_arr[1]], n_sample, n_simulate, eps_percentile, journal_file="tmp.jnl") # 2 steps directly @@ -259,8 +491,9 @@ def setUp(self): self.model = Normal([mu, sigma]) # define a distance function - stat_calc = Identity(degree=2, cross=0) + stat_calc = Identity(degree=2, cross=False) self.dist_calc = Euclidean(stat_calc) + self.dist_calc_mmd = MMD(stat_calc, biased_estimator=False) # create fake observed data # self.observation = self.model.forward_simulate(1, np.random.RandomState(1))[0].tolist() @@ -289,7 +522,8 @@ def test_sample(self): # use the SABC scheme for T = 2 steps, epsilon, n_samples, n_samples_per_param = 2, 10, 10, 1 sampler = SABC([self.model], [self.dist_calc], self.backend, seed=1) - journal = sampler.sample([self.observation], steps, epsilon, n_samples, n_samples_per_param) + journal = sampler.sample([self.observation], steps, epsilon, n_samples, n_samples_per_param, full_output=1, + path_to_save_journal="tmp") mu_post_sample, sigma_post_sample, post_weights = np.array(journal.get_parameters()['mu']), np.array( journal.get_parameters()['sigma']), np.array(journal.get_weights()) @@ -309,6 +543,10 @@ def test_sample(self): self.assertFalse(journal.number_of_simulations == 0) + # check whether it raises the correct error with MMD: + with self.assertRaises(RuntimeError): + sampler = SABC([self.model], [self.dist_calc_mmd], self.backend, seed=1) + class ABCsubsimTests(unittest.TestCase): def setUp(self): @@ -322,7 +560,7 @@ def setUp(self): self.model = Normal([mu, sigma]) # define a distance function - stat_calc = Identity(degree=2, cross=0) + stat_calc = Identity(degree=2, cross=False) self.dist_calc = Euclidean(stat_calc) # create fake observed data @@ -386,14 +624,23 @@ def setUp(self): self.model = Normal([mu, sigma]) # define a distance function - stat_calc = Identity(degree=2, cross=0) + stat_calc = Identity(degree=2, cross=False) self.dist_calc = Euclidean(stat_calc) # create fake observed data - # self.observation = self.model.forward_simulate(1, np.random.RandomState(1))[0].tolist() self.observation = [np.array(9.8)] - def test_sample(self): + # set up for the Bernton et al. implementation: + # define a distance function + stat_calc_2 = Identity(degree=1, cross=False) + self.dist_calc_2 = MMD(stat_calc_2) + # create fake observed data + seed = 12 + self.rng = np.random.RandomState(seed) + self.observation_2 = self.rng.normal(loc=0, scale=1, size=10) + self.observation_2 = [x for x in self.observation_2] + + def test_sample_delmoral(self): # use the SMCABC scheme for T = 1 steps, n_sample, n_simulate = 1, 10, 1 sampler = SMCABC([self.model], [self.dist_calc], self.backend, seed=1) @@ -432,11 +679,159 @@ def test_sample(self): self.assertEqual(mu_sample_shape, (10, 1)) self.assertEqual(sigma_sample_shape, (10, 1)) self.assertEqual(weights_sample_shape, (10, 1)) - self.assertLess(mu_post_mean - (-0.786118677019), 10e-2) - self.assertLess(sigma_post_mean - 4.63324738665, 10e-2) + self.assertAlmostEqual(mu_post_mean, - 0.8888295384029634, delta=10e-3) + self.assertAlmostEqual(sigma_post_mean, 4.299346466029422, delta=10e-3) + + self.assertEqual(journal.number_of_simulations[-1], 19) + + # try now with the r-hit kernel version 1: + T, n_sample, n_simulate = 2, 10, 1 + sampler = SMCABC([self.model], [self.dist_calc], self.backend, seed=1) + journal = sampler.sample([self.observation], T, n_sample, n_simulate, which_mcmc_kernel=1) + mu_post_sample, sigma_post_sample, post_weights = np.array(journal.get_parameters()['mu']), np.array( + journal.get_parameters()['sigma']), np.array(journal.get_weights()) + + # Compute posterior mean + mu_post_mean, sigma_post_mean = journal.posterior_mean()['mu'], journal.posterior_mean()['sigma'] + + # test shape of sample + mu_sample_shape, sigma_sample_shape, weights_sample_shape = (len(mu_post_sample), mu_post_sample[0].shape[1]), \ + (len(sigma_post_sample), + sigma_post_sample[0].shape[1]), post_weights.shape + self.assertEqual(mu_sample_shape, (10, 1)) + self.assertEqual(sigma_sample_shape, (10, 1)) + self.assertEqual(weights_sample_shape, (10, 1)) + self.assertAlmostEqual(mu_post_mean, -0.6507386970288184, delta=10e-3) + self.assertAlmostEqual(sigma_post_mean, 6.253446572247367, delta=10e-3) + + self.assertEqual(journal.number_of_simulations[-1], 56) + + # try now with the r-hit kernel version 2: + T, n_sample, n_simulate = 2, 10, 1 + sampler = SMCABC([self.model], [self.dist_calc], self.backend, seed=1) + journal = sampler.sample([self.observation], T, n_sample, n_simulate, which_mcmc_kernel=2) + mu_post_sample, sigma_post_sample, post_weights = np.array(journal.get_parameters()['mu']), np.array( + journal.get_parameters()['sigma']), np.array(journal.get_weights()) + + # Compute posterior mean + mu_post_mean, sigma_post_mean = journal.posterior_mean()['mu'], journal.posterior_mean()['sigma'] + + # test shape of sample + mu_sample_shape, sigma_sample_shape, weights_sample_shape = (len(mu_post_sample), mu_post_sample[0].shape[1]), \ + (len(sigma_post_sample), + sigma_post_sample[0].shape[1]), post_weights.shape + self.assertEqual(mu_sample_shape, (10, 1)) + self.assertEqual(sigma_sample_shape, (10, 1)) + self.assertEqual(weights_sample_shape, (10, 1)) + self.assertAlmostEqual(mu_post_mean, -0.5486451602421536, delta=10e-3) + self.assertAlmostEqual(sigma_post_mean, 3.633148439032683, delta=10e-3) self.assertFalse(journal.number_of_simulations == 0) + def test_sample_bernton(self): + # check using the standard MCMC kernel raises error: + sampler = SMCABC([self.model], [self.dist_calc_2], self.backend, seed=1, version="Bernton") + with self.assertRaises(RuntimeError): + journal = sampler.sample([self.observation_2], 1, 10, 10, which_mcmc_kernel=0, alpha=0.5, + epsilon_final=0) + + # try now with the r-hit kernel version 1: + T, n_sample, n_simulate = 3, 10, 10 + sampler = SMCABC([self.model], [self.dist_calc_2], self.backend, seed=1, version="Bernton") + journal = sampler.sample([self.observation_2], T, n_sample, n_simulate, which_mcmc_kernel=1, alpha=0.5, + epsilon_final=0) + mu_post_sample, sigma_post_sample, post_weights = np.array(journal.get_parameters()['mu']), np.array( + journal.get_parameters()['sigma']), np.array(journal.get_weights()) + + # Compute posterior mean + mu_post_mean, sigma_post_mean = journal.posterior_mean()['mu'], journal.posterior_mean()['sigma'] + + # test shape of sample + mu_sample_shape, sigma_sample_shape, weights_sample_shape = (len(mu_post_sample), mu_post_sample[0].shape[1]), \ + (len(sigma_post_sample), + sigma_post_sample[0].shape[1]), post_weights.shape + self.assertEqual(mu_sample_shape, (10, 1)) + self.assertEqual(sigma_sample_shape, (10, 1)) + self.assertEqual(weights_sample_shape, (10, 1)) + self.assertAlmostEqual(mu_post_mean, -0.7294075767448996, delta=10e-3) + self.assertAlmostEqual(sigma_post_mean, 3.406347345226374, delta=10e-3) + + self.assertEqual(journal.number_of_simulations[-1], 2286) + + # try now with the r-hit kernel version 2: + T, n_sample, n_simulate = 3, 10, 10 + sampler = SMCABC([self.model], [self.dist_calc_2], self.backend, seed=1, version="Bernton") + journal = sampler.sample([self.observation_2], T, n_sample, n_simulate, which_mcmc_kernel=2, + epsilon_final=0.1) + mu_post_sample, sigma_post_sample, post_weights = np.array(journal.get_parameters()['mu']), np.array( + journal.get_parameters()['sigma']), np.array(journal.get_weights()) + + # Compute posterior mean + mu_post_mean, sigma_post_mean = journal.posterior_mean()['mu'], journal.posterior_mean()['sigma'] + + # test shape of sample + mu_sample_shape, sigma_sample_shape, weights_sample_shape = (len(mu_post_sample), mu_post_sample[0].shape[1]), \ + (len(sigma_post_sample), + sigma_post_sample[0].shape[1]), post_weights.shape + self.assertEqual(mu_sample_shape, (10, 1)) + self.assertEqual(sigma_sample_shape, (10, 1)) + self.assertEqual(weights_sample_shape, (10, 1)) + self.assertAlmostEqual(mu_post_mean, -2.1412732987491303, delta=10e-3) + self.assertAlmostEqual(sigma_post_mean, 5.146988585331478, delta=10e-3) + + self.assertEqual(journal.number_of_simulations[-1], 127) + + def test_restart_from_journal_delmoral(self): + n_sample, n_simulate = 10, 1 + # loop over standard MCMC kernel, r-hit kernel version 1 and r-hit kernel version 2 + for which_mcmc_kernel in [0, 1, 2]: + # 2 steps with intermediate journal: + sampler = SMCABC([self.model], [self.dist_calc], self.backend, seed=1) + journal_intermediate = sampler.sample([self.observation], 2, n_sample, n_simulate, + which_mcmc_kernel=which_mcmc_kernel) + journal_intermediate.save("tmp.jnl") + journal_final_1 = sampler.sample([self.observation], 1, n_sample, n_simulate, + which_mcmc_kernel=which_mcmc_kernel, + journal_file="tmp.jnl") + + # 2 steps directly + sampler = SMCABC([self.model], [self.dist_calc], self.backend, seed=1) + journal_final_2 = sampler.sample([self.observation], 3, n_sample, n_simulate, full_output=1, + which_mcmc_kernel=which_mcmc_kernel) + self.assertEqual(journal_final_1.configuration["epsilon_arr"], journal_final_2.configuration["epsilon_arr"]) + self.assertEqual(journal_final_1.posterior_mean()['mu'], journal_final_2.posterior_mean()['mu']) + + def test_restart_from_journal_bernton(self): + n_sample, n_simulate = 10, 10 + # loop over standard MCMC kernel, r-hit kernel version 1 and r-hit kernel version 2 + for which_mcmc_kernel in [1, 2]: + # 2 steps with intermediate journal: + sampler = SMCABC([self.model], [self.dist_calc_2], self.backend, seed=1, version="Bernton") + journal_intermediate = sampler.sample([self.observation_2], 1, n_sample, n_simulate, + which_mcmc_kernel=which_mcmc_kernel) + journal_intermediate.save("tmp.jnl") + journal_final_1 = sampler.sample([self.observation_2], 1, n_sample, n_simulate, + which_mcmc_kernel=which_mcmc_kernel, + journal_file="tmp.jnl") + + # 2 steps directly + sampler = SMCABC([self.model], [self.dist_calc_2], self.backend, seed=1, version="Bernton") + journal_final_2 = sampler.sample([self.observation_2], 2, n_sample, n_simulate, full_output=1, + which_mcmc_kernel=which_mcmc_kernel) + self.assertEqual(journal_final_1.configuration["epsilon_arr"], journal_final_2.configuration["epsilon_arr"]) + self.assertEqual(journal_final_1.posterior_mean()['mu'], journal_final_2.posterior_mean()['mu']) + + def test_errors(self): + with self.assertRaises(RuntimeError): + sampler = SMCABC([self.model], [self.dist_calc_2], self.backend, seed=1, version="DelMoral") + with self.assertRaises(RuntimeError): + sampler = SMCABC([self.model], [self.dist_calc], self.backend, seed=1, version="Ciao") + + sampler = SMCABC([self.model], [self.dist_calc], self.backend, seed=1, version="Bernton") + with self.assertRaises(RuntimeError): + journal = sampler.sample([self.observation], 1, 10, 10, which_mcmc_kernel=4, alpha=0.5, + epsilon_final=0) + class APMCABCTests(unittest.TestCase): def setUp(self): @@ -450,7 +845,7 @@ def setUp(self): self.model = Normal([mu, sigma]) # define a distance function - stat_calc = Identity(degree=2, cross=0) + stat_calc = Identity(degree=2, cross=False) self.dist_calc = Euclidean(stat_calc) # create fake observed data @@ -513,7 +908,7 @@ def setUp(self): self.model = Normal([mu, sigma]) # define a distance function - stat_calc = Identity(degree=2, cross=0) + stat_calc = Identity(degree=2, cross=False) self.dist_calc = Euclidean(stat_calc) # create fake observed data diff --git a/tests/output_tests.py b/tests/output_tests.py index 12e7fe14..7ebfd882 100644 --- a/tests/output_tests.py +++ b/tests/output_tests.py @@ -135,7 +135,7 @@ def test_plot_wass_dist(self): journal_3.add_weights(weights_identical) self.assertRaises(RuntimeError, journal_3.Wass_convergence_plot) journal_4 = Journal(1) - journal_4.add_accepted_parameters(np.array([np.array([1]), np.array([1, 2])])) + journal_4.add_accepted_parameters(np.array([np.array([1]), np.array([1, 2])], dtype="object")) print(len(journal_4.accepted_parameters)) self.assertRaises(RuntimeError, journal_4.Wass_convergence_plot) @@ -169,6 +169,21 @@ def test_plot_post_distr(self): with self.assertRaises(TypeError): journal.plot_posterior_distr(ranges_parameters={"par1": np.zeros(1)}) + def test_traceplot(self): + rng = np.random.RandomState(1) + weights_identical = np.ones((100, 1)) + params = rng.randn(100).reshape(-1, 1) + journal = Journal(1) + journal.add_weights(weights_identical) + journal.add_accepted_parameters(params) + journal.add_user_parameters([("mu", params[:, 0])]) + self.assertRaises(RuntimeError, journal.traceplot) # as it does not have "acceptance_rates" in configuration + journal.configuration["acceptance_rates"] = [0.3] + with self.assertRaises(KeyError): + journal.traceplot(parameters_to_show=["sigma"]) + # now try correctly: + fig, ax = journal.traceplot() + if __name__ == '__main__': unittest.main() diff --git a/tests/statisticslearning_tests.py b/tests/statisticslearning_tests.py index 77823cdb..8247bd40 100644 --- a/tests/statisticslearning_tests.py +++ b/tests/statisticslearning_tests.py @@ -65,11 +65,12 @@ def setUp(self): if has_torch: # Initialize statistics learning self.statisticslearning = SemiautomaticNN([self.Y], self.statistics_cal, self.backend, n_samples=100, - n_samples_per_param=1, seed=1, n_epochs=10, scale_samples=False) + n_samples_per_param=1, seed=1, n_epochs=10, scale_samples=False, + use_tqdm=False) # with sample scaler: self.statisticslearning_with_scaler = SemiautomaticNN([self.Y], self.statistics_cal, self.backend, n_samples=100, n_samples_per_param=1, seed=1, - n_epochs=10, scale_samples=True) + n_epochs=10, scale_samples=True, use_tqdm=False) def test_initialization(self): if not has_torch: @@ -164,12 +165,12 @@ def setUp(self): # Initialize statistics learning self.statisticslearning = ContrastiveDistanceLearning([self.Y], self.statistics_cal, self.backend, n_samples=100, n_samples_per_param=1, seed=1, - n_epochs=10, scale_samples=False) + n_epochs=10, scale_samples=False, use_tqdm=False) # with sample scaler: self.statisticslearning_with_scaler = ContrastiveDistanceLearning([self.Y], self.statistics_cal, self.backend, n_samples=100, n_samples_per_param=1, seed=1, - n_epochs=10, scale_samples=True) + n_epochs=10, scale_samples=True, use_tqdm=False) def test_initialization(self): if not has_torch: @@ -212,11 +213,11 @@ def setUp(self): if has_torch: # Initialize statistics learning self.statisticslearning = TripletDistanceLearning([self.Y], self.statistics_cal, self.backend, - scale_samples=False, + scale_samples=False, use_tqdm=False, n_samples=100, n_samples_per_param=1, seed=1, n_epochs=10) # with sample scaler: self.statisticslearning_with_scaler = TripletDistanceLearning([self.Y], self.statistics_cal, self.backend, - scale_samples=True, + scale_samples=True, use_tqdm=False, n_samples=100, n_samples_per_param=1, seed=1, n_epochs=10) diff --git a/tests/test_examples.py b/tests/test_examples.py index c5658d39..2fb97d6d 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -14,7 +14,7 @@ def test_pmc(self): from examples.approx_lhd.pmc_hierarchical_models import infer_parameters journal = infer_parameters(steps=1, n_sample=50) test_result = journal.posterior_mean()["school_location"] - expected_result = 0.2566394909510058 + expected_result = 0.2569671145768137 self.assertAlmostEqual(test_result, expected_result) diff --git a/tests/test_examples_mpi.py b/tests/test_examples_mpi.py index 1f148ac6..e7f6af83 100644 --- a/tests/test_examples_mpi.py +++ b/tests/test_examples_mpi.py @@ -27,7 +27,7 @@ def test_example(self): from examples.backends.mpi.pmcabc_gaussian import infer_parameters journal = infer_parameters(backend_mpi, steps=3, n_sample=50) test_result = journal.posterior_mean()['mu'] - expected_result = 174.94735914414332 + expected_result = 174.94717012502286 self.assertAlmostEqual(test_result, expected_result) diff --git a/tests/transformers_tests.py b/tests/transformers_tests.py new file mode 100644 index 00000000..9dfa3d66 --- /dev/null +++ b/tests/transformers_tests.py @@ -0,0 +1,57 @@ +import unittest + +import numpy as np + +from abcpy.transformers import DummyTransformer, BoundedVarTransformer + + +class DummyTransformerTests(unittest.TestCase): + def test(self): + transformer = DummyTransformer() + x = [np.array([3.2]), np.array([2.4])] + self.assertEqual(x, transformer.transform(x)) + self.assertEqual(x, transformer.inverse_transform(x)) + self.assertEqual(0, transformer.jac_log_det_inverse_transform(x)) + self.assertEqual(0, transformer.jac_log_det(x)) + + +class BoundedVarTransformerTests(unittest.TestCase): + def setUp(self): + self.transformer_lower_bounded = BoundedVarTransformer(lower_bound=np.array([0, 0]), + upper_bound=np.array([None, None])) + self.transformer_two_sided = BoundedVarTransformer(lower_bound=np.array([0, 0]), upper_bound=np.array([10, 10])) + self.transformer_mixed = BoundedVarTransformer(lower_bound=np.array([0, 0]), upper_bound=np.array([10, None])) + self.transformer_dummy = BoundedVarTransformer(lower_bound=np.array([None, None]), + upper_bound=np.array([None, None])) + self.list_transformers = [self.transformer_dummy, self.transformer_mixed, + self.transformer_two_sided, self.transformer_lower_bounded] + + def test(self): + x = [np.array([3.2]), np.array([2.4])] + for transformer in self.list_transformers: + self.assertEqual(len(x), len(transformer.inverse_transform(transformer.transform(x)))) + self.assertTrue(np.allclose(np.array(x), np.array(transformer.inverse_transform(transformer.transform(x))))) + self.assertAlmostEqual(transformer.jac_log_det(x), + transformer.jac_log_det_inverse_transform(transformer.transform(x)), delta=1e-7) + + # test dummy transformer actually does nothing: + self.assertEqual(x, self.transformer_dummy.transform(x)) + self.assertEqual(x, self.transformer_dummy.inverse_transform(x)) + self.assertEqual(0, self.transformer_dummy.jac_log_det_inverse_transform(x)) + self.assertEqual(0, self.transformer_dummy.jac_log_det(x)) + + def test_errors(self): + with self.assertRaises(RuntimeError): + transformer = BoundedVarTransformer(lower_bound=[0, 0], upper_bound=[10, 10]) + with self.assertRaises(RuntimeError): + transformer = BoundedVarTransformer(lower_bound=np.array([0, 0]), upper_bound=np.array([10])) + with self.assertRaises(NotImplementedError): + transformer = BoundedVarTransformer(lower_bound=np.array([None, 0]), upper_bound=np.array([10, 10])) + with self.assertRaises(RuntimeError): + self.transformer_lower_bounded.transform(x=[np.array([3.2]), np.array([-2.4])]) + with self.assertRaises(RuntimeError): + self.transformer_two_sided.transform(x=[np.array([13.2]), np.array([2.4])]) + + +if __name__ == '__main__': + unittest.main()