This example notebook illustrates some of the functionality that is available for LogPrior
objects that are currently available within PINTS.
from __future__ import print_function
import pints
import numpy as np
import matplotlib.pyplot as plt
The uniform prior, $U\sim(a, b)$, here with $a=-10$ and $b=15$. When this function is called, its log density is returned.
uniform_log_prior = pints.UniformLogPrior(-10, 15)
print('U(0|a=-10, b=15) = ' + str(uniform_log_prior([0])))
U(0|a=-10, b=15) = -3.2188758248682006
To plot the density, we take the exponential of its log density.
values = np.linspace(-20, 20, 1000)
log_prob = [uniform_log_prior([x]) for x in values]
prob = np.exp(log_prob)
plt.figure(figsize=(10,4))
plt.xlabel('theta')
plt.ylabel('Density')
plt.plot(values, prob)
plt.show()
To specify a multidimensional uniform prior, use the same function. Here we specify, $\theta_1\sim U(2, 4)$ and $\theta_2\sim U(-7,-5)$.
uniform_log_prior = pints.UniformLogPrior([2, -7], [4, -5])
Plot $p(\theta_2|\theta_1 = 3)$.
values = np.linspace(-10, -4, 1000)
log_prob = [uniform_log_prior([3, x]) for x in values]
prob = np.exp(log_prob)
plt.figure(figsize=(10,4))
plt.xlabel('theta[2]')
plt.ylabel('Density')
plt.plot(values, prob)
plt.show()
If you have a prior constrained to lie $\in[0,1]$, you can use a beta prior.
beta_log_prior1 = pints.BetaLogPrior(1, 1)
beta_log_prior2 = pints.BetaLogPrior(5, 3)
beta_log_prior3 = pints.BetaLogPrior(3, 5)
beta_log_prior4 = pints.BetaLogPrior(10, 10)
values = np.linspace(0, 1, 1000)
prob1 = np.exp([beta_log_prior1([x]) for x in values])
prob2 = np.exp([beta_log_prior2([x]) for x in values])
prob3 = np.exp([beta_log_prior3([x]) for x in values])
prob4 = np.exp([beta_log_prior4([x]) for x in values])
plt.figure(figsize=(10,4))
plt.xlabel('theta')
plt.ylabel('Density')
plt.plot(values, prob1)
plt.plot(values, prob2)
plt.plot(values, prob3)
plt.plot(values, prob4)
plt.legend(['beta(1, 1)', 'beta(5, 3)', 'beta(3, 5)', 'beta(10, 10)'])
plt.show()
Specifying a value outside the support of the distribution returns $-\infty$ for the log density.
print('beta(-0.5|a=1, b=1) = ' + str(beta_log_prior1([-0.5])))
beta(-0.5|a=1, b=1) = -inf
Each prior has a mean
function that allows you to quickly check what parameterisation is being used.
print('mean = ' + str(beta_log_prior3.mean()))
mean = 0.375
Each prior also has a sample
function which allows generation of independent samples from each distribution. Using this we can sample from a Student-t density, with input dimensions (location, degrees of freedom, scale)
.
n = 10000
student_t_log_prior = pints.StudentTLogPrior(10, 8, 5)
samples = student_t_log_prior.sample(n)
plt.hist(samples, 20)
plt.xlabel('theta')
plt.ylabel('Frequency')
plt.show()
For models with multiple parameters, we can specify different distributions for each dimension using ComposedLogPrior
.
log_prior1 = pints.GaussianLogPrior(6, 3)
log_prior2 = pints.InverseGammaLogPrior(5, 5)
log_prior3 = pints.LogNormalLogPrior(-1, 1)
composed_log_prior = pints.ComposedLogPrior(log_prior1, log_prior2, log_prior3)
# calling
composed_log_prior([-3, 1, 6])
-13.25607355949521
Functions like sample
and mean
also work for ComposedLogPrior
objects.
print('mean = ' + str(composed_log_prior.mean()))
n = 10
samples = composed_log_prior.sample(1000)
plt.hist(samples[:, 0], alpha=0.5)
plt.hist(samples[:, 1], alpha=0.5)
plt.hist(samples[:, 2], alpha=0.5)
plt.legend(['Gaussian(6, 3)', 'InverseGamma(5, 5)', 'LogNormal(-1, 1)'])
plt.xlabel('theta')
plt.ylabel('Frequency')
plt.show()
mean = [6.0, 1.25, 0.6065306597126334]
We also have multivariate priors in PINTS. For example, the multivariate Gaussian.
two_d_gaussian_log_prior = pints.MultivariateGaussianLogPrior([0, 10], [[1, 0.5], [0.5, 3]])
# Contour plot of pdf
x = np.linspace(-3, 3, 100)
y = np.linspace(4, 15, 100)
X, Y = np.meshgrid(x, y)
Z = np.exp([[two_d_gaussian_log_prior([i, j]) for i in x] for j in y])
plt.contour(X, Y, Z)
plt.xlabel('theta[2]')
plt.ylabel('theta[1]')
plt.show()