# Copyright (C) 2019 Steven Reyes
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
"""
This modules provides functions for estimating the marginal
likelihood or evidence of a model.
"""
import numpy
from scipy import integrate
# numpy renamed trapz to trapezoid in 2.0 and removed trapz in 2.x
try:
from numpy import trapezoid
except ImportError: # numpy < 2.0
from numpy import trapz as trapezoid
[docs]
def arithmetic_mean_estimator(log_likelihood):
"""Returns the log evidence via the prior arithmetic mean estimator (AME).
The logarithm form of AME is used. This is the most basic
evidence estimator, and often requires O(billions) of samples
from the prior.
Parameters
----------
log_likelihood : 1d array of floats
The log likelihood of the data sampled from the prior
distribution.
Returns
-------
float :
Estimation of the log of the evidence.
"""
num_samples = len(log_likelihood)
logl_max = numpy.max(log_likelihood)
log_evidence = 0.
for i, _ in enumerate(log_likelihood):
log_evidence += numpy.exp(log_likelihood[i] - logl_max)
log_evidence = numpy.log(log_evidence)
log_evidence += logl_max - numpy.log(num_samples)
return log_evidence
[docs]
def harmonic_mean_estimator(log_likelihood):
"""Returns the log evidence via posterior harmonic mean estimator (HME).
The logarithm form of HME is used. This method is not
recommended for general use. It is very slow to converge,
formally, has infinite variance, and very error prone.
Not recommended for general use.
Parameters
----------
log_likelihood : 1d array of floats
The log likelihood of the data sampled from the posterior
distribution.
Returns
-------
float :
Estimation of the log of the evidence.
"""
num_samples = len(log_likelihood)
logl_max = numpy.max(-1.0*log_likelihood)
log_evidence = 0.
for i, _ in enumerate(log_likelihood):
log_evidence += numpy.exp(-1.0*log_likelihood[i] + logl_max)
log_evidence = -1.0*numpy.log(log_evidence)
log_evidence += logl_max
log_evidence += numpy.log(num_samples)
return log_evidence
# numpy.trapz was renamed to numpy.trapezoid in numpy 2.0.
try:
from numpy import trapezoid as _trapezoid
except ImportError: # numpy < 2.0
from numpy import trapz as _trapezoid
[docs]
def mean_logl_by_temperature(logls, betas):
"""The mean log likelihood at each distinct inverse temperature.
Parameters
----------
logls : numpy.ndarray
Log likelihoods of shape (ntemps, nwalkers, niterations).
betas : numpy.ndarray
The inverse temperatures, of shape (ntemps, niterations); a ladder
that adapts visits more than one temperature per chain.
Returns
-------
betas : numpy.ndarray
Each distinct inverse temperature.
mean_logls : numpy.ndarray
The mean log likelihood at each of them.
"""
mean_logls = []
unique_betas = []
for ti in range(betas.shape[0]):
ubti, idx = numpy.unique(betas[ti, :], return_inverse=True)
unique_idx = numpy.unique(idx)
loglsti = logls[ti, :, :]
for ii in unique_idx:
# average over the walkers and iterations at this temperature
getiters = numpy.where(ii == unique_idx)[0]
mean_logls.append(loglsti[:, getiters].mean())
unique_betas.append(ubti[ii])
return numpy.array(unique_betas), numpy.array(mean_logls)
[docs]
def ladder_thermodynamic_integration(betas, logls):
"""Thermodynamic integration estimate of the evidence.
This is the same estimator used by the ``ptemcee`` sampler; see
:py:func:`pycbc.inference.sampler.ptemcee` for details.
Parameters
----------
betas : array
The inverse temperatures to use for the quadrature.
logls : array
The mean log-likelihoods corresponding to ``betas``.
Returns
-------
logZ : float
Estimate of the log-evidence.
dlogZ : float
The error associated with the finite number of temperatures at which
the posterior has been sampled.
"""
if len(betas) != len(logls):
raise ValueError("Need the same number of log(L) values as "
"temperatures.")
order = numpy.argsort(betas)[::-1]
betas = betas[order]
logls = logls[order]
betas0 = numpy.copy(betas)
if betas[-1] != 0:
betas = numpy.concatenate((betas0, [0]))
betas2 = numpy.concatenate((betas0[::2], [0]))
logls2 = numpy.concatenate((logls[::2], [logls[-1]]))
logls = numpy.concatenate((logls, [logls[-1]]))
else:
betas2 = numpy.concatenate((betas0[:-1:2], [0]))
logls2 = numpy.concatenate((logls[:-1:2], [logls[-1]]))
logZ = -_trapezoid(logls, betas)
logZ2 = -_trapezoid(logls2, betas2)
return logZ, numpy.abs(logZ - logZ2)
[docs]
def thermodynamic_integration(log_likelihood, betas,
method="simpsons"):
"""Returns the log evidence of the model via thermodynamic integration.
Also returns an estimated standard deviation for the log evidence.
Current options are integration through the trapezoid rule, a
first-order corrected trapezoid rule, and Simpson's rule.
Parameters
----------
log_likelihood : numpy.ndarray
3d array of shape (ntemps, nwalkers, niterations)
The log likelihood for each temperature separated by
temperature, walker, and iteration.
betas : numpy.ndarray
The inverse temperatures used in the MCMC.
method : string
Optional. Can be one of {"trapezoid", "trapezoid_corrected", "simpsons"}
The numerical integration method to use for the
thermodynamic integration. Choices include: "trapezoid",
"trapezoid_corrected", "simpsons", for the trapezoid rule,
the first-order correction to the trapezoid rule, and
Simpson's rule. [Default = "simpsons"]
Returns
-------
log_evidence : float
Estimation of the log of the evidence.
mcmc_std : float
The standard deviation of the log evidence estimate from
Monte-Carlo spread.
"""
# Check if the method of integration is in the list of choices
method_list = ["trapezoid", "trapezoid_corrected", "simpsons"]
if method not in method_list:
raise ValueError("Method %s not supported. Expected %s"
% (method, method_list))
# Read in the data and ensure ordering of data.
# Ascending order sort
order = numpy.argsort(betas)
betas = betas[order]
log_likelihood = log_likelihood[order]
# Assume log likelihood is given in shape of beta, walker,
# and iteration.
log_likelihood = numpy.reshape(log_likelihood,
(len(betas),
len(log_likelihood[0].flatten())))
average_logl = numpy.average(log_likelihood, axis=1)
if method in ("trapezoid", "trapezoid_corrected"):
log_evidence = trapezoid(average_logl, betas)
if method == "trapezoid_corrected":
# var_correction holds the derivative correction terms
# See Friel et al. 2014 for expression and derivation.
# https://link.springer.com/article/10.1007/s11222-013-9397-1
var_correction = 0
for i in range(len(betas) - 1):
delta_beta = betas[i+1] - betas[i]
pre_fac_var = (1. / 12.) * (delta_beta ** 2.0)
var_diff = numpy.var(log_likelihood[i+1])
var_diff -= numpy.var(log_likelihood[i])
var_correction -= pre_fac_var * var_diff
# Add the derivative correction term back to the log_evidence
# from the first if statement.
log_evidence += var_correction
elif method == "simpsons":
# beta -> 0 tends to contribute the least to the integral
# so we can sacrifice precision there, rather than near
# beta -> 1. Option even="last" puts trapezoid rule at
# first few points.
log_evidence = integrate.simps(average_logl, betas,
even="last")
# Estimate the Monte Carlo variance of the evidence calculation
# See (Evans, Annis, 2019.)
# https://www.sciencedirect.com/science/article/pii/S0022249617302651
ti_vec = numpy.zeros(len(log_likelihood[0]))
# Get log likelihood chains by sample and not by temperature.
logl_per_samp = []
for i, _ in enumerate(log_likelihood[0]):
logl_per_samp.append([log_likelihood[x][i] for x in range(len(betas))])
if method in ("trapezoid", "trapezoid_corrected"):
for i, _ in enumerate(log_likelihood[0]):
ti_vec[i] = trapezoid(logl_per_samp[i], betas)
elif method == "simpsons":
for i, _ in enumerate(log_likelihood[0]):
ti_vec[i] = integrate.simps(logl_per_samp[i], betas,
even="last")
# Standard error is sample std / sqrt(number of samples)
mcmc_std = numpy.std(ti_vec) / numpy.sqrt(float(len(log_likelihood[0])))
return log_evidence, mcmc_std
[docs]
def stepping_stone_algorithm(log_likelihood, betas):
"""Returns the log evidence of the model via stepping stone algorithm.
Also returns an estimated standard deviation for the log evidence.
Parameters
----------
log_likelihood : numpy.ndarray
3d array of shape (ntemps, nwalkers, niterations)
The log likelihood for each temperature separated by
temperature, walker, and iteration.
betas : numpy.ndarray
The inverse temperatures used in the MCMC.
Returns
-------
log_evidence : float
Estimation of the log of the evidence.
mcmc_std : float
The standard deviation of the log evidence estimate from
Monte-Carlo spread.
"""
# Reverse order sort
order = numpy.argsort(betas)[::-1]
betas = betas[order]
log_likelihood = log_likelihood[order]
# Assume log likelihood is given in shape of beta,
# walker, iteration.
log_likelihood = numpy.reshape(log_likelihood,
(len(betas),
len(log_likelihood[0].flatten())))
log_rk_pb = numpy.zeros(len(betas) - 1)
for i in range(len(betas) - 1):
delta_beta = betas[i] - betas[i+1]
# Max log likelihood for beta [i+1]
max_logl_pb = numpy.max(log_likelihood[i+1])
val_1 = delta_beta * max_logl_pb
val_2 = delta_beta * (log_likelihood[i+1] - max_logl_pb)
val_2 = numpy.log(numpy.average(numpy.exp(val_2)))
log_rk_pb[i] = val_1 + val_2
log_rk = numpy.sum(log_rk_pb)
log_evidence = log_rk
# Calculate the Monte Carlo variation
mcmc_std = 0
for i in range(len(betas) - 1):
delta_beta = betas[i] - betas[i+1]
pre_fact = (delta_beta * log_likelihood[i+1]) - log_rk_pb[i]
pre_fact = numpy.exp(pre_fact) - 1.0
val = numpy.sum(pre_fact ** 2)
mcmc_std += val
mcmc_std /= float(len(log_likelihood[0])) ** 2.0
mcmc_std = numpy.sqrt(mcmc_std)
return log_evidence, mcmc_std