Finding the parameters of a normal distribution

In this tutorial we will use black-it to find the mean and variance of a normal distribution by fitting the 'model' to a dataset. Obviously the parameters of a normal distribution can be obtained more efficiently (such as a by maximum likelihood), yet this example can be useful to understand how black-it works in practice.

# import standard libraries
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

We stard by loading the dataset that we want to fit with our model.

# load a dataset

true_params = [1, 1]  # in general these are not known!
real_data = np.atleast_2d(
    np.loadtxt(f"data/gaussian_mean{true_params[0]}_var{true_params[1]}.txt")
).T

plt.plot(real_data[:, 0])
plt.xlabel("time step")
plt.ylabel("value")

Text(0, 0.5, 'value')

Initialize a calibrator object

Then, to set up a calibration one needs to define the following components first:

1. a model to be calibrated
2. a loss function to measure the distance between the real time series and the simulated time series
3. a set of samplers that iteratively suggest a set of parameter values to explore
4. the parameter space that should be explored

1 Model simulator

# a normal distribution with unknown mean and variance
from models.simple_models import NormalMV

# when called with a set of parameter values, the model provides a simulated time series
NormalMV([1, 1], N=10, seed=0)

array([[2.76405235],
       [1.40015721],
       [1.97873798],
       [3.2408932 ],
       [2.86755799],
       [0.02272212],
       [1.95008842],
       [0.84864279],
       [0.89678115],
       [1.4105985 ]])

2 Loss function

# loss function based on the Method Of Moments
from black_it.loss_functions.msm import MethodOfMomentsLoss

loss = MethodOfMomentsLoss()

3 Samplers

# import some samplers
from black_it.samplers.best_batch import BestBatchSampler
from black_it.samplers.halton import HaltonSampler
from black_it.samplers.r_sequence import RSequenceSampler
from black_it.samplers.random_forest import RandomForestSampler
from black_it.samplers.random_uniform import RandomUniformSampler

# initialize the samplers with their specific 'batch_size', i.e. the number of points
# explored every time they are called
batch_size = 4
random_sampler = RandomUniformSampler(batch_size=batch_size)
rseq_sampler = RSequenceSampler(batch_size=batch_size)
halton_sampler = HaltonSampler(batch_size=batch_size)
best_batch_sampler = BestBatchSampler(batch_size=batch_size)
random_forest_sampler = RandomForestSampler(batch_size=batch_size)

samplers = [
    random_sampler,
    rseq_sampler,
    halton_sampler,
    random_forest_sampler,
    best_batch_sampler,
]

4) Parameter space (bounds and precision)

# the full space of parameters is defined by the lower and upper bounds
# and by the precision of each parameter
bounds = [[0.00, 0.01], [2.00, 2.00]]

precisions = [0.0001, 0.0001]

Finally, initialise a Calibrator object

from black_it.calibrator import Calibrator

# initialize a Calibrator object
cal = Calibrator(
    samplers=samplers,
    real_data=real_data,
    model=NormalMV,
    parameters_bounds=bounds,
    parameters_precision=precisions,
    ensemble_size=5,
    loss_function=loss,
    saving_folder=None,
)

***
Number of free params:       2.
Explorable param space size: 398039901.
***

Selecting 4 processes for the parallel evaluation of the model

Calibration

params, losses = cal.calibrate(10)

BATCH NUMBER:   1
PARAMS SAMPLED: 0

METHOD: RandomUniformSampler
----> sim exec elapsed time: 1.7s
---->   min loss new params: 0.17
---->   avg loss new params: 1.41
----> avg loss exist params: 1.41
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 1.7s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 1.6
---->   avg loss new params: 1.97
----> avg loss exist params: 1.69
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.77
---->   avg loss new params: 1.46
----> avg loss exist params: 1.61
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.2s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.64
---->   avg loss new params: 0.94
----> avg loss exist params: 1.45
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 1.7s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.2s
---->   min loss new params: 0.54
---->   avg loss new params: 0.82
----> avg loss exist params: 1.32
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.2s

BATCH NUMBER:   2
PARAMS SAMPLED: 20

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.65
---->   avg loss new params: 1.76
----> avg loss exist params: 1.39
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.68
---->   avg loss new params: 1.84
----> avg loss exist params: 1.46
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.2s
---->   min loss new params: 0.86
---->   avg loss new params: 1.71
----> avg loss exist params: 1.49
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.2s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.35
---->   avg loss new params: 0.57
----> avg loss exist params: 1.39
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 1.7s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.42
---->   avg loss new params: 0.53
----> avg loss exist params: 1.3
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

BATCH NUMBER:   3
PARAMS SAMPLED: 40

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.96
---->   avg loss new params: 2.14
----> avg loss exist params: 1.38
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.52
---->   avg loss new params: 1.24
----> avg loss exist params: 1.37
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.67
---->   avg loss new params: 1.89
----> avg loss exist params: 1.41
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.34
---->   avg loss new params: 0.57
----> avg loss exist params: 1.35
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 1.7s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.29
---->   avg loss new params: 0.55
----> avg loss exist params: 1.29
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

BATCH NUMBER:   4
PARAMS SAMPLED: 60

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 1.46
---->   avg loss new params: 2.06
----> avg loss exist params: 1.34
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.52
---->   avg loss new params: 1.77
----> avg loss exist params: 1.37
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.92
---->   avg loss new params: 1.45
----> avg loss exist params: 1.37
---->         curr min loss: 0.16632065926583173
====>    total elapsed time: 0.1s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.1
---->   avg loss new params: 0.37
----> avg loss exist params: 1.32
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 1.6s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.23
---->   avg loss new params: 0.29
----> avg loss exist params: 1.27
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

