/usr/lib/python3/dist-packages/emcee/mh.py is in python3-emcee 2.1.0-5.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 | #!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
A vanilla Metropolis-Hastings sampler
"""
from __future__ import (division, print_function, absolute_import,
unicode_literals)
__all__ = ["MHSampler"]
import numpy as np
from . import autocorr
from .sampler import Sampler
# === MHSampler ===
class MHSampler(Sampler):
"""
The most basic possible Metropolis-Hastings style MCMC sampler
:param cov:
The covariance matrix to use for the proposal distribution.
:param dim:
Number of dimensions in the parameter space.
:param lnpostfn:
A function that takes a vector in the parameter space as input and
returns the natural logarithm of the posterior probability for that
position.
:param args: (optional)
A list of extra positional arguments for ``lnpostfn``. ``lnpostfn``
will be called with the sequence ``lnpostfn(p, *args, **kwargs)``.
:param kwargs: (optional)
A list of extra keyword arguments for ``lnpostfn``. ``lnpostfn``
will be called with the sequence ``lnpostfn(p, *args, **kwargs)``.
"""
def __init__(self, cov, *args, **kwargs):
super(MHSampler, self).__init__(*args, **kwargs)
self.cov = cov
def reset(self):
super(MHSampler, self).reset()
self._chain = np.empty((0, self.dim))
self._lnprob = np.empty(0)
def sample(self, p0, lnprob=None, randomstate=None, thin=1,
storechain=True, iterations=1):
"""
Advances the chain ``iterations`` steps as an iterator
:param p0:
The initial position vector.
:param lnprob0: (optional)
The log posterior probability at position ``p0``. If ``lnprob``
is not provided, the initial value is calculated.
:param rstate0: (optional)
The state of the random number generator. See the
:func:`random_state` property for details.
:param iterations: (optional)
The number of steps to run.
:param thin: (optional)
If you only want to store and yield every ``thin`` samples in the
chain, set thin to an integer greater than 1.
:param storechain: (optional)
By default, the sampler stores (in memory) the positions and
log-probabilities of the samples in the chain. If you are
using another method to store the samples to a file or if you
don't need to analyse the samples after the fact (for burn-in
for example) set ``storechain`` to ``False``.
At each iteration, this generator yields:
* ``pos`` - The current positions of the chain in the parameter
space.
* ``lnprob`` - The value of the log posterior at ``pos`` .
* ``rstate`` - The current state of the random number generator.
"""
self.random_state = randomstate
p = np.array(p0)
if lnprob is None:
lnprob = self.get_lnprob(p)
# Resize the chain in advance.
if storechain:
N = int(iterations / thin)
self._chain = np.concatenate((self._chain,
np.zeros((N, self.dim))), axis=0)
self._lnprob = np.append(self._lnprob, np.zeros(N))
i0 = self.iterations
# Use range instead of xrange for python 3 compatability
for i in range(int(iterations)):
self.iterations += 1
# Calculate the proposal distribution.
q = self._random.multivariate_normal(p, self.cov)
newlnprob = self.get_lnprob(q)
diff = newlnprob - lnprob
# M-H acceptance ratio
if diff < 0:
diff = np.exp(diff) - self._random.rand()
if diff > 0:
p = q
lnprob = newlnprob
self.naccepted += 1
if storechain and i % thin == 0:
ind = i0 + int(i / thin)
self._chain[ind, :] = p
self._lnprob[ind] = lnprob
# Heavy duty iterator action going on right here...
yield p, lnprob, self.random_state
@property
def acor(self):
"""
An estimate of the autocorrelation time for each parameter (length:
``dim``).
"""
return self.get_autocorr_time()
def get_autocorr_time(self, window=50):
"""
Compute an estimate of the autocorrelation time for each parameter
(length: ``dim``).
:param window: (optional)
The size of the windowing function. This is equivalent to the
maximum number of lags to use. (default: 50)
"""
return autocorr.integrated_time(self.chain, axis=0, window=window)
|