diff --git a/tupak/core/likelihood.py b/tupak/core/likelihood.py
index c6048b40d38622a9f38ce62d5177201ca297cc1c..5784c859b53feee1ae811cea259aa0e634572187 100644
--- a/tupak/core/likelihood.py
+++ b/tupak/core/likelihood.py
@@ -137,7 +137,10 @@ class GaussianLikelihood(Likelihood):
model = self.function(self.x, **sampler.pymc3_priors)
- return sampler.external_sampler.Normal('likelihood', mu=model, sd=self.sigma, observed=self.y)
+ # set the distribution
+ sampler.external_sampler.Normal('likelihood', mu=model,
+ sd=self.sigma, observed=self.y)
+
class PoissonLikelihood(Likelihood):
def __init__(self, x, counts, func):
@@ -258,5 +261,6 @@ class PoissonLikelihood(Likelihood):
# get rate function
rate = self.function(self.x, **sampler.pymc3_priors)
- return sampler.external_sampler.Poisson('likelihood', mu=rate,
- observed=self.y)
+ # set the distribution
+ sampler.external_sampler.Poisson('likelihood', mu=rate,
+ observed=self.y)
diff --git a/tupak/core/prior.py b/tupak/core/prior.py
index d48617970b81683ca4f6fcc8001e47d91be63706..0715476f1be98c69868be17f26595d91467a39ba 100644
--- a/tupak/core/prior.py
+++ b/tupak/core/prior.py
@@ -541,10 +541,6 @@ class DeltaFunction(Prior):
else:
return 0
- def pymc3_prior(self, sampler):
- # just return the value
- return self.peak
-
def __repr__(self):
"""Call to helper method in the super class."""
return Prior._subclass_repr_helper(self, subclass_args=['peak'])
@@ -651,10 +647,6 @@ class Uniform(Prior):
"""
Prior.__init__(self, name, latex_label, minimum, maximum)
- # set PyMC3 Uniform distribution attributes
- self.lower = self.minimum
- self.upper = self.maximum
-
def rescale(self, val):
Prior.test_valid_for_rescaling(val)
return self.minimum + val * (self.maximum - self.minimum)
@@ -860,17 +852,11 @@ class Gaussian(Prior):
class Normal(Gaussian):
def __init__(self, mu, sigma, name=None, latex_label=None):
- """A copy of the Gaussian prior, but with "Normal" name to copy that
- used for the distribution in PyMC3.
-
+ """A synonym for the Gaussian prior
"""
Gaussian.__init__(self, mu, sigma, name, latex_label)
- # set argument names used in PyMC3 distribution
- self.mu = mu
- self.sd = sigma
-
class TruncatedGaussian(Prior):
@@ -939,6 +925,14 @@ class TruncatedGaussian(Prior):
return Prior._subclass_repr_helper(self, subclass_args=['mu', 'sigma'])
+class TruncatedNormal(Gaussian):
+ def __init__(self, mu, sigma, minimum, maximum, name=None, latex_label=None):
+ """A synonym for the Truncated Gaussian prior
+ """
+
+ TruncatedGaussian.__init__(self, mu, sigma, minimum, maximum, name, latex_label)
+
+
class Interped(Prior):
def __init__(self, xx, yy, minimum=np.nan, maximum=np.nan, name=None, latex_label=None):
diff --git a/tupak/core/sampler.py b/tupak/core/sampler.py
index b3cee6884c4acfd5f010de928d323883c3b8263d..35eea0c933b45673987acd57b091235603ddc404 100644
--- a/tupak/core/sampler.py
+++ b/tupak/core/sampler.py
@@ -11,7 +11,7 @@ from collections import OrderedDict
from tupak.core.utils import logger
from tupak.core.result import Result, read_in_result
-from tupak.core.prior import Prior, DeltaFunction
+from tupak.core.prior import Prior
from tupak.core import utils
import tupak
@@ -944,6 +944,144 @@ class Ptemcee(Emcee):
class Pymc3(Sampler):
""" https://docs.pymc.io/ """
+ def setup_prior_mapping(self):
+ """
+ Set the mapping between predefined tupak priors and the equivalent
+ PyMC3 distributions.
+ """
+
+ self.prior_map = {}
+
+ # predefined PyMC3 distributions
+
+ # Gaussian
+ self.prior_map['Gaussian'] = {'pymc3': 'Normal',
+ 'argmap': {'mu': 'mu', 'sigma': 'sd'}}
+ self.prior_map['Normal'] = self.prior_map['Gaussian']
+
+ # Truncated Gaussian
+ self.prior_map['TruncatedGaussian'] = {'pymc3': 'TruncatedNormal',
+ 'argmap': {'mu': 'mu', 'sigma': 'sd', 'minimum': 'lower', 'maximum': 'upper'}}
+ self.prior_map['TruncatedNormal'] = self.prior_map['TruncatedGaussian']
+
+ # Uniform
+ self.prior_map['Uniform'] = {'pymc3': 'Uniform',
+ 'argmap': {'minimum': 'lower', 'maximum': 'upper'}}
+
+ # Interpolated
+ self.prior_map['Interped'] = {'pymc3': 'Interpolated',
+ 'argmap': {'xx': 'x_points', 'yy': 'pdf_points'}}
+
+ # internally defined mappings for tupak priors
+
+ # DeltaFunction
+ self.prior_map['DeltaFunction'] = {'internal': self._deltafunction_prior}
+
+ # Sine
+ self.prior_map['Sine'] = {'internal': self._sine_prior}
+
+ # Cosine
+ self.prior_map['Cosine'] = {'internal': self._cosine_prior}
+
+ def _deltafunction_prior(self, key):
+ """
+ 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
+
+ # 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 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.norm = (np.sin(upper) - np.sin(lower))
+ self.mean = (upper*np.sin(upper)+np.cos(upper)-lower*np.sin(lower)-np.cos(lower))/self.norm
+
+ 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 Sine".format(key))
+
def _run_external_sampler(self):
pymc3 = self.external_sampler
@@ -954,20 +1092,21 @@ class Pymc3(Sampler):
self.tune = self.kwargs.get('tune', 1000) # burn in samples
self.discard_tuned_samples = self.kwargs.get('discard_tuned_samples',
True)
+
# set the model
model = pymc3.Model()
with model:
self.set_prior()
- likelihood = self.likelihood.pymc3_likelihood(self)
+ self.likelihood.pymc3_likelihood(self)
# perform the sampling
trace = pymc3.sample(self.draws, tune=self.tune, cores=self.cores,
chains=self.chains,
discard_tuned_samples=self.discard_tuned_samples)
- nparams = len(trace.varnames)
+ nparams = len(self.priors.keys())
nsamples = len(trace)*self.chains
self.result.samples = np.zeros((nsamples, nparams))
@@ -985,8 +1124,12 @@ class Pymc3(Sampler):
"""
+ self.setup_prior_mapping()
+
self.pymc3_priors = dict()
+ pymc3 = self.external_sampler
+
# set the parameter prior distributions
for key in self.priors:
# if the prior contains a pymc3_prior method use that otherwise try
@@ -997,24 +1140,32 @@ class Pymc3(Sampler):
# use Prior distribution name
distname = self.priors[key].__class__.__name__
- # check whether name is a PyMC3 distribution
- if distname in self.external_sampler.__dict__:
- # check the required arguments for the PyMC3 distribution
- reqargs = inspect.getargspec(self.external_sampler.__dict__[distname].__init__).args[1:]
-
- priorkwargs = dict()
-
- # check whether the Prior class has required attributes
- for arg in reqargs:
- if hasattr(self.priors[key], arg):
- priorkwargs[arg] = getattr(self.priors[key], arg)
- else:
- priorkwargs[arg] = None
-
- # set the prior
- self.pymc3_priors[key] = self.external_sampler.__dict__[distname](key, **priorkwargs)
+ 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] and distname in pymc3.__dict__:
+ # check the required arguments for the PyMC3 distribution
+ pymc3distname = self.prior_map[distname]['pymc3']
+
+ 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:
+ priorkwargs[parg] = 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 PyMC3 distribution.".format(distname))
+ raise ValueError("Prior '{}' is not a known distribution.".format(distname))
def calculate_autocorrelation(self, samples, c=3):
""" Uses the `emcee.autocorr` module to estimate the autocorrelation