Skip to content

Commit af2e68f

Browse files
committed
Fix sample_ppc
1 parent a764d6e commit af2e68f

File tree

7 files changed

+210
-71
lines changed

7 files changed

+210
-71
lines changed

pymc3/distributions/continuous.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -298,7 +298,7 @@ def __init__(self, sd=None, tau=None, *args, **kwargs):
298298
assert_negative_support(sd, 'sd', 'HalfNormal')
299299

300300
def random(self, point=None, size=None, repeat=None):
301-
sd = draw_values([self.sd], point=point)
301+
sd = draw_values([self.sd], point=point)[0]
302302
return generate_samples(stats.halfnorm.rvs, loc=0., scale=sd,
303303
dist_shape=self.shape,
304304
size=size)
@@ -578,7 +578,7 @@ def __init__(self, lam, *args, **kwargs):
578578
assert_negative_support(lam, 'lam', 'Exponential')
579579

580580
def random(self, point=None, size=None, repeat=None):
581-
lam = draw_values([self.lam], point=point)
581+
lam = draw_values([self.lam], point=point)[0]
582582
return generate_samples(np.random.exponential, scale=1. / lam,
583583
dist_shape=self.shape,
584584
size=size)
@@ -962,7 +962,7 @@ def _random(self, beta, size=None):
962962
return beta * np.abs(np.tan(np.pi * (u - 0.5)))
963963

964964
def random(self, point=None, size=None, repeat=None):
965-
beta = draw_values([self.beta], point=point)
965+
beta = draw_values([self.beta], point=point)[0]
966966
return generate_samples(self._random, beta,
967967
dist_shape=self.shape,
968968
size=size)

pymc3/distributions/discrete.py

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,7 @@ def __init__(self, p, *args, **kwargs):
170170
self.mode = tt.cast(tround(p), 'int8')
171171

172172
def random(self, point=None, size=None, repeat=None):
173-
p = draw_values([self.p], point=point)
173+
p = draw_values([self.p], point=point)[0]
174174
return generate_samples(stats.bernoulli.rvs, p,
175175
dist_shape=self.shape,
176176
size=size)
@@ -286,7 +286,7 @@ def __init__(self, mu, *args, **kwargs):
286286
self.mode = tt.floor(mu).astype('int32')
287287

288288
def random(self, point=None, size=None, repeat=None):
289-
mu = draw_values([self.mu], point=point)
289+
mu = draw_values([self.mu], point=point)[0]
290290
return generate_samples(stats.poisson.rvs, mu,
291291
dist_shape=self.shape,
292292
size=size)
@@ -398,7 +398,7 @@ def __init__(self, p, *args, **kwargs):
398398
self.mode = 1
399399

400400
def random(self, point=None, size=None, repeat=None):
401-
p = draw_values([self.p], point=point)
401+
p = draw_values([self.p], point=point)[0]
402402
return generate_samples(np.random.geometric, p,
403403
dist_shape=self.shape,
404404
size=size)
@@ -559,7 +559,7 @@ def __init__(self, c, *args, **kwargs):
559559
self.mean = self.median = self.mode = self.c = c = tt.as_tensor_variable(c)
560560

561561
def random(self, point=None, size=None, repeat=None):
562-
c = draw_values([self.c], point=point)
562+
c = draw_values([self.c], point=point)[0]
563563
dtype = np.array(c).dtype
564564

565565
def _random(c, dtype=dtype, size=None):

pymc3/distributions/distribution.py

Lines changed: 56 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
import numbers
12
import numpy as np
23
import theano.tensor as tt
34
from theano import function
@@ -48,7 +49,8 @@ def dist(cls, *args, **kwargs):
4849
dist.__init__(*args, **kwargs)
4950
return dist
5051

