# Source code for orbitize.priors

import numpy as np
import sys
import abc

"""
This module defines priors with methods to draw samples and compute log(probability)
"""

# Python 2 & 3 handle ABCs differently
if sys.version_info[0] < 3:
ABC = abc.ABCMeta('ABC', (), {})
else:
ABC = abc.ABC

[docs]class Prior(ABC):
"""
Abstract base class for prior objects.
All prior objects should inherit from this class.

Written: Sarah Blunt, 2018
"""

@abc.abstractmethod
def draw_samples(self, num_samples):
pass

@abc.abstractmethod
def compute_lnprob(self, element_array):
pass

[docs]class GaussianPrior(Prior):
"""Gaussian prior.

.. math::

log(p(x|\\sigma, \\mu)) \\propto \\frac{(x - \\mu)}{\\sigma}

Args:
mu (float): mean of the distribution
sigma (float): standard deviation of the distribution
no_negatives (bool): if True, only positive values will be drawn from
this prior, and the probability of negative values will be 0 (default:True).

(written) Sarah Blunt, 2018
"""
def __init__(self, mu, sigma, no_negatives=True):
self.mu = mu
self.sigma = sigma
self.no_negatives = no_negatives

def __repr__(self):
return "Gaussian"

[docs]    def draw_samples(self, num_samples):
"""
Draw positive samples from a Gaussian distribution.
Negative samples will not be returned.

Args:
num_samples (float): the number of samples to generate

Returns:
numpy array of float: samples drawn from the appropriate
Gaussian distribution. Array has length num_samples.
"""

samples = np.random.normal(
loc=self.mu, scale=self.sigma, size=num_samples
)

if self.no_negatives:

)

return samples

[docs]    def compute_lnprob(self, element_array):
"""
Compute log(probability) of an array of numbers wrt a Gaussian distibution.
Negative numbers return a probability of -inf.

Args:
element_array (float or np.array of float): array of numbers. We want the
probability of drawing each of these from the appopriate Gaussian
distribution

Returns:
numpy array of float: array of log(probability) values,
corresponding to the probability of drawing each of the numbers
in the input element_array.
"""
lnprob = -0.5*np.log(2.*np.pi*self.sigma) - 0.5*((element_array - self.mu) / self.sigma)**2

if self.no_negatives:

return lnprob

[docs]class LogUniformPrior(Prior):
"""
This is the probability distribution :math:p(x) \\propto 1/x

The __init__ method should take in a "min" and "max" value
of the distribution, which correspond to the domain of the prior.
(If this is not implemented, the prior has a singularity at 0 and infinite
integrated probability).

Args:
minval (float): the lower bound of this distribution
maxval (float): the upper bound of this distribution

"""
def __init__(self, minval, maxval):
self.minval = minval
self.maxval = maxval

self.logmin = np.log(minval)
self.logmax = np.log(maxval)

def __repr__(self):
return "Log Uniform"

[docs]    def draw_samples(self, num_samples):
"""
Draw samples from this 1/x distribution.

Args:
num_samples (float): the number of samples to generate

Returns:
np.array:  samples ranging from [minval, maxval) as floats.
"""
# sample from a uniform distribution in log space
samples = np.random.uniform(self.logmin, self.logmax, num_samples)

# convert from log space to linear space
samples = np.exp(samples)

return samples

[docs]    def compute_lnprob(self, element_array):
"""
Compute the prior probability of each element given that its drawn from a Log-Uniofrm  prior

Args:
element_array (float or np.array of float): array of paramters to compute the prior probability of

Returns:
np.array: array of prior probabilities
"""
normalizer = self.logmax - self.logmin

lnprob = -np.log((element_array*normalizer))

# account for scalar inputs
if np.shape(lnprob) == ():
if (element_array > self.maxval) or (element_array < self.minval):
lnprob = -np.inf
else:
lnprob[(element_array > self.maxval) | (element_array < self.minval)] = -np.inf

return lnprob

[docs]class UniformPrior(Prior):
"""
This is the probability distribution p(x) propto constant.

Args:
minval (float): the lower bound of the uniform prior
maxval (float): the upper bound of the uniform prior

"""
def __init__(self, minval, maxval):
self.minval = minval
self.maxval = maxval

def __repr__(self):
return "Uniform"

