Part I: Known number of Partitions
In this section of the tutorial, we demonstrate how BayesBay can be utilized to infer the piecewise function given prior knowledge of the number of segments composing it. This is achieved by dividing the domain into nine partitions, defined by a 1-D Voronoi tessellation, while treating the positions of the Voronoi sites as unknown. For a trans-dimensional approach, where the number of partitions is an unknown itself, please refer to the Part II: Trans-dimensional Modeling.
Import libraries and define constants
import numpy as np
import matplotlib.pyplot as plt
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=9,
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: 153966/400000 (38.49 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 19269/80243 (24.01%)
ParamPerturbation(voronoi.discretization): 12175/79861 (15.25%)
ParamPerturbation(voronoi.y): 122522/239896 (51.07%)
Chain ID: 1
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 143310/400000 (35.83 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 18280/79988 (22.85%)
ParamPerturbation(voronoi.discretization): 9919/79973 (12.40%)
ParamPerturbation(voronoi.y): 115111/240039 (47.96%)
Chain ID: 2
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 159921/400000 (39.98 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 20610/80260 (25.68%)
ParamPerturbation(voronoi.discretization): 12862/80084 (16.06%)
ParamPerturbation(voronoi.y): 126449/239656 (52.76%)
Chain ID: 3
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 150881/400000 (37.72 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 19355/80481 (24.05%)
ParamPerturbation(voronoi.discretization): 10625/79992 (13.28%)
ParamPerturbation(voronoi.y): 120901/239527 (50.47%)
Chain ID: 4
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 159046/400000 (39.76 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 20601/79811 (25.81%)
ParamPerturbation(voronoi.discretization): 11366/80018 (14.20%)
ParamPerturbation(voronoi.y): 127079/240171 (52.91%)
Chain ID: 5
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 154239/400000 (38.56 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 19358/80072 (24.18%)
ParamPerturbation(voronoi.discretization): 11810/80058 (14.75%)
ParamPerturbation(voronoi.y): 123071/239870 (51.31%)
Chain ID: 6
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 181101/400000 (45.28 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 23905/79877 (29.93%)
ParamPerturbation(voronoi.discretization): 16126/79799 (20.21%)
ParamPerturbation(voronoi.y): 141070/240324 (58.70%)
Chain ID: 7
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 169359/400000 (42.34 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 21456/79912 (26.85%)
ParamPerturbation(voronoi.discretization): 14087/79997 (17.61%)
ParamPerturbation(voronoi.y): 133816/240091 (55.74%)
Chain ID: 8
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 157328/400000 (39.33 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 20149/79854 (25.23%)
ParamPerturbation(voronoi.discretization): 11834/80028 (14.79%)
ParamPerturbation(voronoi.y): 125345/240118 (52.20%)
Chain ID: 9
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 186764/400000 (46.69 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 24787/80076 (30.95%)
ParamPerturbation(voronoi.discretization): 16562/79654 (20.79%)
ParamPerturbation(voronoi.y): 145415/240270 (60.52%)
Chain ID: 10
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 158838/400000 (39.71 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 20573/80036 (25.70%)
ParamPerturbation(voronoi.discretization): 11640/80183 (14.52%)
ParamPerturbation(voronoi.y): 126625/239781 (52.81%)
Chain ID: 11
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 158359/400000 (39.59 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 20340/79590 (25.56%)
ParamPerturbation(voronoi.discretization): 11011/80016 (13.76%)
ParamPerturbation(voronoi.y): 127008/240394 (52.83%)
Chain ID: 12
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 153008/400000 (38.25 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 19364/79847 (24.25%)
ParamPerturbation(voronoi.discretization): 12118/79829 (15.18%)
ParamPerturbation(voronoi.y): 121526/240324 (50.57%)
Chain ID: 13
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 167548/400000 (41.89 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 21406/79910 (26.79%)
ParamPerturbation(voronoi.discretization): 14232/79884 (17.82%)
ParamPerturbation(voronoi.y): 131910/240206 (54.92%)
Chain ID: 14
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 145866/400000 (36.47 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 18482/80208 (23.04%)
ParamPerturbation(voronoi.discretization): 10681/79856 (13.38%)
ParamPerturbation(voronoi.y): 116703/239936 (48.64%)
Chain ID: 15
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 154991/400000 (38.75 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 19356/80119 (24.16%)
ParamPerturbation(voronoi.discretization): 11800/79897 (14.77%)
ParamPerturbation(voronoi.y): 123835/239984 (51.60%)
Chain ID: 16
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 159376/400000 (39.84 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 20595/79893 (25.78%)
ParamPerturbation(voronoi.discretization): 12060/80299 (15.02%)
ParamPerturbation(voronoi.y): 126721/239808 (52.84%)
Chain ID: 17
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 145072/400000 (36.27 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 18355/79477 (23.09%)
ParamPerturbation(voronoi.discretization): 10184/79967 (12.74%)
ParamPerturbation(voronoi.y): 116533/240556 (48.44%)
Chain ID: 18
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 155754/400000 (38.94 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 19762/79960 (24.71%)
ParamPerturbation(voronoi.discretization): 11494/80245 (14.32%)
ParamPerturbation(voronoi.y): 124498/239795 (51.92%)
Chain ID: 19
TEMPERATURE: 1
EXPLORED MODELS: 400000
ACCEPTANCE RATE: 162794/400000 (40.70 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(d_obs): 21207/80118 (26.47%)
ParamPerturbation(voronoi.discretization): 12100/80006 (15.12%)
ParamPerturbation(voronoi.y): 129487/239876 (53.98%)
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, (ax1, ax2) = plt.subplots(1,2, figsize=(16, 6), gridspec_kw={'width_ratios': [1.5, 1]})
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')
ax2.fill_between([target.std_min, target.std_max],
1 / (target.std_max - target.std_min),
alpha=0.3,
color='orange',
label='Prior')
ax2.hist(results['d_obs.std'],
color='b',
alpha=0.8,
density=True,
bins=10,
ec='w',
label='Posterior')
ax2.axvline(x=NOISE_STD, color='r', lw=2, alpha=1, label='True', zorder=100)
ax2.set_xlabel('Noise standard deviation')
ax2.legend(framealpha=0.9)
ax2.set_title('Data noise')
plt.tight_layout()
plt.show()
