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
|
#!/usr/bin/env python
"""
Example model with strong correlations between the fitted parameters.
We use a*x = y + N(0,1) made complicated by defining a=p1+p2.
The expected distribution for p1 and p2 will be uniform, with p2 = a-p1 in
each sample. Because this distribution is inherently unbounded, artificial
bounds are required on a least one of the parameters for finite duration
simulations.
The expected distribution for p1+p2 can be determined from the linear model
y = a*x. This is reported along with the values estimated from MCMC.
"""
from dream import * # sampler functions
from pylab import * # Numeric functions and plotting
# Create the correlation function and generate some fake data
x = linspace(-1.0, 1, 40)
fn = lambda p: sum(p) * x
bounds = (-20, -inf), (40, inf)
sigma = 1
data = fn((1, 1)) + randn(*x.shape) * sigma # Fake data
# Sample from the posterior density function
n = 2
model = Simulation(f=fn, data=data, sigma=sigma, bounds=bounds, labels=["x", "y"])
sampler = Dream(
model=model,
population=randn(5 * n, 4, n),
thinning=1,
draws=20000,
)
mc = sampler.sample()
mc.title = "Strong anti-correlation"
# Create a derived parameter without the correlation
mc.derive_vars(lambda p: (p[0] + p[1]), labels=["x+y"])
# Compare the MCMC estimate for the derived parameter to a least squares fit
from bumps.wsolve import wpolyfit
poly = wpolyfit(x, data, degree=1, origin=True)
print("x+y from linear fit", poly.coeff[0], poly.std[0])
points, logp = mc.sample(portion=0.5)
print("x+y from MCMC", mean(points[:, 2]), std(points[:, 2], ddof=1))
# Plot the samples
plot_all(mc, portion=0.5)
show()
|