[docs]    def draw_samples(self, num_samples):
"""
Draw samples from this uniform distribution.

Args:
num_samples (float): the number of samples to generate

Returns:
np.array:  samples ranging from [0, pi) as floats.
"""
# sample from a uniform distribution in log space
samples = np.random.uniform(self.minval, self.maxval, num_samples)

return samples

[docs]    def compute_lnprob(self, element_array):
"""
Compute the prior probability of each element given that its drawn from this uniform prior

Args:
element_array (float or np.array of float): array of paramters to compute the prior probability of

Returns:
np.array: array of prior probabilities
"""
lnprob = np.log(np.ones(np.size(element_array))/(self.maxval - self.minval))

# account for scalar inputs
if np.shape(lnprob) == ():
if (element_array > self.maxval) or (element_array < self.minval):
lnprob = -np.inf
else:
lnprob[(element_array > self.maxval) | (element_array < self.minval)] = -np.inf

return lnprob

[docs]class SinPrior(Prior):
"""
This is the probability distribution :math:p(x) \\propto sin(x)

The domain of this prior is [0,pi].
"""

def __init__(self):
pass

def __repr__(self):
return "Sine"

[docs]    def draw_samples(self, num_samples):
"""
Draw samples from a Sine distribution.

Args:
num_samples (float): the number of samples to generate

Returns:
np.array:  samples ranging from [0, pi) as floats.
"""

# draw uniform from -1 to 1
samples = np.random.uniform(-1, 1, num_samples)

samples = np.arccos(samples) % np.pi

return samples

[docs]    def compute_lnprob(self, element_array):
"""
Compute the prior probability of each element given that its drawn from a sine prior

Args:
element_array (float or np.array of float): array of paramters to compute the prior probability of

Returns:
np.array: array of prior probabilities
"""
normalization = 2.

lnprob = np.log(np.sin(element_array)/normalization)

# account for scalar inputs
if np.shape(lnprob) == ():
if (element_array>=np.pi) or (element_array<=0):
lnprob = -np.inf
else:
lnprob[(element_array>=np.pi) | (element_array<=0)] = -np.inf

return lnprob

[docs]class LinearPrior(Prior):
"""
Draw samples from the probability distribution:

.. math::

p(x) \\propto mx+b

where m is negative, b is positive, and the
range is [0,-b/m].

Args:
m (float): slope of line. Must be negative.
b (float): y intercept of line. Must be positive.

"""
def __init__(self, m, b):
self.m = m
self.b = b

def __repr__(self):
return "Linear"

[docs]    def draw_samples(self, num_samples):
"""
Draw samples from a descending linear distribution.

Args:
num_samples (float): the number of samples to generate

Returns:
samples (np.array):  samples ranging from [0, -b/m) as floats.
"""
norm = -0.5*self.b**2/self.m

# draw uniform from 0 to 1
samples = np.random.uniform(0, 1, num_samples)

# generate samples following a linear distribution
linear_samples = -np.sqrt(2.*norm*samples/self.m + (self.b/self.m)**2) - (self.b/self.m)

return linear_samples

def compute_lnprob(self, element_array):

x_intercept = -self.b/self.m
normalizer = -0.5*self.b**2/self.m

lnprob = np.log((self.m*element_array + self.b)/normalizer)

# account for scalar inputs
if np.shape(lnprob) == ():
if (element_array>=x_intercept) or (element_array<0):
lnprob = -np.inf
else:
lnprob[(element_array>=x_intercept) | (element_array<0)] = -np.inf

return lnprob

[docs]def all_lnpriors(params, priors):
"""
Calculates log(prior probability) of a set of parameters and a list of priors

Args:
params (np.array): size of N parameters
priors (list): list of N prior objects corresponding to each parameter

Returns:
float: prior probability of this set of parameters
"""
logp = 0.
for param, prior in zip(params, priors):
param = np.array([param])

logp += prior.compute_lnprob(param) # retrun a float

return logp

if __name__ == '__main__':

myPrior = LinearPrior(-1., 1.)
mySamples = myPrior.draw_samples(1000)
print(mySamples)
myProbs = myPrior.compute_lnprob(mySamples)
print(myProbs)

myPrior = GaussianPrior(1.3, 0.2)
mySamples = myPrior.draw_samples(1)
print(mySamples)

myProbs = myPrior.compute_lnprob(mySamples)
print(myProbs)