Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"""
if self.accepted_weights_bds is None:
return 1.0 / self.n_samples
else:
prior_prob = self.pdf_of_prior(self.model, theta, 0)[0][0] #first [0] because it returns a list, second [0] because if we have multiple models
denominator = 0.0
for i in range(0, self.n_samples):
pdf_value = self.kernel.pdf(self.accepted_parameters_bds.value()[i,:], self.accepted_cov_mat_bds.value(), theta)
denominator += self.accepted_weights_bds.value()[i, 0] * pdf_value
return 1.0 * prior_prob / denominator
class PMC(BasePMC, 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
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
inference schemes are based on Andrieu and Roberts [2].
[1] Cappé, O., Guillin, A., Marin, J.-M., and Robert, C. P. (2004). Population Monte Carlo.
Journal of Computational and Graphical Statistics, 13(4), 907–929.
[2] C. Andrieu and G. O. Roberts. The pseudo-marginal approach for efficient Monte Carlo computations.
Annals of Statistics, 37(2):697–725, 04 2009.
"""
raise NotImplementedError
@abstractproperty
def n_samples_per_param(self):
"""To be overwritten by any sub-class: an attribute specifying the number of data points in each simulated data set."""
raise NotImplementedError
@abstractproperty
def observations_bds(self):
"""To be overwritten by any sub-class: an attribute saving the observations as bds
"""
raise NotImplementedError
class BasePMC(InferenceMethod, metaclass = ABCMeta):
"""
This abstract base class represents inference methods that use Population Monte Carlo.
"""
@abstractmethod
def _update_broadcasts(self, accepted_parameters, accepted_weights, accepted_cov_mat):
"""
To be overwritten by any sub-class: broadcasts visited values
Parameters
----------
accepted_parameters: numpy.array
Contains all new accepted parameters.
accepted_weights: numpy.array
Contains all the new accepted weights.
accepted_cov_mat: numpy.ndarray
y_sim = self.model[0].sample_from_distribution(self.n_samples_per_param, rng=rng).tolist()
distance = self.distance.distance(self.observations_bds.value(), y_sim)
ratio_prior_prob = self.pdf_of_prior(self.model, new_theta, 0)[0][0] / self.pdf_of_prior(self.model, theta, 0)[0][0]
kernel_numerator = self.kernel.pdf(new_theta, self.accepted_cov_mat_bds.value(), theta)
kernel_denominator = self.kernel.pdf(theta, self.accepted_cov_mat_bds.value(), new_theta)
ratio_kernel_prob = kernel_numerator / kernel_denominator
probability_acceptance = min(1, ratio_prior_prob * ratio_kernel_prob)
if distance < self.epsilon[-1] and rng.binomial(1, probability_acceptance) == 1:
index_accept += 1
else:
self.set_parameters(self.model, theta, 0)
distance = self.accepted_dist_bds.value()[index[0]]
return (self.get_parameters(self.model), distance, index_accept)
class APMCABC(BaseAdaptivePopulationMC, InferenceMethod):
"""This base class implements Adaptive Population Monte Carlo Approximate Bayesian computation of
M. Lenormand et al. [1].
[1] M. Lenormand, F. Jabot and G. Deffuant, Adaptive approximate Bayesian computation
for complex models. Computational Statistics, 28:2777–2796, 2013.
Parameters
----------
model : abcpy.models.Model
Model object defining the model to be used.
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.
backend : abcpy.backends.Backend
Backend object defining the backend to be used.
## Calculate acceptance probability:
ratio_prior_prob = self.pdf_of_prior(self.model, new_theta, 0)[0][0] / self.pdf_of_prior(self.model,
self.accepted_parameters_bds.value()[index, :], 0)[0][0]
ratio_likelihood_prob = np.exp((self.smooth_distances_bds.value()[index] - smooth_distance) / self.epsilon)
acceptance_prob = ratio_prior_prob * ratio_likelihood_prob
## If accepted
if rng.rand(1) < acceptance_prob:
acceptance = 1
else:
distance = np.inf
return (new_theta, distance, all_parameters, all_distances, index, acceptance)
class ABCsubsim(BaseAnnealing, InferenceMethod):
"""This base class implements Approximate Bayesian Computation by subset simulation (ABCsubsim) algorithm of [1].
[1] M. Chiachio, J. L. Beck, J. Chiachio, and G. Rus., Approximate Bayesian computation by subset
simulation. SIAM J. Sci. Comput., 36(3):A1339–A1358, 2014/10/03 2014.
Parameters
----------
model : abcpy.models.Model
Model object defining the model to be used.
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.
backend : abcpy.backends.Backend
Backend object defining the backend to be used.
seed : integer, optional
"""
raise NotImplementedError
@abstractproperty
def accepted_parameters_bds(self):
"""To be overwritten by any sub-class: an attribute saving the accepted parameters as bds
"""
raise NotImplementedError
@abstractproperty
def accepted_cov_mat_bds(self):
"""To be overwritten by any sub-class: an attribute saving the accepted covariance matrix as bds
"""
raise NotImplementedError
class RejectionABC(InferenceMethod):
"""This base class implements the rejection algorithm based inference scheme [1] for
Approximate Bayesian Computation.
[1] Tavaré, S., Balding, D., Griffith, R., Donnelly, P.: Inferring coalescence
times from DNA sequence data. Genetics 145(2), 505–518 (1997).
Parameters
----------
model: abcpy.models.Model
Model object defining the model to be used.
distance: abcpy.distances.Distance
Distance object defining the distance measure to compare simulated and observed data sets.
backend: abcpy.backends.Backend
Backend object defining the backend to be used.
seed: integer, optional
Optional initial seed for the random number generator. The default value is generated randomly.
if self.accepted_weights_bds is None:
return 1.0 / self.n_samples
else:
#TODO MULTIPLE MODELS
prior_prob = self.pdf_of_prior(self.model, theta, 0)[0][0]
denominator = 0.0
for i in range(0, self.n_samples):
pdf_value = self.kernel.pdf(self.accepted_parameters_bds.value()[i,:], self.accepted_cov_mat_bds.value(), theta)
denominator += self.accepted_weights_bds.value()[i, 0] * pdf_value
return 1.0 * prior_prob / denominator
class SABC(BaseAnnealing, InferenceMethod):
"""
This base class implements a modified version of Simulated Annealing Approximate Bayesian Computation (SABC) of [1] when the prior is non-informative.
[1] C. Albert, H. R. Kuensch and A. Scheidegger. A Simulated Annealing Approach to
Approximate Bayes Computations. Statistics and Computing, (2014).
Parameters
----------
model : abcpy.models.Model
Model object defining the model to be used.
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.
backend : abcpy.backends.Backend
Backend object defining the backend to be used.
"""
raise NotImplementedError
@abstractproperty
def accepted_parameters_bds(self):
"""To be overwritten by any sub-class: an attribute saving the accepted parameters as bds
"""
raise NotImplementedError
@abstractproperty
def accepted_cov_mat_bds(self):
"""To be overwritten by any sub-class: an attribute saving the accepted covariance matrix as bds
"""
raise NotImplementedError
class BaseAdaptivePopulationMC(InferenceMethod, metaclass = ABCMeta):
"""
This abstract base class represents inference methods that use Adaptive Population Monte Carlo.
"""
@abstractmethod
def _update_broadcasts(self):
"""
To be overwritten by any sub-class: broadcasts visited values
Parameters
----------
accepted_parameters: numpy.array
Contains all new accepted parameters.
accepted_weights: numpy.array
Contains all the new accepted weights.
kernel_numerator = self.kernel.pdf(new_theta, accepted_cov_mat_transformed, theta)
kernel_denominator = self.kernel.pdf(theta, accepted_cov_mat_transformed, new_theta)
ratio_likelihood_prob = kernel_numerator / kernel_denominator
acceptance_prob = min(1, ratio_prior_prob * ratio_likelihood_prob) * (
new_distance < self.anneal_parameter)
## If accepted
if rng.binomial(1, acceptance_prob) == 1:
theta = new_theta
acceptance = acceptance + 1
if acceptance / 10 <= 0.5 and acceptance / 10 >= 0.3:
return (accepted_cov_mat_transformed, t, 1)
else:
return (accepted_cov_mat_transformed, t, 0)
#NOTE when testing with example values -> raises singular matrix error during calculation of pdf of kernel. I think this might happen in general, not a mistake of the code ---> desired behavior?
class RSMCABC(BaseAdaptivePopulationMC, InferenceMethod):
"""This base class implements Adaptive Population Monte Carlo Approximate Bayesian computation of
Drovandi and Pettitt [1].
[1] CC. Drovandi CC and AN. Pettitt, Estimation of parameters for macroparasite population evolution using
approximate Bayesian computation. Biometrics 67(1):225–233, 2011.
Parameters
----------
model : abcpy.models.Model
Model object defining the model to be used.
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.
backend : abcpy.backends.Backend
Backend object defining the backend to be used.
#todo multiple models
y_sim = self.model[0].sample_from_distribution(self.n_samples_per_param, rng=rng).tolist()
dist = self.distance.distance(self.observations_bds.value(), y_sim)
prior_prob = self.pdf_of_prior(self.model, new_theta, 0)[0][0]
denominator = 0.0
for i in range(0, len(self.accepted_weights_bds.value())):
pdf_value = self.kernel.pdf(self.accepted_parameters_bds.value()[i,:], self.accepted_cov_mat_bds.value(), new_theta)
denominator += self.accepted_weights_bds.value()[i, 0] * pdf_value
weight = 1.0 * prior_prob / denominator
return (self.get_parameters(self.model), dist, weight)
#NOTE takes long time to test with example, but no obvious mistakes so far
class SMCABC(BaseAdaptivePopulationMC, InferenceMethod):
"""This base class implements Adaptive Population Monte Carlo Approximate Bayesian computation of
Del Moral et al. [1].
[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.
Parameters
----------
model : abcpy.models.Model
Model object defining the model to be used.
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.
backend : abcpy.backends.Backend
Backend object defining the backend to be used.
@abstractproperty
def accepted_weights_bds(self):
"""To be overwritten by any sub-class: an attribute saving the accepted weights as bds
"""
raise NotImplementedError
@abstractproperty
def accepted_cov_mat_bds(self):
"""To be overwritten by any sub-class: an attribute saving the accepted covariance matrix as bds
"""
raise NotImplementedError
class BaseAnnealing(InferenceMethod, metaclass = ABCMeta):
"""
This abstract base class represents inference methods that use annealing.
"""
@abstractmethod
def _update_broadcasts(self):
raise NotImplementedError
@abstractmethod
def _accept_parameter(self):
raise NotImplementedError
@abstractproperty
def distance(self):
"""To be overwritten by any sub-class: an attribute specifying the distance measure to be used