Part II: Trans-dimensional Modeling
In this section, we discretize the domain of the sought piecewise function through a trans-dimensional 1-D Voronoi tessellation. This approach allows us to treat the number of segments in the piecewise function as unknown.
Import libraries and define constants
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from bayesbay.discretization import Voronoi1D
from bayesbay.prior import UniformPrior
from bayesbay.parameterization import Parameterization
from bayesbay import Target, LogLikelihood, BayesianInversion
np.random.seed(30)
X_DATA = np.linspace(0, 10, 100) # 100 data points
NOISE_STD = 5
@np.vectorize
def piecewise_function(x):
if x <= 1:
return 1
elif 1 < x <= 2.5:
return 20
elif 2.5 < x <= 3:
return 0
elif 3 < x <= 4:
return -3
elif 4 < x <= 6:
return -10
elif 6 < x <= 6.5:
return -20
elif 6.5 < x <= 8:
return 25
elif 8 < x <= 9:
return 0
elif 9 < x <= 10:
return 10
Observed data
d_obs = piecewise_function(X_DATA) + np.random.normal(0, NOISE_STD, X_DATA.size)
fig, ax = plt.subplots(figsize=(10, 6))
ax.step(X_DATA, piecewise_function(X_DATA), 'gray', lw=3, label='True model')
ax.plot(X_DATA, d_obs, 'ro', markeredgecolor='k', label='Obs. data')
ax.legend()
ax.grid()
plt.show()
Inference parameterization
y = UniformPrior('y', vmin=-35, vmax=35, perturb_std=3.5)
voronoi = Voronoi1D(
name="voronoi",
vmin=0,
vmax=10,
perturb_std=0.75,
n_dimensions=None,
n_dimensions_min=2,
n_dimensions_max=40,
parameters=[y],
birth_from='prior'
)
parameterization = Parameterization(voronoi)
Forward and Log Likelihood
def fwd_function(state):
voro = state['voronoi']
x_nuclei = voro['discretization']
x_extents = Voronoi1D.compute_cell_extents(x_nuclei)
x1 = x_extents[0]
indexes = []
indexes.extend([0] * np.flatnonzero(X_DATA <= x1).size)
for i, extent in enumerate(x_extents[1:], 1):
if extent == 0:
idx_size = np.flatnonzero(X_DATA > x1).size
else:
x2 = x1 + x_extents[i]
idx_size = np.flatnonzero((X_DATA > x1) & (X_DATA <= x2)).size
x1 = x2
indexes.extend([i] * idx_size)
d_pred = voro['y'][indexes]
return d_pred
target = Target("d_obs",
d_obs,
std_min=0,
std_max=20,
std_perturb_std=2)
log_likelihood = LogLikelihood(targets=target, fwd_functions=fwd_function)
Bayesian Inference
inversion = BayesianInversion(
parameterization=parameterization,
log_likelihood=log_likelihood,
n_chains=20
)
inversion.run(
sampler=None,
n_iterations=400_000,
burnin_iterations=50_000,
save_every=40,
verbose=False
)
for chain in inversion.chains:
chain.print_statistics()
Chain ID: 0
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 118583/400000 (29.65 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3968/56793 (6.99%)
DeathPerturbation(voronoi): 3963/57442 (6.90%)
NoisePerturbation(d_obs): 12926/57028 (22.67%)
ParamPerturbation(voronoi.discretization): 8612/57154 (15.07%)
ParamPerturbation(voronoi.y): 89114/171583 (51.94%)
Chain ID: 1
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 120724/400000 (30.18 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 4272/57112 (7.48%)
DeathPerturbation(voronoi): 4268/56877 (7.50%)
NoisePerturbation(d_obs): 13128/57034 (23.02%)
ParamPerturbation(voronoi.discretization): 8906/57257 (15.55%)
ParamPerturbation(voronoi.y): 90150/171720 (52.50%)
Chain ID: 2
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 120321/400000 (30.08 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 4424/57297 (7.72%)
DeathPerturbation(voronoi): 4414/57590 (7.66%)
NoisePerturbation(d_obs): 12873/56661 (22.72%)
ParamPerturbation(voronoi.discretization): 8443/57497 (14.68%)
ParamPerturbation(voronoi.y): 90167/170955 (52.74%)
Chain ID: 3
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 115114/400000 (28.78 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3581/57452 (6.23%)
DeathPerturbation(voronoi): 3580/57442 (6.23%)
NoisePerturbation(d_obs): 12601/56914 (22.14%)
ParamPerturbation(voronoi.discretization): 8029/57195 (14.04%)
ParamPerturbation(voronoi.y): 87323/170997 (51.07%)
Chain ID: 4
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 120892/400000 (30.22 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 4330/57341 (7.55%)
DeathPerturbation(voronoi): 4332/57305 (7.56%)
NoisePerturbation(d_obs): 13043/57306 (22.76%)
ParamPerturbation(voronoi.discretization): 8907/57110 (15.60%)
ParamPerturbation(voronoi.y): 90280/170938 (52.81%)
Chain ID: 5
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 117688/400000 (29.42 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3806/56983 (6.68%)
DeathPerturbation(voronoi): 3807/56770 (6.71%)
NoisePerturbation(d_obs): 12972/57316 (22.63%)
ParamPerturbation(voronoi.discretization): 8301/56964 (14.57%)
ParamPerturbation(voronoi.y): 88802/171967 (51.64%)
Chain ID: 6
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 118733/400000 (29.68 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 4039/57463 (7.03%)
DeathPerturbation(voronoi): 4039/57116 (7.07%)
NoisePerturbation(d_obs): 13072/57042 (22.92%)
ParamPerturbation(voronoi.discretization): 8363/57199 (14.62%)
ParamPerturbation(voronoi.y): 89220/171180 (52.12%)
Chain ID: 7
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 119023/400000 (29.76 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3823/56872 (6.72%)
DeathPerturbation(voronoi): 3812/57199 (6.66%)
NoisePerturbation(d_obs): 13932/57662 (24.16%)
ParamPerturbation(voronoi.discretization): 8713/57156 (15.24%)
ParamPerturbation(voronoi.y): 88743/171111 (51.86%)
Chain ID: 8
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 119750/400000 (29.94 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3860/57147 (6.75%)
DeathPerturbation(voronoi): 3854/57152 (6.74%)
NoisePerturbation(d_obs): 13688/57461 (23.82%)
ParamPerturbation(voronoi.discretization): 8691/56999 (15.25%)
ParamPerturbation(voronoi.y): 89657/171241 (52.36%)
Chain ID: 9
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 115797/400000 (28.95 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3468/57252 (6.06%)
DeathPerturbation(voronoi): 3466/57124 (6.07%)
NoisePerturbation(d_obs): 13151/57113 (23.03%)
ParamPerturbation(voronoi.discretization): 8184/57566 (14.22%)
ParamPerturbation(voronoi.y): 87528/170945 (51.20%)
Chain ID: 10
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 114926/400000 (28.73 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3482/57347 (6.07%)
DeathPerturbation(voronoi): 3475/57055 (6.09%)
NoisePerturbation(d_obs): 12852/57407 (22.39%)
ParamPerturbation(voronoi.discretization): 7947/57032 (13.93%)
ParamPerturbation(voronoi.y): 87170/171159 (50.93%)
Chain ID: 11
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 116873/400000 (29.22 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3734/56715 (6.58%)
DeathPerturbation(voronoi): 3730/56988 (6.55%)
NoisePerturbation(d_obs): 12932/56963 (22.70%)
ParamPerturbation(voronoi.discretization): 8099/57340 (14.12%)
ParamPerturbation(voronoi.y): 88378/171994 (51.38%)
Chain ID: 12
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 115765/400000 (28.94 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3612/57074 (6.33%)
DeathPerturbation(voronoi): 3610/57020 (6.33%)
NoisePerturbation(d_obs): 13002/57100 (22.77%)
ParamPerturbation(voronoi.discretization): 8127/56956 (14.27%)
ParamPerturbation(voronoi.y): 87414/171850 (50.87%)
Chain ID: 13
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 118379/400000 (29.59 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3853/56860 (6.78%)
DeathPerturbation(voronoi): 3856/57189 (6.74%)
NoisePerturbation(d_obs): 13166/57171 (23.03%)
ParamPerturbation(voronoi.discretization): 8653/57350 (15.09%)
ParamPerturbation(voronoi.y): 88851/171430 (51.83%)
Chain ID: 14
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 120377/400000 (30.09 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 4006/57082 (7.02%)
DeathPerturbation(voronoi): 4000/56921 (7.03%)
NoisePerturbation(d_obs): 13422/57399 (23.38%)
ParamPerturbation(voronoi.discretization): 8938/57484 (15.55%)
ParamPerturbation(voronoi.y): 90011/171114 (52.60%)
Chain ID: 15
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 115726/400000 (28.93 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3652/57076 (6.40%)
DeathPerturbation(voronoi): 3651/56984 (6.41%)
NoisePerturbation(d_obs): 12954/57009 (22.72%)
ParamPerturbation(voronoi.discretization): 8277/57349 (14.43%)
ParamPerturbation(voronoi.y): 87192/171582 (50.82%)
Chain ID: 16
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 115693/400000 (28.92 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3508/57124 (6.14%)
DeathPerturbation(voronoi): 3502/57076 (6.14%)
NoisePerturbation(d_obs): 13320/57379 (23.21%)
ParamPerturbation(voronoi.discretization): 8113/57198 (14.18%)
ParamPerturbation(voronoi.y): 87250/171223 (50.96%)
Chain ID: 17
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 119872/400000 (29.97 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 4153/57012 (7.28%)
DeathPerturbation(voronoi): 4153/57117 (7.27%)
NoisePerturbation(d_obs): 13036/57231 (22.78%)
ParamPerturbation(voronoi.discretization): 8656/57354 (15.09%)
ParamPerturbation(voronoi.y): 89874/171286 (52.47%)
Chain ID: 18
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 116698/400000 (29.17 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3559/57163 (6.23%)
DeathPerturbation(voronoi): 3556/56850 (6.26%)
NoisePerturbation(d_obs): 13328/57319 (23.25%)
ParamPerturbation(voronoi.discretization): 8204/56948 (14.41%)
ParamPerturbation(voronoi.y): 88051/171720 (51.28%)
Chain ID: 19
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 116328/400000 (29.08 %)
PARTIAL ACCEPTANCE RATES:
BirthPerturbation(voronoi): 3532/57066 (6.19%)
DeathPerturbation(voronoi): 3530/57417 (6.15%)
NoisePerturbation(d_obs): 13133/57395 (22.88%)
ParamPerturbation(voronoi.discretization): 8044/57256 (14.05%)
ParamPerturbation(voronoi.y): 88089/170866 (51.55%)
Results
results = inversion.get_results(concatenate_chains=True)
percentiles = 10, 90
statistics = Voronoi1D.get_tessellation_statistics(
results['voronoi.discretization'],
results['voronoi.y'],
X_DATA,
percentiles=percentiles
)
fig = plt.figure(figsize=(10, 8))
gs = gridspec.GridSpec(2, 2, height_ratios=[1.5, 1])
ax1 = plt.subplot(gs[0, :]) # First row, spans all columns
ax2 = plt.subplot(gs[1, 0]) # Second row, first column
ax3 = plt.subplot(gs[1, 1]) # Second row, second column
ax1.step(X_DATA, piecewise_function(X_DATA), 'gray', lw=3, label='True model')
ax1.plot(X_DATA, statistics['median'], 'b', label='Ensemble median')
ax1.fill_between(X_DATA,
*statistics['percentiles'],
color='b',
alpha=0.2,
label='Uncertainty (%s-%sth perc.)'%(percentiles))
ax1.plot(X_DATA, d_obs, 'ro', markeredgecolor='k', label='Obs. data')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.legend(framealpha=0.5)
ax1.grid()
ax1.set_title('Inferred model')
ndim_min, ndim_max = voronoi._n_dimensions_min, voronoi._n_dimensions_max
ax2.fill_between([ndim_min, ndim_max],
1 / (ndim_max - ndim_min),
alpha=0.3,
color='orange',
label='Prior')
ax2.set_xlabel('No. partitions')
ax2.hist(results['voronoi.n_dimensions'],
bins=np.arange(ndim_min-0.5, ndim_max+0.5),
color='b',
alpha=0.8,
density=True,
ec='w',
label='Posterior')
ax2.axvline(x=9, color='r', lw=2, alpha=1, label='True', zorder=100)
ax2.legend(framealpha=0)
ax2.set_title('Partitions')
ax2.legend(framealpha=0.9)
ax3.fill_between([target.std_min, target.std_max],
1 / (target.std_max - target.std_min),
alpha=0.3,
color='orange',
label='Prior')
ax3.hist(results['d_obs.std'],
color='b',
alpha=0.8,
density=True,
bins=10,
ec='w',
label='Posterior')
ax3.axvline(x=NOISE_STD, color='r', lw=2, alpha=1, label='True', zorder=100)
ax3.set_xlabel('Noise standard deviation')
ax3.legend(framealpha=0.9)
ax3.set_title('Data noise')
plt.tight_layout()
plt.show()