51-
def __init__(self, shape, dtype, testval=None, defaults=(), transform=None, broadcastable=None):
52+
def __init__(self, shape, dtype, testval=None, defaults=(),
53+
transform=None, broadcastable=None):
5254
self.shape = np.atleast_1d(shape)
5355
if False in (np.floor(self.shape) == self.shape):
5456
raise TypeError("Expected int elements in shape")
@@ -70,8 +72,10 @@ def get_test_val(self, val, defaults):
7072
return self.getattr_value(val)
7173

7274
if val is None:
73-
raise AttributeError(str(self) + " has no finite default value to use, checked: " +
74-
str(defaults) + " pass testval argument or adjust so value is finite.")
75+
raise AttributeError("%s has no finite default value to use, "
76+
"checked: %s. Pass testval argument or "
77+
"adjust so value is finite."
78+
% (self, str(defaults)))
7579

7680
def getattr_value(self, val):
7781
if isinstance(val, string_types):
@@ -97,7 +101,8 @@ def TensorType(dtype, shape, broadcastable=None):
97101

98102
class NoDistribution(Distribution):
99103

100-
def __init__(self, shape, dtype, testval=None, defaults=(), transform=None, parent_dist=None, *args, **kwargs):
104+
def __init__(self, shape, dtype, testval=None, defaults=(),
105+
transform=None, parent_dist=None, *args, **kwargs):
101106
super(NoDistribution, self).__init__(shape=shape, dtype=dtype,
102107
testval=testval, defaults=defaults,
103108
*args, **kwargs)
@@ -137,7 +142,8 @@ def __init__(self, shape=(), dtype=None, defaults=('mode',),
137142
class Continuous(Distribution):
138143
"""Base class for continuous distributions"""
139144

140-
def __init__(self, shape=(), dtype=None, defaults=('median', 'mean', 'mode'), *args, **kwargs):
145+
def __init__(self, shape=(), dtype=None, defaults=('median', 'mean', 'mode'),
146+
*args, **kwargs):
141147
if dtype is None:
142148
dtype = theano.config.floatX
143149
super(Continuous, self).__init__(
@@ -190,18 +196,13 @@ def draw_values(params, point=None):
190196
if param.name in named_nodes:
191197
named_nodes.pop(param.name)
192198
for name, node in named_nodes.items():
193-
if not isinstance(node, (tt.sharedvar.TensorSharedVariable,
199+
if not isinstance(node, (tt.sharedvar.SharedVariable,
194200
tt.TensorConstant)):
195-
givens[name] = (node, draw_value(node, point=point))
196-
values = [None for _ in params]
197-
for i, param in enumerate(params):
198-
# "Homogonise" output
199-
values[i] = np.atleast_1d(draw_value(
200-
param, point=point, givens=givens.values()))
201-
if len(values) == 1:
202-
return values[0]
203-
else:
204-
return values
201+
givens[name] = (node, _draw_value(node, point=point))
202+
values = []
203+
for param in params:
204+
values.append(_draw_value(param, point=point, givens=givens.values()))
205+
return values
205206

206207

207208
@memoize
@@ -229,43 +230,45 @@ def _compile_theano_function(param, vars, givens=None):
229230
allow_input_downcast=True)
230231

231232

232-
def draw_value(param, point=None, givens=()):
233-
if hasattr(param, 'name'):
234-
if hasattr(param, 'model'):
235-
if point is not None and param.name in point:
236-
value = point[param.name]
237-
elif hasattr(param, 'random') and param.random is not None:
238-
value = param.random(point=point, size=None)
239-
else:
240-
value = param.tag.test_value
241-
else:
242-
input_pairs = ([g[0] for g in givens],
243-
[g[1] for g in givens])
233+
def _draw_value(param, point=None, givens=None):
234+
"""Draw a random value from a distribution or return a constant.
244235
245-
value = _compile_theano_function(param,
246-
input_pairs[0])(*input_pairs[1])
236+
Parameters
237+
----------
238+
param : number, array like, theano variable or pymc3 random variable
239+
The value or distribution. Constants or shared variables
240+
will be converted to an array and returned. Theano variables
241+
are evaluated. If `param` is a pymc3 random variables, draw
242+
a new value from it and return that, unless a value is specified
243+
in `point`.
244+
point : dict, optional
245+
A dictionary from pymc3 variable names to their values.
246+
givens : dict, optional
247+
A dictionary from theano variables to their values. These values
248+
are used to evaluate `param` if it is a theano variable.
249+
"""
250+
if isinstance(param, numbers.Number):
251+
return param
252+
elif isinstance(param, np.ndarray):
253+
return param
254+
elif isinstance(param, tt.TensorConstant):
255+
return param.value
256+
elif isinstance(param, tt.sharedvar.SharedVariable):
257+
return param.get_value()
258+
elif isinstance(param, tt.TensorVariable):
259+
if point and hasattr(param, 'model') and param.name in point:
260+
return point[param.name]
261+
elif hasattr(param, 'random') and param.random is not None:
262+
return param.random(point=point, size=None)
263+
else:
264+
if givens:
265+
variables, values = list(zip(*givens))
266+
else:
267+
variables = values = []
268+
func = _compile_theano_function(param, variables)
269+
return func(*values)
247270
else:
248-
value = param
249-
250-
# Sanitising values may be necessary.
251-
if hasattr(value, 'value'):
252-
value = value.value
253-
elif hasattr(value, 'get_value'):
254-
value = value.get_value()
255-
256-
if hasattr(param, 'dtype'):
257-
value = np.atleast_1d(value).astype(param.dtype)
258-
if hasattr(param, 'shape'):
259-
try:
260-
shape = param.shape.tag.test_value
261-
except AttributeError:
262-
try:
263-
shape = param.shape.eval()
264-
except AttributeError:
265-
shape = param.shape
266-
if len(shape) == 0 and len(value) == 1:
267-
value = value[0]
268-
return value
271+
raise ValueError('Unexpected type in draw_value: %s' % type(param))
269272

270273

271274
def broadcast_shapes(*args):
@@ -490,8 +493,8 @@ class Bound(object):
490493
# or within the model context
491494
NegativeNormal = pm.Bound(pm.Normal, upper=0.0)
492495
par2 = NegativeNormal('par2', mu=0.0, sd=1.0, testval=1.0)
493-
494-
# or you can define it implicitly within the model context
496+
497+
# or you can define it implicitly within the model context
495498
par3 = pm.Bound(pm.Normal, lower=-1.0, upper=1.0)(
496499
'par3', mu=0.0, sd=1.0, testval=1.0)
497500
"""

pymc3/distributions/multivariate.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -410,7 +410,7 @@ def __init__(self, a, transform=transforms.stick_breaking,
410410
np.nan)
411411

412412
def random(self, point=None, size=None):
413-
a = draw_values([self.a], point=point)
413+
a = draw_values([self.a], point=point)[0]
414414

415415
def _random(a, size=None):
416416
return stats.dirichlet.rvs(a, None if size == a.shape else size)

pymc3/step_methods/elliptical_slice.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ def astep(self, q0, logp):
8888

8989
# Draw from the normal prior by multiplying the Cholesky decomposition
9090
# of the covariance with draws from a standard normal
91-
chol = draw_values([self.prior_chol])
91+
chol = draw_values([self.prior_chol])[0]
9292
nu = np.dot(chol, nr.randn(chol.shape[0]))
9393
y = logp(q0) - nr.standard_exponential()
9494

pymc3/tests/test_random.py

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
import pymc3 as pm
2+
import numpy as np
3+
import numpy.testing as npt
4+
import pytest
5+
import theano.tensor as tt
6+
import theano
7+
8+
from pymc3.distributions.distribution import _draw_value, draw_values
9+
10+
11+
def test_draw_value():
12+
npt.assert_equal(_draw_value(np.array([5, 6])), [5, 6])
13+
npt.assert_equal(_draw_value(np.array(5.)), 5)
14+
15+
npt.assert_equal(_draw_value(tt.constant([5., 6.])), [5, 6])
16+
assert _draw_value(tt.constant(5)) == 5
17+
npt.assert_equal(_draw_value(2 * tt.constant([5., 6.])), [10, 12])
18+
19+
val = theano.shared(np.array([5., 6.]))
20+
npt.assert_equal(_draw_value(val), [5, 6])
21+
npt.assert_equal(_draw_value(2 * val), [10, 12])
22+
23+
a = tt.scalar('a')
24+
a.tag.test_value = 6
25+
npt.assert_equal(_draw_value(2 * a, givens=[(a, 1)]), 2)
26+
27+
assert _draw_value(5) == 5
28+
assert _draw_value(5.) == 5
29+
assert isinstance(_draw_value(5.), type(5.))
30+
assert isinstance(_draw_value(5), type(5))
31+
32+
with pm.Model():
33+
mu = 2 * tt.constant(np.array([5., 6.])) + theano.shared(np.array(5))
34+
a = pm.Normal('a', mu=mu, sd=5, shape=2)
35+
36+
val1 = _draw_value(a)
37+
val2 = _draw_value(a)
38+
assert np.all(val1 != val2)
39+
40+
with pytest.raises(ValueError) as err:
41+
_draw_value([])
42+
err.match('Unexpected type')
43+
44+
45+
class TestDrawValues(object):
46+
def test_empty(self):
47+
assert draw_values([]) == []
48+
49+
def test_vals(self):
50+
npt.assert_equal(draw_values([np.array([5, 6])])[0], [5, 6])
51+
npt.assert_equal(draw_values([np.array(5.)])[0], 5)
52+
53+
npt.assert_equal(draw_values([tt.constant([5., 6.])])[0], [5, 6])
54+
assert draw_values([tt.constant(5)])[0] == 5
55+
npt.assert_equal(draw_values([2 * tt.constant([5., 6.])])[0], [10, 12])
56+
57+
val = theano.shared(np.array([5., 6.]))
58+
npt.assert_equal(draw_values([val])[0], [5, 6])
59+
npt.assert_equal(draw_values([2 * val])[0], [10, 12])
60+
61+
def test_simple_model(self):
62+
with pm.Model():
63+
mu = 2 * tt.constant(np.array([5., 6.])) + theano.shared(np.array(5))
64+
a = pm.Normal('a', mu=mu, sd=5, shape=2)
65+
66+
val1 = draw_values([a])
67+
val2 = draw_values([a])
68+
assert np.all(val1[0] != val2[0])
69+
70+
point = {'a': np.array([3., 4.])}
71+
npt.assert_equal(draw_values([a], point=point), [point['a']])
72+
73+
def test_dep_vars(self):
74+
with pm.Model():
75+
mu = 2 * tt.constant(np.array([5., 6.])) + theano.shared(np.array(5))
76+
sd = pm.HalfNormal('sd', shape=2)
77+
tau = 1 / sd ** 2
78+
a = pm.Normal('a', mu=mu, tau=tau, shape=2)
79+
80+
point = {'a': np.array([1., 2.])}
81+
npt.assert_equal(draw_values([a], point=point), [point['a']])
82+
83+
with pytest.raises(theano.gof.MissingInputError):
84+
draw_values([a])
85+
86+
# We need the untransformed vars
87+
with pytest.raises(theano.gof.MissingInputError):
88+
draw_values([a], point={'sd': np.array([2., 3.])})
89+
90+
val1 = draw_values([a], point={'sd_log__': np.array([2., 3.])})[0]
91+
val2 = draw_values([a], point={'sd_log__': np.array([2., 3.])})[0]
92+
assert np.all(val1 != val2)

0 commit comments

Comments
 (0)