/usr/lib/python3/dist-packages/csb/statistics/ars.py is in python3-csb 1.2.2+dfsg-2ubuntu1.
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 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 | """
Adaptive Rejection Sampling (ARS)
The ARS class generates a single random sample from a
univariate distribution specified by an instance of the
LogProb class, implemented by the user. An instance of
LogProb returns the log of the probability density and
its derivative. The log probability function passed must
be concave.
The user must also supply initial guesses. It is not
essential that these values be very accurate, but performance
will generally depend on their accuracy.
"""
from numpy import exp, log
class Envelope(object):
"""
Envelope function for adaptive rejection sampling.
The envelope defines a piecewise linear upper and lower
bounding function of the concave log-probability.
"""
def __init__(self, x, h, dh):
from numpy import array, inf
self.x = array(x)
self.h = array(h)
self.dh = array(dh)
self.z0 = -inf
self.zk = inf
def z(self):
"""
Support intervals for upper bounding function.
"""
from numpy import concatenate
h = self.h
dh = self.dh
x = self.x
z = (h[1:] - h[:-1] + x[:-1] * dh[:-1] - x[1:] * dh[1:]) / \
(dh[:-1] - dh[1:])
return concatenate(([self.z0], z, [self.zk]))
def u(self, x):
"""
Piecewise linear upper bounding function.
"""
z = self.z()[1:-1]
j = (x > z).sum()
return self.h[j] + self.dh[j] * (x - self.x[j])
def u_max(self):
z = self.z()[1:-1]
return (self.h + self.dh * (z - self.x)).max()
def l(self, x):
"""
Piecewise linear lower bounding function.
"""
from numpy import inf
j = (x > self.x).sum()
if j == 0 or j == len(self.x):
return -inf
else:
j -= 1
return ((self.x[j + 1] - x) * self.h[j] + (x - self.x[j]) * self.h[j + 1]) / \
(self.x[j + 1] - self.x[j])
def insert(self, x, h, dh):
"""
Insert new support point for lower bounding function
(and indirectly for upper bounding function).
"""
from numpy import concatenate
j = (x > self.x).sum()
self.x = concatenate((self.x[:j], [x], self.x[j:]))
self.h = concatenate((self.h[:j], [h], self.h[j:]))
self.dh = concatenate((self.dh[:j], [dh], self.dh[j:]))
def log_masses(self):
from numpy import abs, putmask
z = self.z()
b = self.h - self.x * self.dh
a = abs(self.dh)
m = (self.dh > 0)
q = self.x * 0.
putmask(q, m, z[1:])
putmask(q, 1 - m, z[:-1])
log_M = b - log(a) + log(1 - exp(-a * (z[1:] - z[:-1]))) + \
self.dh * q
return log_M
def masses(self):
z = self.z()
b = self.h - self.x * self.dh
a = self.dh
return exp(b) * (exp(a * z[1:]) - exp(a * z[:-1])) / a
def sample(self):
from numpy.random import random
from numpy import add
from csb.numeric import log_sum_exp
log_m = self.log_masses()
log_M = log_sum_exp(log_m)
c = add.accumulate(exp(log_m - log_M))
u = random()
j = (u > c).sum()
a = self.dh[j]
z = self.z()
xmin, xmax = z[j], z[j + 1]
u = random()
if a > 0:
return xmax + log(u + (1 - u) * exp(-a * (xmax - xmin))) / a
else:
return xmin + log(u + (1 - u) * exp(a * (xmax - xmin))) / a
class LogProb(object):
def __call__(self, x):
raise NotImplementedError()
class Gauss(LogProb):
def __init__(self, mu, sigma=1.):
self.mu = float(mu)
self.sigma = float(sigma)
def __call__(self, x):
return -0.5 * (x - self.mu) ** 2 / self.sigma ** 2, \
- (x - self.mu) / self.sigma ** 2
class ARS(object):
from numpy import inf
def __init__(self, logp):
self.logp = logp
def initialize(self, x, z0=-inf, zmax=inf):
from numpy import array
self.hull = Envelope(array(x), *self.logp(array(x)))
self.hull.z0 = z0
self.hull.zk = zmax
def sample(self, maxiter=100):
from numpy.random import random
for i in range(maxiter):
x = self.hull.sample()
l = self.hull.l(x)
u = self.hull.u(x)
w = random()
if w <= exp(l - u): return x
h, dh = self.logp(x)
if w <= exp(h - u): return x
self.hull.insert(x, h, dh)
|