Note
Click here to download the full example code
Loglikelihood of timeseries with Exponential correlation and vector noise.
This example shows how to use tripy.chol_loglike_1D and how it compares to
scipy.stats.multivariate_normal for different numbers of observations
under the assumption of Exponential correlation.
We assume that the measurements are obtained from the following probabilistic model:
\[\mathbf{X}_{\mathrm{model}}(\mathbf{\theta}) =
\mathbf{K}(\boldsymbol{\theta}) \cdot f(\mathbf{t})
+ \mathbf{E}_{\mathrm{meas}},\]
where:
\(\mathbf{X}_{\mathrm{model}}\) is the vector of model predictions.
\(f(\mathbf{t})\) is a vector valued function.
\(\mathbf{K}(\boldsymbol{\theta})\) is the correlated multiplicative
uncertainty factor.
\(\mathbf{E}_{\mathrm{meas}}\) i.i.d. Gaussian random variables.
\(\boldsymbol{\theta}\) is a set of parameters of the probabilistic model.
- The multiplicative factor \(\mathbf{K}\) and additive factor \(\mathbf{E}\)
are distributed as
:math:`mathbf{K}(boldsymbol{theta}) sim mathcal{N} left[ 1.0, boldsymbol{Sigma}_{mathrm{K}}(boldsymbol{theta})
right]`
- and :math:`mathbf{E}_{mathrm{meas}} sim mathcal{N}(0,
boldsymbol{Sigma}_{mathrm{E}}`,
respectively. This yields:
\[\mathbf{X}_{\mathrm{model}}(\boldsymbol{\theta}) \sim \mathcal{N}
\left[ f(\mathbf{t}), \mathbf{\Sigma}(\boldsymbol{\theta}) \right],\]
with \(\boldsymbol{\Sigma}(\boldsymbol{\theta}) = f(\mathbf{t}) \cdot \boldsymbol{\Sigma}_{\mathrm{K}} + \boldsymbol{\Sigma}_{\mathrm{E}}\)
Import packages:
from timeit import default_timer as timer
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
from tripy.kernels import Exponential
from tripy.loglikelihood import chol_loglike_1D
Problem setup
N = 100 # Number of points
std_meas = np.repeat(
2.0, N
) # Standard deviation of the additive measurement uncertainty
std_model = np.repeat(0.5, N) # Standard deviation of the modeling uncertainty
l_corr = 10.0 # Correlation lengthscale
# Coordinate vector
t = np.sort(np.linspace(0, 10, N) + 0.1 * np.random.rand(N))
# Vector valued function evaluation
y_func = np.cos(t) + 5
# Observations
v_obs = np.random.rand(N) + y_func
Reference solution using scipy:
e_cov_mx = np.diag(std_meas**2)
kernel = Exponential(np.reshape(t, (-1, 1)))
corr_mx = kernel.eval(std_model, length_scale=l_corr)
kph_cov_mx = np.matmul(np.diag(y_func), np.matmul(corr_mx, np.diag(y_func)))
cov_mx = kph_cov_mx + e_cov_mx
logL_ref = stats.multivariate_normal.logpdf(v_obs, cov=cov_mx)
Timing comparison between conventional and efficient solutions:
N_iter = 10
# Reference
N_vec = np.arange(100, 1000, 100)
t_ref_list = []
for i, N in enumerate(N_vec):
# Vectors of measurement and modeling uncertainty std. dev.
std_meas = np.repeat(2.0, N)
std_model = np.repeat(0.5, N)
# Coord vector and a test function
t = np.sort(np.linspace(0, 10, N) + 0.1 * np.random.rand(N))
y_func = np.cos(t) + 5
# Observations
v_obs = np.random.rand(N) + y_func
e_cov_mx = np.diag(std_meas**2)
kernel = Exponential(np.reshape(t, (-1, 1)))
# Reference
t1 = timer()
for _j in range(N_iter):
corr_mx = kernel.eval(std_model, length_scale=l_corr)
kph_cov_mx = np.matmul(np.diag(y_func), np.matmul(corr_mx, np.diag(y_func)))
cov_mx = kph_cov_mx + e_cov_mx
test = stats.multivariate_normal.logpdf(v_obs, cov=cov_mx)
t2 = timer()
t_ref_list.append((t2 - t1) / N_iter)
print("=============================")
print(f"Iter {i + 1} / {len(N_vec)}")
print(f"N = {N}")
print(f"t = {t_ref_list[-1]}")
print("=============================")
# Efficient solution
p_vec = np.linspace(2, 6, 50)
t_list = []
for i, p in enumerate(p_vec):
N = int(10**p)
# Vectors of measurement and modeling uncertainty std. dev.
std_meas = np.repeat(2.0, N)
std_model = np.repeat(0.5, N)
# Coord vector and a test function
t = np.sort(np.linspace(0, 10, N) + 0.1 * np.random.rand(N))
y_func = np.cos(t) + 5
# Observations
v_obs = np.random.rand(N) + y_func
# Initial evaluation to perform jit compilation using numba.
_ = chol_loglike_1D(v_obs, t, l_corr, std_model, std_meas=std_meas, y_model=y_func)
# Call function
t1 = timer()
for _j in range(N_iter):
test = chol_loglike_1D(
v_obs, t, l_corr, std_model, std_meas=std_meas, y_model=y_func
)
t2 = timer()
t_list.append((t2 - t1) / N_iter)
print("=============================")
print(f"Iter {i + 1} / {len(p_vec)}")
print(f"p = {p}")
print(f"t = {t_list[-1]}")
print("=============================")
=============================
Iter 1 / 9
N = 100
t = 0.0036246815000595236
=============================
=============================
Iter 2 / 9
N = 200
t = 0.02045987489991603
=============================
=============================
Iter 3 / 9
N = 300
t = 0.044014394900023034
=============================
=============================
Iter 4 / 9
N = 400
t = 0.09252052540005025
=============================
=============================
Iter 5 / 9
N = 500
t = 0.12648369179996735
=============================
=============================
Iter 6 / 9
N = 600
t = 0.21510887879994697
=============================
=============================
Iter 7 / 9
N = 700
t = 0.27731878619997585
=============================
=============================
Iter 8 / 9
N = 800
t = 0.34173714859998655
=============================
=============================
Iter 9 / 9
N = 900
t = 0.44427265500007707
=============================
=============================
Iter 1 / 50
p = 2.0
t = 0.00012304689998927644
=============================
=============================
Iter 2 / 50
p = 2.0816326530612246
t = 0.00011380969999663648
=============================
=============================
Iter 3 / 50
p = 2.163265306122449
t = 0.00012194929995530401
=============================
=============================
Iter 4 / 50
p = 2.2448979591836733
t = 0.00012368229999992764
=============================
=============================
Iter 5 / 50
p = 2.326530612244898
t = 0.0001289259000259335
=============================
=============================
Iter 6 / 50
p = 2.4081632653061225
t = 0.00012926529998367187
=============================
=============================
Iter 7 / 50
p = 2.489795918367347
t = 0.0001337468000201625
=============================
=============================
Iter 8 / 50
p = 2.571428571428571
t = 0.00013729530001000968
=============================
=============================
Iter 9 / 50
p = 2.6530612244897958
t = 0.00013877950004825834
=============================
=============================
Iter 10 / 50
p = 2.7346938775510203
t = 0.00015546220001851906
=============================
=============================
Iter 11 / 50
p = 2.816326530612245
t = 0.000162010999974882
=============================
=============================
Iter 12 / 50
p = 2.8979591836734695
t = 0.00016884699998627185
=============================
=============================
Iter 13 / 50
p = 2.979591836734694
t = 0.00017914220006787218
=============================
=============================
Iter 14 / 50
p = 3.061224489795918
t = 0.00018904159996964155
=============================
=============================
Iter 15 / 50
p = 3.142857142857143
t = 0.0002114910999807762
=============================
=============================
Iter 16 / 50
p = 3.224489795918367
t = 0.00022334539999064874
=============================
=============================
Iter 17 / 50
p = 3.3061224489795915
t = 0.00024360530005651526
=============================
=============================
Iter 18 / 50
p = 3.387755102040816
t = 0.0002743873999861535
=============================
=============================
Iter 19 / 50
p = 3.4693877551020407
t = 0.00030056140003580366
=============================
=============================
Iter 20 / 50
p = 3.5510204081632653
t = 0.0003395527999600745
=============================
=============================
Iter 21 / 50
p = 3.63265306122449
t = 0.00038942380006119494
=============================
=============================
Iter 22 / 50
p = 3.7142857142857144
t = 0.0004407431999425171
=============================
=============================
Iter 23 / 50
p = 3.7959183673469385
t = 0.0005189471000448976
=============================
=============================
Iter 24 / 50
p = 3.877551020408163
t = 0.0005905097000322712
=============================
=============================
Iter 25 / 50
p = 3.9591836734693877
t = 0.0006981109999287582
=============================
=============================
Iter 26 / 50
p = 4.040816326530612
t = 0.0008135963000313496
=============================
=============================
Iter 27 / 50
p = 4.122448979591836
t = 0.0009729644000799453
=============================
=============================
Iter 28 / 50
p = 4.204081632653061
t = 0.0011343375000251398
=============================
=============================
Iter 29 / 50
p = 4.285714285714286
t = 0.0013835774000654055
=============================
=============================
Iter 30 / 50
p = 4.36734693877551
t = 0.001651009499983047
=============================
=============================
Iter 31 / 50
p = 4.448979591836734
t = 0.0019680749999679394
=============================
=============================
Iter 32 / 50
p = 4.530612244897959
t = 0.0023962232000485527
=============================
=============================
Iter 33 / 50
p = 4.612244897959183
t = 0.002935126099964691
=============================
=============================
Iter 34 / 50
p = 4.6938775510204085
t = 0.003544928900009836
=============================
=============================
Iter 35 / 50
p = 4.775510204081632
t = 0.004340908700032742
=============================
=============================
Iter 36 / 50
p = 4.857142857142857
t = 0.005365556199922139
=============================
=============================
Iter 37 / 50
p = 4.938775510204081
t = 0.0064890468000157854
=============================
=============================
Iter 38 / 50
p = 5.020408163265306
t = 0.007919170000059239
=============================
=============================
Iter 39 / 50
p = 5.1020408163265305
t = 0.016144610100036517
=============================
=============================
Iter 40 / 50
p = 5.183673469387754
t = 0.01915125449995685
=============================
=============================
Iter 41 / 50
p = 5.26530612244898
t = 0.023485571399942274
=============================
=============================
Iter 42 / 50
p = 5.346938775510203
t = 0.028397505200064187
=============================
=============================
Iter 43 / 50
p = 5.428571428571429
t = 0.03439178850003373
=============================
=============================
Iter 44 / 50
p = 5.5102040816326525
t = 0.04267896829996971
=============================
=============================
Iter 45 / 50
p = 5.591836734693877
t = 0.05176645910005391
=============================
=============================
Iter 46 / 50
p = 5.673469387755102
t = 0.06384881459998723
=============================
=============================
Iter 47 / 50
p = 5.755102040816326
t = 0.06806369309997536
=============================
=============================
Iter 48 / 50
p = 5.836734693877551
t = 0.0820311565999873
=============================
=============================
Iter 49 / 50
p = 5.918367346938775
t = 0.11317631760002769
=============================
=============================
Iter 50 / 50
p = 6.0
t = 0.1331412384000032
=============================
Plot a comparison of the wall clock times:
fig = plt.figure()
plt.plot([int(10**p) for p in p_vec], t_list, label="Efficient evaluation")
plt.plot(N_vec, t_ref_list, label="Naive evaluation")
plt.xlabel("No. of points")
plt.ylabel("Wall clock time [s]")
plt.xscale("log")
plt.legend(bbox_to_anchor=(0.8, 0.85), bbox_transform=fig.transFigure)
plt.grid()
plt.show()
Total running time of the script: ( 0 minutes 26.065 seconds)