BATCH NUMBER:   5
PARAMS SAMPLED: 80

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.63
---->   avg loss new params: 1.1
----> avg loss exist params: 1.26
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 1.17
---->   avg loss new params: 1.67
----> avg loss exist params: 1.28
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.7
---->   avg loss new params: 1.99
----> avg loss exist params: 1.31
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.25
---->   avg loss new params: 0.41
----> avg loss exist params: 1.27
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 1.7s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.34
---->   avg loss new params: 0.54
----> avg loss exist params: 1.24
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

BATCH NUMBER:   6
PARAMS SAMPLED: 100

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.2s
---->   min loss new params: 1.06
---->   avg loss new params: 1.68
----> avg loss exist params: 1.26
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.2s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.53
---->   avg loss new params: 1.91
----> avg loss exist params: 1.28
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.38
---->   avg loss new params: 1.69
----> avg loss exist params: 1.3
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.3
---->   avg loss new params: 0.65
----> avg loss exist params: 1.28
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 1.6s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.2
---->   avg loss new params: 0.53
----> avg loss exist params: 1.25
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

BATCH NUMBER:   7
PARAMS SAMPLED: 120

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.7
---->   avg loss new params: 1.3
----> avg loss exist params: 1.25
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.36
---->   avg loss new params: 1.58
----> avg loss exist params: 1.26
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.51
---->   avg loss new params: 2.12
----> avg loss exist params: 1.29
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.43
---->   avg loss new params: 0.54
----> avg loss exist params: 1.27
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 1.9s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.29
---->   avg loss new params: 0.48
----> avg loss exist params: 1.24
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

BATCH NUMBER:   8
PARAMS SAMPLED: 140

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 1.54
---->   avg loss new params: 2.54
----> avg loss exist params: 1.28
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 1.05
---->   avg loss new params: 1.86
----> avg loss exist params: 1.3
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 1.0
---->   avg loss new params: 1.73
----> avg loss exist params: 1.31
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RandomForestSampler
Warning: Repeated samples still found after 5 duplication passes. This is probably due to a small search space.
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.54
---->   avg loss new params: 0.64
----> avg loss exist params: 1.29
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 1.7s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.21
---->   avg loss new params: 0.32
----> avg loss exist params: 1.27
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

BATCH NUMBER:   9
PARAMS SAMPLED: 160

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.81
---->   avg loss new params: 2.2
----> avg loss exist params: 1.29
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.88
---->   avg loss new params: 2.12
----> avg loss exist params: 1.31
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.63
---->   avg loss new params: 2.04
----> avg loss exist params: 1.33
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.31
---->   avg loss new params: 0.6
----> avg loss exist params: 1.31
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 1.4s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.2
---->   avg loss new params: 0.29
----> avg loss exist params: 1.29
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

BATCH NUMBER:   10
PARAMS SAMPLED: 180

METHOD: RandomUniformSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.67
---->   avg loss new params: 1.12
----> avg loss exist params: 1.28
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: RSequenceSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.58
---->   avg loss new params: 1.31
----> avg loss exist params: 1.28
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

METHOD: HaltonSampler
----> sim exec elapsed time: 0.2s
---->   min loss new params: 0.52
---->   avg loss new params: 1.26
----> avg loss exist params: 1.28
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.2s

METHOD: RandomForestSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.35
---->   avg loss new params: 0.44
----> avg loss exist params: 1.27
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 1.6s

METHOD: BestBatchSampler
----> sim exec elapsed time: 0.1s
---->   min loss new params: 0.31
---->   avg loss new params: 0.53
----> avg loss exist params: 1.25
---->         curr min loss: 0.09902649522714026
====>    total elapsed time: 0.1s

# best parameters obtained so far
params[0]

array([1.0002, 0.9791])

Plots

# index of mimumum loss
idxmin = np.argmin(cal.losses_samp)
param_min = cal.params_samp[idxmin]

# scatter plot of losses
plt.figure(figsize=(6, 5))

plt.scatter(
    cal.params_samp[:, 0], cal.params_samp[:, 1], c=cal.losses_samp, edgecolor="k"
)

plt.scatter(param_min[0], param_min[1], marker="x", s=500, color="y")
plt.scatter(true_params[0], true_params[1], marker="x", s=500, color="red")


plt.xlabel("mean")
plt.ylabel("variance")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fa43a70f890>

from black_it.utils.time_series import get_mom_ts

real_moments = get_mom_ts(cal.real_data)
simulated_moments = get_mom_ts(cal.series_samp[idxmin, 0, :, :])

# agreement of moments
plt.figure()
plt.plot(real_moments, "-o", label="real series")
plt.plot(simulated_moments, "-o", label="simulated series")
plt.xlabel("moment index")
plt.xticks(np.arange(18))
plt.ylabel("moments")
plt.legend()

<matplotlib.legend.Legend at 0x7fa43a589510>

# convergence

losses_per_batch = [
    cal.losses_samp[cal.batch_num_samp == i]
    for i in range(int(max(cal.batch_num_samp)) + 1)
]
mins_per_batch = np.array([np.min(l) for l in losses_per_batch])
cummin_per_batch = [
    np.min(mins_per_batch[: i + 1]) for i in range(mins_per_batch.shape[0])
]

plt.figure()
plt.plot(mins_per_batch, "-o", label="batch minimum")
plt.plot(cummin_per_batch, "-o", label="cumulative minimum")
plt.ylabel("loss")
plt.xlabel("batch number")
plt.legend()

<matplotlib.legend.Legend at 0x7fa4334355d0>