Skip to content

Commit 96d0eeb

Browse files
author
Junpeng Lao
committed
update examples
* add example into pm.Bound warning, also improve pm.Bound docstring
1 parent 39ed379 commit 96d0eeb

File tree

9 files changed

+60
-85
lines changed

9 files changed

+60
-85
lines changed

pymc3/distributions/distribution.py

Lines changed: 10 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -476,8 +476,14 @@ class Bound(object):
476476
477477
Example
478478
-------
479-
boundedNormal = pymc3.Bound(pymc3.Normal, lower=0.0)
480-
par = boundedNormal(mu=0.0, sd=1.0, testval=1.0)
479+
# In general Bounded distribution need to be defined befor model context
480+
boundedNormal = pm.Bound(pm.Normal, lower=0.0)
481+
with pm.Model():
482+
par1 = boundedNormal('par1', mu=0.0, sd=1.0, testval=1.0)
483+
484+
# or you can define it implicitly within the model context
485+
par2 = pm.Bound(pm.Normal, lower=0.0)(
486+
'par2', mu=0.0, sd=1.0, testval=1.0)
481487
"""
482488

483489
def __init__(self, distribution, lower=-np.inf, upper=np.inf):
@@ -490,7 +496,8 @@ def __call__(self, *args, **kwargs):
490496
raise ValueError('Observed Bound distributions are not allowed. '
491497
'If you want to model truncated data '
492498
'you can use a pm.Potential in combination '
493-
'with the cumulative probability function.')
499+
'with the cumulative probability function. See '
500+
'pymc3/examples/censored_data.py for an example.')
494501
first, args = args[0], args[1:]
495502

496503
return Bounded(first, self.distribution, self.lower, self.upper,

pymc3/examples/baseball.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ def build_model():
1010
data = np.loadtxt(pm.get_data('efron-morris-75-data.tsv'), delimiter="\t",
1111
skiprows=1, usecols=(2,3))
1212

13-
atBats = pm.floatX(data[:,0])
13+
atbats = pm.floatX(data[:,0])
1414
hits = pm.floatX(data[:,1])
1515

1616
N = len(hits)
@@ -22,7 +22,7 @@ def build_model():
2222
phi = pm.Uniform('phi', lower=0.0, upper=1.0)
2323
kappa = BoundedKappa('kappa', alpha=1.0001, m=1.5)
2424
thetas = pm.Beta('thetas', alpha=phi*kappa, beta=(1.0-phi)*kappa, shape=N)
25-
ys = pm.Binomial('ys', n=atBats, p=thetas, observed=hits)
25+
ys = pm.Binomial('ys', n=atbats, p=thetas, observed=hits)
2626
return model
2727

2828
def run(n=2000):

pymc3/examples/censored_data.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ def normal_lcdf(mu, sigma, x):
3838
z = (x - mu) / sigma
3939
return tt.switch(
4040
tt.lt(z, -1.0),
41-
tt.log(tt.erfcx(-z / tt.sqrt(2.)) / 2.) - tt.sqr(z) / 2,
41+
tt.log(tt.erfcx(-z / tt.sqrt(2.)) / 2.) - tt.sqr(z) / 2.,
4242
tt.log1p(-tt.erfc(z / tt.sqrt(2.)) / 2.)
4343
)
4444

@@ -47,7 +47,7 @@ def normal_lccdf(mu, sigma, x):
4747
z = (x - mu) / sigma
4848
return tt.switch(
4949
tt.gt(z, 1.0),
50-
tt.log(tt.erfcx(z / tt.sqrt(2.)) / 2) - tt.sqr(z) / 2.,
50+
tt.log(tt.erfcx(z / tt.sqrt(2.)) / 2.) - tt.sqr(z) / 2.,
5151
tt.log1p(-tt.erfc(-z / tt.sqrt(2.)) / 2.)
5252
)
5353

pymc3/examples/custom_dists.py

Lines changed: 14 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@
1111
import matplotlib.pyplot as plt
1212
import numpy as np
1313

14-
import pymc3
14+
import pymc3 as pm
1515
import theano.tensor as tt
1616

1717
np.random.seed(42)
@@ -24,20 +24,23 @@
2424
ydata = np.random.normal(ydata, 10)
2525
data = {'x': xdata, 'y': ydata}
2626

27-
with pymc3.Model() as model:
28-
alpha = pymc3.Uniform('intercept', -100, 100)
27+
# define loglikelihood outside of the model context, otherwise njobs wont work:
28+
# Lambdas defined in local namespace are not picklable (see issue #1995)
29+
def loglike1(value):
30+
return -1.5 * tt.log(1 + value**2)
31+
def loglike2(value):
32+
return -tt.log(tt.abs_(value))
33+
34+
with pm.Model() as model:
35+
alpha = pm.Normal('intercept', mu=0, sd=100)
2936
# Create custom densities
30-
beta = pymc3.DensityDist('slope', lambda value: -
31-
1.5 * tt.log(1 + value**2), testval=0)
32-
sigma = pymc3.DensityDist(
33-
'sigma', lambda value: -tt.log(tt.abs_(value)), testval=1)
37+
beta = pm.DensityDist('slope', loglike1, testval=0)
38+
sigma = pm.DensityDist('sigma', loglike2, testval=1)
3439
# Create likelihood
35-
like = pymc3.Normal('y_est', mu=alpha + beta *
40+
like = pm.Normal('y_est', mu=alpha + beta *
3641
xdata, sd=sigma, observed=ydata)
3742

38-
start = pymc3.find_MAP()
39-
step = pymc3.NUTS(scaling=start) # Instantiate sampler
40-
trace = pymc3.sample(10000, step, start=start)
43+
trace = pm.sample(2000, njobs=2)
4144

4245

4346
#################################################

pymc3/examples/disaster_model_arbitrary_deterministic.py

Lines changed: 8 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,7 @@
77

88
import pymc3 as pm
99
import theano.tensor as tt
10-
from theano import as_op
11-
from numpy import arange, array, empty
10+
from numpy import arange, array
1211

1312
__all__ = ['disasters_data', 'switchpoint', 'early_mean', 'late_mean', 'rate',
1413
'disasters']
@@ -23,18 +22,6 @@
2322
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
2423
years = len(disasters_data)
2524

26-
# here is the trick
27-
28-
29-
@as_op(itypes=[tt.lscalar, tt.dscalar, tt.dscalar], otypes=[tt.dvector])
30-
def rateFunc(switchpoint, early_mean, late_mean):
31-
''' Concatenate Poisson means '''
32-
out = empty(years)
33-
out[:switchpoint] = early_mean
34-
out[switchpoint:] = late_mean
35-
return out
36-
37-
3825
with pm.Model() as model:
3926

4027
# Prior for distribution of switchpoint location
@@ -46,23 +33,18 @@ def rateFunc(switchpoint, early_mean, late_mean):
4633
# Allocate appropriate Poisson rates to years before and after current
4734
# switchpoint location
4835
idx = arange(years)
49-
# theano style:
50-
# rate = switch(switchpoint >= idx, early_mean, late_mean)
51-
# non-theano style
52-
rate = rateFunc(switchpoint, early_mean, late_mean)
53-
36+
rate = tt.switch(switchpoint >= idx, early_mean, late_mean)
37+
5438
# Data likelihood
5539
disasters = pm.Poisson('disasters', rate, observed=disasters_data)
5640

57-
# Initial values for stochastic nodes
58-
start = {'early_mean': 2., 'late_mean': 3.}
59-
6041
# Use slice sampler for means
6142
step1 = pm.Slice([early_mean, late_mean])
6243
# Use Metropolis for switchpoint, since it accomodates discrete variables
6344
step2 = pm.Metropolis([switchpoint])
64-
65-
# njobs>1 works only with most recent (mid August 2014) Theano version:
66-
# https://github.com/Theano/Theano/pull/2021
67-
tr = pm.sample(1000, tune=500, start=start, step=[step1, step2], njobs=1)
45+
46+
# Initial values for stochastic nodes
47+
start = {'early_mean': 2., 'late_mean': 3.}
48+
49+
tr = pm.sample(1000, tune=500, start=start, step=[step1, step2], njobs=2)
6850
pm.traceplot(tr)

pymc3/examples/factor_potential.py

Lines changed: 14 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,24 @@
11
import pymc3 as pm
22

3-
with pm.Model() as model:
4-
x = pm.Normal('x', 1, 1)
5-
x2 = pm.Potential('x2', -x ** 2)
3+
"""
4+
You can add an arbitrary factor potential to the model likelihood using
5+
pm.Potential. For example you can added Jacobian Adjustment using pm.Potential
6+
when you do model reparameterization. It's similar to `increment_log_prob` in
7+
STAN.
8+
"""
69

7-
start = model.test_point
8-
h = pm.find_hessian(start)
9-
step = pm.Metropolis(model.vars, h)
10+
def build_model():
11+
with pm.Model() as model:
12+
x = pm.Normal('x', 1, 1)
13+
x2 = pm.Potential('x2', -x ** 2)
14+
return model
1015

11-
12-
def run(n=3000):
16+
def run(n=1000):
17+
model = build_model()
1318
if n == "short":
1419
n = 50
1520
with model:
16-
pm.sample(n, step=step, start=start)
21+
pm.sample(n)
1722

1823
if __name__ == '__main__':
1924
run()

pymc3/examples/garch_example.py

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,13 +38,15 @@ def get_garch_model():
3838
shape = r.shape
3939

4040
with Model() as garch:
41-
alpha1 = Normal('alpha1', mu=np.zeros(shape=shape), sd=np.ones(shape=shape), shape=shape)
41+
alpha1 = Normal('alpha1', mu=np.zeros(shape=shape),
42+
sd=np.ones(shape=shape), shape=shape)
4243
BoundedNormal = Bound(Normal, upper=(1 - alpha1))
4344
beta1 = BoundedNormal('beta1',
4445
mu=np.zeros(shape=shape),
4546
sd=1e6 * np.ones(shape=shape),
4647
shape=shape)
47-
mu = Normal('mu', mu=np.zeros(shape=shape), sd=1e6 * np.ones(shape=shape), shape=shape)
48+
mu = Normal('mu', mu=np.zeros(shape=shape),
49+
sd=1e6 * np.ones(shape=shape), shape=shape)
4850
theta = tt.sqrt(alpha0 + alpha1 * tt.pow(r - mu, 2) +
4951
beta1 * tt.pow(sigma1, 2))
5052
Normal('obs', mu, sd=theta, observed=r)

pymc3/examples/lightspeed_example.py

Lines changed: 1 addition & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -23,21 +23,10 @@
2323
# define likelihood
2424
y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=light_speed)
2525

26-
# inference fitting the model
27-
28-
# I have to use slice because the following command
29-
# trace = pm.sample(5000)
30-
# produce the error
31-
# ValueError: Cannot construct a ufunc with more than 32 operands
32-
# (requested number were: inputs = 51 and outputs = 1)valueerror
33-
3426

3527
def run(n=5000):
3628
with model_1:
37-
xstart = pm.find_MAP()
38-
xstep = pm.Slice()
39-
trace = pm.sample(5000, xstep, start=xstart,
40-
random_seed=123, progressbar=True)
29+
trace = pm.sample(n)
4130

4231
pm.summary(trace)
4332

pymc3/examples/simpletest.py

Lines changed: 5 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -9,31 +9,18 @@
99
data = np.random.normal(size=(2, 20))
1010

1111

12-
model = pm.Model()
13-
14-
with model:
15-
x = pm.Normal('x', mu=.5, tau=2. ** -2, shape=(2, 1))
12+
with pm.Model() as model:
13+
x = pm.Normal('x', mu=.5, sd=2., shape=(2, 1))
1614
z = pm.Beta('z', alpha=10, beta=5.5)
17-
d = pm.Normal('data', mu=x, tau=.75 ** -2, observed=data)
18-
step = pm.NUTS()
15+
d = pm.Normal('data', mu=x, sd=.75, observed=data)
1916

2017

2118
def run(n=1000):
2219
if n == "short":
2320
n = 50
2421
with model:
25-
trace = pm.sample(n, step)
26-
27-
plt.subplot(2, 2, 1)
28-
plt.plot(trace[x][:, 0, 0])
29-
plt.subplot(2, 2, 2)
30-
plt.hist(trace[x][:, 0, 0])
31-
32-
plt.subplot(2, 2, 3)
33-
plt.plot(trace[x][:, 1, 0])
34-
plt.subplot(2, 2, 4)
35-
plt.hist(trace[x][:, 1, 0])
36-
plt.show()
22+
trace = pm.sample(n)
23+
pm.traceplot(trace, varnames=['x'])
3724

3825
if __name__ == '__main__':
3926
run()

0 commit comments

Comments
 (0)