diff --git a/examples/other_examples/linear_regression_pymc3.py b/examples/other_examples/linear_regression_pymc3.py
new file mode 100644
index 0000000000000000000000000000000000000000..ccd4f48eb40d2416e273dd7f75536abbeb49495d
--- /dev/null
+++ b/examples/other_examples/linear_regression_pymc3.py
@@ -0,0 +1,63 @@
+#!/bin/python
+"""
+An example of how to use tupak to perform paramater estimation for
+non-gravitational wave data. In this case, fitting a linear function to
+data with background Gaussian noise
+
+"""
+from __future__ import division
+import tupak
+import numpy as np
+import matplotlib.pyplot as plt
+import inspect
+
+from tupak.core.likelihood import GaussianLikelihood
+
+# A few simple setup steps
+label = 'linear_regression_pymc3'
+outdir = 'outdir'
+tupak.utils.check_directory_exists_and_if_not_mkdir(outdir)
+
+# First, we define our "signal model", in this case a simple linear function
+def model(time, m, c):
+ return time * m + c
+
+# Now we define the injection parameters which we make simulated data with
+injection_parameters = dict(m=0.5, c=0.2)
+
+# For this example, we'll use standard Gaussian noise
+sigma = 1
+
+# These lines of code generate the fake data. Note the ** just unpacks the
+# contents of the injection_paramsters when calling the model function.
+sampling_frequency = 10
+time_duration = 10
+time = np.arange(0, time_duration, 1/sampling_frequency)
+N = len(time)
+data = model(time, **injection_parameters) + np.random.normal(0, sigma, N)
+
+# We quickly plot the data to check it looks sensible
+fig, ax = plt.subplots()
+ax.plot(time, data, 'o', label='data')
+ax.plot(time, model(time, **injection_parameters), '--r', label='signal')
+ax.set_xlabel('time')
+ax.set_ylabel('y')
+ax.legend()
+fig.savefig('{}/{}_data.png'.format(outdir, label))
+
+# Now lets instantiate a version of our GaussianLikelihood, giving it
+# the time, data and signal model
+likelihood = GaussianLikelihood(time, data, model, sigma=sigma)
+
+# From hereon, the syntax is exactly equivalent to other tupak examples
+# We make a prior
+priors = {}
+priors['m'] = tupak.core.prior.Uniform(0, 5, 'm')
+priors['c'] = tupak.core.prior.Uniform(-2, 2, 'c')
+
+# And run sampler
+result = tupak.run_sampler(
+ likelihood=likelihood, priors=priors, sampler='pymc3',
+ injection_parameters=injection_parameters, outdir=outdir,
+ draws=2000, label=label)
+result.plot_corner()
diff --git a/examples/other_examples/linear_regression_pymc3_custom_likelihood.py b/examples/other_examples/linear_regression_pymc3_custom_likelihood.py
new file mode 100644
index 0000000000000000000000000000000000000000..3b3ab8af138ead0dc1428e1737eb46ed1ea255ca
--- /dev/null
+++ b/examples/other_examples/linear_regression_pymc3_custom_likelihood.py
@@ -0,0 +1,149 @@
+#!/bin/python
+"""
+An example of how to use tupak to perform paramater estimation for
+non-gravitational wave data. In this case, fitting a linear function to
+data with background Gaussian noise. This example uses a custom
+likelihood function to show how it should be defined, although this
+would give equivalent results as using the pre-defined 'Gaussian Likelihood'
+
+"""
+
+from __future__ import division
+import tupak
+import numpy as np
+import matplotlib.pyplot as plt
+import inspect
+import pymc3 as pm
+
+# A few simple setup steps
+label = 'linear_regression_pymc3_custom_likelihood'
+outdir = 'outdir'
+tupak.utils.check_directory_exists_and_if_not_mkdir(outdir)
+
+# First, we define our "signal model", in this case a simple linear function
+def model(time, m, c):
+ return time * m + c
+
+# Now we define the injection parameters which we make simulated data with
+injection_parameters = dict(m=0.5, c=0.2)
+
+# For this example, we'll use standard Gaussian noise
+sigma = 1
+
+# These lines of code generate the fake data. Note the ** just unpacks the
+# contents of the injection_paramsters when calling the model function.
+sampling_frequency = 10
+time_duration = 10
+time = np.arange(0, time_duration, 1/sampling_frequency)
+N = len(time)
+data = model(time, **injection_parameters) + np.random.normal(0, sigma, N)
+
+# We quickly plot the data to check it looks sensible
+fig, ax = plt.subplots()
+ax.plot(time, data, 'o', label='data')
+ax.plot(time, model(time, **injection_parameters), '--r', label='signal')
+ax.set_xlabel('time')
+ax.set_ylabel('y')
+ax.legend()
+fig.savefig('{}/{}_data.png'.format(outdir, label))
+
+
+# Parameter estimation: we now define a Gaussian Likelihood class relevant for
+# our model.
+class GaussianLikelihoodPyMC3(tupak.Likelihood):
+ def __init__(self, x, y, sigma, function):
+ """
+ A general Gaussian likelihood - the parameters are inferred from the
+ arguments of function
+
+ Parameters
+ ----------
+ x, y: array_like
+ The data to analyse
+ sigma: float
+ The standard deviation of the noise
+ function:
+ The python function to fit to the data. Note, this must take the
+ dependent variable as its first argument. The other arguments are
+ will require a prior and will be sampled over (unless a fixed
+ value is given).
+ """
+ self.x = x
+ self.y = y
+ self.sigma = sigma
+ self.N = len(x)
+ self.function = function
+
+ # These lines of code infer the parameters from the provided function
+ parameters = inspect.getargspec(function).args
+ parameters.pop(0)
+ self.parameters = dict.fromkeys(parameters)
+
+ def log_likelihood(self, sampler=None):
+ """
+ Parameters
+ ----------
+ sampler: :class:`tupak.core.sampler.Pymc3`
+ A Sampler object must be passed containing the prior distributions
+ and PyMC3 :class:`~pymc3.Model` to use as a context manager.
+ """
+
+ from tupak.core.sampler import Pymc3
+
+ if not isinstance(sampler, Pymc3):
+ raise ValueError("Sampler is not a tupak Pymc3 sampler object")
+
+ if not hasattr(sampler, 'pymc3_model'):
+ raise AttributeError("Sampler has not PyMC3 model attribute")
+
+ with sampler.pymc3_model:
+ mdist = sampler.pymc3_priors['m']
+ cdist = sampler.pymc3_priors['c']
+
+ mu = model(time, mdist, cdist)
+
+ # set the likelihood distribution
+ pm.Normal('likelihood', mu=mu, sd=self.sigma, observed=self.y)
+
+# Now lets instantiate a version of our GaussianLikelihood, giving it
+# the time, data and signal model
+likelihood = GaussianLikelihoodPyMC3(time, data, sigma, model)
+
+
+# Define a custom prior for one of the parameter for use with PyMC3
+class PriorPyMC3(tupak.core.prior.Prior):
+ def __init__(self, minimum, maximum, name=None, latex_label=None):
+ """
+ Uniform prior with bounds (should be equivalent to tupak.prior.Uniform)
+ """
+
+ tupak.core.prior.Prior.__init__(self, name, latex_label,
+ minimum=minimum,
+ maximum=maximum)
+
+ def ln_prob(self, sampler=None):
+ """
+ Change ln_prob method to take in a Sampler and return a PyMC3
+ distribution.
+ """
+
+ from tupak.core.sampler import Pymc3
+
+ if not isinstance(sampler, Pymc3):
+ raise ValueError("Sampler is not a tupak Pymc3 sampler object")
+
+ return pm.Uniform(self.name, lower=self.minimum,
+ upper=self.maximum)
+
+# From hereon, the syntax is exactly equivalent to other tupak examples
+# We make a prior
+priors = {}
+priors['m'] = tupak.core.prior.Uniform(0, 5, 'm')
+priors['c'] = PriorPyMC3(-2, 2, 'c')
+
+# And run sampler
+result = tupak.run_sampler(
+ likelihood=likelihood, priors=priors, sampler='pymc3', draws=1000,
+ tune=1000, discard_tuned_samples=True, injection_parameters=injection_parameters,
+ outdir=outdir, label=label)
+result.plot_corner()
diff --git a/examples/other_examples/radioactive_decay.py b/examples/other_examples/radioactive_decay.py
index 156406da22ca603927d662dacc3b1317c0f12bbd..bcc8fbc7c9ed77b95e9c5e24eda0026a5d6c4f76 100644
--- a/examples/other_examples/radioactive_decay.py
+++ b/examples/other_examples/radioactive_decay.py
@@ -20,6 +20,8 @@ tupak.utils.check_directory_exists_and_if_not_mkdir(outdir)
# generate a set of counts per minute for Polonium 214 with a half-life of 20 mins
halflife = 20
N0 = 1.e-19 # initial number of radionucleotides in moles
+atto = 1e-18
+N0 /= atto
def decayrate(deltat, halflife, N0):
"""
@@ -33,7 +35,9 @@ def decayrate(deltat, halflife, N0):
ln2 = np.log(2.)
NA = 6.02214078e23 # Avagadro's number
- rates = (N0*NA*(np.exp(-ln2*(times[0:-1]/halflife))
+ N0a = N0*atto*NA # number of nucleotides
+
+ rates = (N0a*(np.exp(-ln2*(times[0:-1]/halflife))
- np.exp(-ln2*(times[1:]/halflife)))/deltat)
return rates
@@ -70,8 +74,8 @@ likelihood = PoissonLikelihood(deltat, counts, decayrate)
priors = {}
priors['halflife'] = tupak.core.prior.LogUniform(1e-5, 1e5, 'halflife',
latex_label='$t_{1/2}$ (min)')
-priors['N0'] = tupak.core.prior.LogUniform(1e-25, 1e-10, 'N0',
- latex_label='$N_0$ (mole)')
+priors['N0'] = tupak.core.prior.LogUniform(1e-25/atto, 1e-10/atto, 'N0',
+ latex_label='$N_0$ (attomole)')
# And run sampler
result = tupak.run_sampler(
diff --git a/optional_requirements.txt b/optional_requirements.txt
index 9519ca78e7c2e4c511949f0556308ca6f36f20d5..bb4d2aadd2451de8043928164baa0c8e428fdf45 100644
--- a/optional_requirements.txt
+++ b/optional_requirements.txt
@@ -1,3 +1,4 @@
astropy
lalsuite
gwpy
+theano
\ No newline at end of file
diff --git a/tupak/core/prior.py b/tupak/core/prior.py
index eee8f98a8b763968164b6cee8ade061ef93d62cb..151883c6b66f850a818b60709e09fe9fe361df64 100644
--- a/tupak/core/prior.py
+++ b/tupak/core/prior.py
@@ -1056,7 +1056,6 @@ class Exponential(Prior):
"""
Prior.__init__(self, name=name, minimum=0., latex_label=latex_label)
self.mu = mu
- self.logmu = np.log(mu)
def rescale(self, val):
"""
diff --git a/tupak/core/sampler.py b/tupak/core/sampler.py
index aadbe4a0ebaf70f8ffc49ed033ed6238a2b7cfd0..fe558a4de1871b9ed93f8c026a9dd38fddeed406 100644
--- a/tupak/core/sampler.py
+++ b/tupak/core/sampler.py
@@ -941,6 +941,531 @@ class Ptemcee(Emcee):
return self.result
+class Pymc3(Sampler):
+ """ https://docs.pymc.io/ """
+
+ def _verify_parameters(self):
+ """
+ Change `_verify_parameters()` to just pass, i.e., don't try and
+ evaluate the likelihood for PyMC3.
+ """
+ pass
+
+ def _verify_use_ratio(self):
+ """
+ Change `_verify_use_ratio() to just pass.
+ """
+ pass
+
+ def _initialise_parameters(self):
+ """
+ Change `_initialise_parameters()`, so that it does call the `sample`
+ method in the Prior class.
+
+ """
+
+ self.__search_parameter_keys = []
+ self.__fixed_parameter_keys = []
+
+ for key in self.priors:
+ if isinstance(self.priors[key], Prior) \
+ and self.priors[key].is_fixed is False:
+ self.__search_parameter_keys.append(key)
+ elif isinstance(self.priors[key], Prior) \
+ and self.priors[key].is_fixed is True:
+ self.__fixed_parameter_keys.append(key)
+
+ logger.info("Search parameters:")
+ for key in self.__search_parameter_keys:
+ logger.info(' {} = {}'.format(key, self.priors[key]))
+ for key in self.__fixed_parameter_keys:
+ logger.info(' {} = {}'.format(key, self.priors[key].peak))
+
+ def _initialise_result(self):
+ """
+ Initialise results within Pymc3 subclass.
+ """
+ result = Result()
+ result.sampler = self.__class__.__name__.lower()
+ result.search_parameter_keys = self.__search_parameter_keys
+ result.fixed_parameter_keys = self.__fixed_parameter_keys
+ result.parameter_labels = [
+ self.priors[k].latex_label for k in
+ self.__search_parameter_keys]
+ result.label = self.label
+ result.outdir = self.outdir
+ result.kwargs = self.kwargs
+ return result
+
+ @property
+ def kwargs(self):
+ """ Ensures that proper keyword arguments are used for the Pymc3 sampler.
+
+ Returns
+ -------
+ dict: Keyword arguments used for the Nestle Sampler
+
+ """
+ return self.__kwargs
+
+ @kwargs.setter
+ def kwargs(self, kwargs):
+ self.__kwargs = dict()
+ self.__kwargs.update(kwargs)
+
+ # set some defaults
+
+ # set the number of draws
+ self.draws = 1000 if 'draws' not in self.__kwargs else self.__kwargs.pop('draws')
+
+ if 'chains' not in self.__kwargs:
+ self.__kwargs['chains'] = 2
+ self.chains = self.__kwargs['chains']
+
+ if 'cores' not in self.__kwargs:
+ self.__kwargs['cores'] = 1
+
+ def setup_prior_mapping(self):
+ """
+ Set the mapping between predefined tupak priors and the equivalent
+ PyMC3 distributions.
+ """
+
+ prior_map = {}
+ self.prior_map = prior_map
+
+ # predefined PyMC3 distributions
+ prior_map['Gaussian'] = {'pymc3': 'Normal',
+ 'argmap': {'mu': 'mu', 'sigma': 'sd'}}
+ prior_map['TruncatedGaussian'] = {'pymc3': 'TruncatedNormal',
+ 'argmap': {'mu': 'mu', 'sigma': 'sd', 'minimum': 'lower', 'maximum': 'upper'}}
+ prior_map['HalfGaussian'] = {'pymc3': 'HalfNormal',
+ 'argmap': {'sigma': 'sd'}}
+ prior_map['Uniform'] = {'pymc3': 'Uniform',
+ 'argmap': {'minimum': 'lower', 'maximum': 'upper'}}
+ prior_map['LogNormal'] = {'pymc3': 'Lognormal',
+ 'argmap': {'mu': 'mu', 'sigma': 'sd'}}
+ prior_map['Exponential'] = {'pymc3': 'Exponential',
+ 'argmap': {'mu': 'lam'},
+ 'argtransform': {'mu': lambda mu: 1./mu}}
+ prior_map['StudentT'] = {'pymc3': 'StudentT',
+ 'argmap': {'df': 'nu', 'mu': 'mu', 'scale': 'sd'}}
+ prior_map['Beta'] = {'pymc3': 'Beta',
+ 'argmap': {'alpha': 'alpha', 'beta': 'beta'}}
+ prior_map['Logistic'] = {'pymc3': 'Logistic',
+ 'argmap': {'mu': 'mu', 'scale': 's'}}
+ prior_map['Cauchy'] = {'pymc3': 'Cauchy',
+ 'argmap': {'alpha': 'alpha', 'beta': 'beta'}}
+ prior_map['Gamma'] = {'pymc3': 'Gamma',
+ 'argmap': {'k': 'alpha', 'theta': 'beta'},
+ 'argtransform': {'theta': lambda theta: 1./theta}}
+ prior_map['ChiSquared'] = {'pymc3': 'ChiSquared',
+ 'argmap': {'nu': 'nu'}}
+ prior_map['Interped'] = {'pymc3': 'Interpolated',
+ 'argmap': {'xx': 'x_points', 'yy': 'pdf_points'}}
+ prior_map['Normal'] = prior_map['Gaussian']
+ prior_map['TruncatedNormal'] = prior_map['TruncatedGaussian']
+ prior_map['HalfNormal'] = prior_map['HalfGaussian']
+ prior_map['LogGaussian'] = prior_map['LogNormal']
+ prior_map['Lorentzian'] = prior_map['Cauchy']
+ prior_map['FromFile'] = prior_map['Interped']
+
+ # internally defined mappings for tupak priors
+ prior_map['DeltaFunction'] = {'internal': self._deltafunction_prior}
+ prior_map['Sine'] = {'internal': self._sine_prior}
+ prior_map['Cosine'] = {'internal': self._cosine_prior}
+ prior_map['PowerLaw'] = {'internal': self._powerlaw_prior}
+ prior_map['LogUniform'] = {'internal': self._powerlaw_prior}
+
+ def _deltafunction_prior(self, key, **kwargs):
+ """
+ Map the tupak delta function prior to a single value for PyMC3.
+ """
+
+ from tupak.core.prior import DeltaFunction
+
+ # check prior is a DeltaFunction
+ if isinstance(self.priors[key], DeltaFunction):
+ return self.priors[key].peak
+ else:
+ raise ValueError("Prior for '{}' is not a DeltaFunction".format(key))
+
+ def _sine_prior(self, key):
+ """
+ Map the tupak Sine prior to a PyMC3 style function
+ """
+
+ from tupak.core.prior import Sine
+
+ # check prior is a Sine
+ if isinstance(self.priors[key], Sine):
+ pymc3 = self.external_sampler
+
+ try:
+ import theano.tensor as tt
+ from pymc3.theanof import floatX
+ except ImportError:
+ raise ImportError("You must have Theano installed to use PyMC3")
+
+ class Pymc3Sine(pymc3.Continuous):
+ def __init__(self, lower=0., upper=np.pi):
+ if lower >= upper:
+ raise ValueError("Lower bound is above upper bound!")
+
+ # set the mode
+ self.lower = lower = tt.as_tensor_variable(floatX(lower))
+ self.upper = upper = tt.as_tensor_variable(floatX(upper))
+ self.norm = (tt.cos(lower) - tt.cos(upper))
+ self.mean = (tt.sin(upper)+lower*tt.cos(lower) - tt.sin(lower) - upper*tt.cos(upper))/self.norm
+
+ transform = pymc3.distributions.transforms.interval(lower, upper)
+
+ super(Pymc3Sine, self).__init__(transform=transform)
+
+ def logp(self, value):
+ upper = self.upper
+ lower = self.lower
+ return pymc3.distributions.dist_math.bound(tt.log(tt.sin(value)/self.norm), lower <= value, value <= upper)
+
+ return Pymc3Sine(key, lower=self.priors[key].minimum, upper=self.priors[key].maximum)
+ else:
+ raise ValueError("Prior for '{}' is not a Sine".format(key))
+
+ def _cosine_prior(self, key):
+ """
+ Map the tupak Cosine prior to a PyMC3 style function
+ """
+
+ from tupak.core.prior import Cosine
+
+ # check prior is a Cosine
+ if isinstance(self.priors[key], Cosine):
+ pymc3 = self.external_sampler
+
+ # import theano
+ try:
+ import theano.tensor as tt
+ from pymc3.theanof import floatX
+ except ImportError:
+ raise ImportError("You must have Theano installed to use PyMC3")
+
+ class Pymc3Cosine(pymc3.Continuous):
+ def __init__(self, lower=-np.pi/2., upper=np.pi/2.):
+ if lower >= upper:
+ raise ValueError("Lower bound is above upper bound!")
+
+ self.lower = lower = tt.as_tensor_variable(floatX(lower))
+ self.upper = upper = tt.as_tensor_variable(floatX(upper))
+ self.norm = (tt.sin(upper) - tt.sin(lower))
+ self.mean = (upper*tt.sin(upper) + tt.cos(upper)-lower*tt.sin(lower)-tt.cos(lower))/self.norm
+
+ transform = pymc3.distributions.transforms.interval(lower, upper)
+
+ super(Pymc3Cosine, self).__init__(transform=transform)
+
+ def logp(self, value):
+ upper = self.upper
+ lower = self.lower
+ return pymc3.distributions.dist_math.bound(tt.log(tt.cos(value)/self.norm), lower <= value, value <= upper)
+
+ return Pymc3Cosine(key, lower=self.priors[key].minimum, upper=self.priors[key].maximum)
+ else:
+ raise ValueError("Prior for '{}' is not a Cosine".format(key))
+
+ def _powerlaw_prior(self, key):
+ """
+ Map the tupak PowerLaw prior to a PyMC3 style function
+ """
+
+ from tupak.core.prior import PowerLaw
+
+ # check prior is a PowerLaw
+ if isinstance(self.priors[key], PowerLaw):
+ pymc3 = self.external_sampler
+
+ # check power law is set
+ if not hasattr(self.priors[key], 'alpha'):
+ raise AttributeError("No 'alpha' attribute set for PowerLaw prior")
+
+ # import theano
+ try:
+ import theano.tensor as tt
+ from pymc3.theanof import floatX
+ except ImportError:
+ raise ImportError("You must have Theano installed to use PyMC3")
+
+ if self.priors[key].alpha < -1.:
+ # use Pareto distribution
+ palpha = -(1. + self.priors[key].alpha)
+
+ return pymc3.Bound(pymc3.Pareto, upper=self.priors[key].minimum)(key, alpha=palpha, m=self.priors[key].maximum)
+ else:
+ class Pymc3PowerLaw(pymc3.Continuous):
+ def __init__(self, lower, upper, alpha, testval=1):
+ falpha = alpha
+ self.lower = lower = tt.as_tensor_variable(floatX(lower))
+ self.upper = upper = tt.as_tensor_variable(floatX(upper))
+ self.alpha = alpha = tt.as_tensor_variable(floatX(alpha))
+
+ if falpha == -1:
+ self.norm = 1./(tt.log(self.upper/self.lower))
+ else:
+ beta = (1. + self.alpha)
+ self.norm = 1. /(beta * (tt.pow(self.upper, beta)
+ - tt.pow(self.lower, beta)))
+
+ transform = pymc3.distributions.transforms.interval(lower, upper)
+
+ super(Pymc3PowerLaw, self).__init__(transform=transform, testval=testval)
+
+ def logp(self, value):
+ upper = self.upper
+ lower = self.lower
+ alpha = self.alpha
+
+ return pymc3.distributions.dist_math.bound(self.alpha*tt.log(value) + tt.log(self.norm), lower <= value, value <= upper)
+
+ return Pymc3PowerLaw(key, lower=self.priors[key].minimum, upper=self.priors[key].maximum, alpha=self.priors[key].alpha)
+ else:
+ raise ValueError("Prior for '{}' is not a Power Law".format(key))
+
+ def _run_external_sampler(self):
+ pymc3 = self.external_sampler
+
+ # set the step method
+ from pymc3.sampling import STEP_METHODS
+
+ step_methods = {m.__name__.lower(): m.__name__ for m in STEP_METHODS}
+ if 'step' in self.__kwargs:
+ step_method = self.__kwargs.pop('step').lower()
+
+ if step_method not in step_methods:
+ raise ValueError("Using invalid step method '{}'".format(step_method))
+ else:
+ step_method = None
+
+ # initialise the PyMC3 model
+ self.pymc3_model = pymc3.Model()
+
+ # set the step method
+ sm = None if step_method is None else pymc3.__dict__[step_methods[step_method]]()
+
+ # set the prior
+ self.set_prior()
+
+ # if a custom log_likelihood function requires a `sampler` argument
+ # then use that log_likelihood function, with the assumption that it
+ # takes in a Pymc3 Sampler, with a pymc3_model attribute, and defines
+ # the likelihood within that context manager
+ likeargs = inspect.getargspec(self.likelihood.log_likelihood).args
+ if 'sampler' in likeargs:
+ self.likelihood.log_likelihood(sampler=self)
+ else:
+ # set the likelihood function from predefined functions
+ self.set_likelihood()
+
+ with self.pymc3_model:
+ # perform the sampling
+ trace = pymc3.sample(self.draws, step=sm, **self.kwargs)
+
+ nparams = len([key for key in self.priors.keys() if self.priors[key].__class__.__name__ != 'DeltaFunction'])
+ nsamples = len(trace)*self.chains
+
+ self.result.samples = np.zeros((nsamples, nparams))
+ count = 0
+ for key in self.priors.keys():
+ if self.priors[key].__class__.__name__ != 'DeltaFunction': # ignore DeltaFunction variables
+ self.result.samples[:,count] = trace[key]
+ count += 1
+
+ self.result.sampler_output = np.nan
+ self.calculate_autocorrelation(self.result.samples)
+ self.result.log_evidence = np.nan
+ self.result.log_evidence_err = np.nan
+ return self.result
+
+ def set_prior(self):
+ """
+ Set the PyMC3 prior distributions.
+ """
+
+ self.setup_prior_mapping()
+
+ self.pymc3_priors = dict()
+
+ pymc3 = self.external_sampler
+
+ # set the parameter prior distributions (in the model context manager)
+ with self.pymc3_model:
+ for key in self.priors:
+ # if the prior contains ln_prob method that takes a 'sampler' argument
+ # then try using that
+ lnprobargs = inspect.getargspec(self.priors[key].ln_prob).args
+ if 'sampler' in lnprobargs:
+ try:
+ self.pymc3_priors[key] = self.priors[key].ln_prob(sampler=self)
+ except RuntimeError:
+ raise RuntimeError(("Problem setting PyMC3 prior for ",
+ "'{}'".format(key)))
+ else:
+ # use Prior distribution name
+ distname = self.priors[key].__class__.__name__
+
+ if distname in self.prior_map:
+ # check if we have a predefined PyMC3 distribution
+ if 'pymc3' in self.prior_map[distname] and 'argmap' in self.prior_map[distname]:
+ # check the required arguments for the PyMC3 distribution
+ pymc3distname = self.prior_map[distname]['pymc3']
+
+ if pymc3distname not in pymc3.__dict__:
+ raise ValueError("Prior '{}' is not a known PyMC3 distribution.".format(pymc3distname))
+
+ reqargs = inspect.getargspec(pymc3.__dict__[pymc3distname].__init__).args[1:]
+
+ # set keyword arguments
+ priorkwargs = {}
+ for (targ, parg) in self.prior_map[distname]['argmap'].items():
+ if hasattr(self.priors[key], targ):
+ if parg in reqargs:
+ if 'argtransform' in self.prior_map[distname]:
+ if targ in self.prior_map[distname]['argtransform']:
+ tfunc = self.prior_map[distname]['argtransform'][targ]
+ else:
+ tfunc = lambda x: x
+ else:
+ tfunc = lambda x: x
+
+ priorkwargs[parg] = tfunc(getattr(self.priors[key], targ))
+ else:
+ raise ValueError("Unknown argument {}".format(parg))
+ else:
+ if parg in reqargs:
+ priorkwargs[parg] = None
+ self.pymc3_priors[key] = pymc3.__dict__[pymc3distname](key, **priorkwargs)
+ elif 'internal' in self.prior_map[distname]:
+ self.pymc3_priors[key] = self.prior_map[distname]['internal'](key)
+ else:
+ raise ValueError("Prior '{}' is not a known distribution.".format(distname))
+ else:
+ raise ValueError("Prior '{}' is not a known distribution.".format(distname))
+
+ def set_likelihood(self):
+ """
+ Convert any tupak likelihoods to PyMC3 distributions.
+ """
+
+ pymc3 = self.external_sampler
+
+ with self.pymc3_model:
+ # check if it is a predefined likelhood function
+ if self.likelihood.__class__.__name__ == 'GaussianLikelihood':
+ # check required attributes exist
+ if (not hasattr(self.likelihood, 'sigma') or
+ not hasattr(self.likelihood, 'x') or
+ not hasattr(self.likelihood, 'y') or
+ not hasattr(self.likelihood, 'function') or
+ not hasattr(self.likelihood, 'function_keys')):
+ raise ValueError("Gaussian Likelihood does not have all the correct attributes!")
+
+ if 'sigma' in self.pymc3_priors:
+ # if sigma is suppled use that value
+ if self.likelihood.sigma is None:
+ self.likelihood.sigma = self.pymc3_priors.pop('sigma')
+ else:
+ del self.pymc3_priors['sigma']
+
+ for key in self.pymc3_priors:
+ if key not in self.likelihood.function_keys:
+ raise ValueError("Prior key '{}' is not a function key!".format(key))
+
+ model = self.likelihood.function(self.likelihood.x, **self.pymc3_priors)
+
+ # set the distribution
+ pymc3.Normal('likelihood', mu=model, sd=self.likelihood.sigma,
+ observed=self.likelihood.y)
+ elif self.likelihood.__class__.__name__ == 'PoissonLikelihood':
+ # check required attributes exist
+ if (not hasattr(self.likelihood, 'x') or
+ not hasattr(self.likelihood, 'y') or
+ not hasattr(self.likelihood, 'function') or
+ not hasattr(self.likelihood, 'function_keys')):
+ raise ValueError("Poisson Likelihood does not have all the correct attributes!")
+
+ for key in self.pymc3_priors:
+ if key not in self.likelihood.function_keys:
+ raise ValueError("Prior key '{}' is not a function key!".format(key))
+
+ # get rate function
+ model = self.likelihood.function(self.likelihood.x, **self.pymc3_priors)
+
+ # set the distribution
+ pymc3.Poisson('likelihood', mu=model, observed=self.likelihood.y)
+ elif self.likelihood.__class__.__name__ == 'ExponentialLikelihood':
+ # check required attributes exist
+ if (not hasattr(self.likelihood, 'x') or
+ not hasattr(self.likelihood, 'y') or
+ not hasattr(self.likelihood, 'function') or
+ not hasattr(self.likelihood, 'function_keys')):
+ raise ValueError("Exponential Likelihood does not have all the correct attributes!")
+
+ for key in self.pymc3_priors:
+ if key not in self.likelihood.function_keys:
+ raise ValueError("Prior key '{}' is not a function key!".format(key))
+
+ # get mean function
+ model = self.likelihood.function(self.likelihood.x, **self.pymc3_priors)
+
+ # set the distribution
+ pymc3.Exponential('likelihood', lam=1./model, observed=self.likelihood.y)
+ elif self.likelihood.__class__.__name__ == 'StudentTLikelihood':
+ # check required attributes exist
+ if (not hasattr(self.likelihood, 'x') or
+ not hasattr(self.likelihood, 'y') or
+ not hasattr(self.likelihood, 'nu') or
+ not hasattr(self.likelihood, 'sigma') or
+ not hasattr(self.likelihood, 'function') or
+ not hasattr(self.likelihood, 'function_keys')):
+ raise ValueError("StudentT Likelihood does not have all the correct attributes!")
+
+ if 'nu' in self.pymc3_priors:
+ # if nu is suppled use that value
+ if self.likelihood.nu is None:
+ self.likelihood.nu = self.pymc3_priors.pop('nu')
+ else:
+ del self.pymc3_priors['nu']
+
+ for key in self.pymc3_priors:
+ if key not in self.likelihood.function_keys:
+ raise ValueError("Prior key '{}' is not a function key!".format(key))
+
+ model = self.likelihood.function(self.likelihood.x, **self.pymc3_priors)
+
+ # set the distribution
+ pymc3.StudentT('likelihood', nu=self.likelihood.nu, mu=model, sd=self.likelihood.sigma, observed=self.likelihood.y)
+ else:
+ raise ValueError("Unknown likelihood has been provided")
+
+ def calculate_autocorrelation(self, samples, c=3):
+ """ Uses the `emcee.autocorr` module to estimate the autocorrelation
+
+ Parameters
+ ----------
+ c: float
+ The minimum number of autocorrelation times needed to trust the
+ estimate (default: `3`). See `emcee.autocorr.integrated_time`.
+ """
+
+ import emcee
+ try:
+ self.result.max_autocorrelation_time = int(np.max(
+ emcee.autocorr.integrated_time(samples, c=c)))
+ logger.info("Max autocorr time = {}".format(
+ self.result.max_autocorrelation_time))
+ except emcee.autocorr.AutocorrError as e:
+ self.result.max_autocorrelation_time = None
+ logger.info("Unable to calculate autocorr time: {}".format(e))
+
+
def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
sampler='dynesty', use_ratio=None, injection_parameters=None,
conversion_function=None, plot=False, default_priors_file=None,