From 3f271e24b73e790de4abd7023da8c317fb3633ea Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 16 Jun 2024 18:01:59 +0300 Subject: [PATCH] doc: updated bayesian tutorial --- tutorials/bayesian.py | 151 ++++++++++++++++++++++++++++-------------- 1 file changed, 102 insertions(+), 49 deletions(-) diff --git a/tutorials/bayesian.py b/tutorials/bayesian.py index 42e16a6f..653d6b3d 100755 --- a/tutorials/bayesian.py +++ b/tutorials/bayesian.py @@ -21,19 +21,26 @@ Based on the above definition, we construct some prior models in the frequency domain, convert each of them to the time domain and use such an ensemble -to estimate the prior mean :math:`\mu_\mathbf{x}` and model -covariance :math:`\mathbf{C_x}`. +to estimate the prior mean :math:`\mathbf{x}_0` and model +covariance :math:`\mathbf{C}_{x_0}`. -We then create our data by sampling the true signal at certain locations +We then create our data by sampling the true signal at certain locations and solve the resconstruction problem within a Bayesian framework. Since we are assuming gaussianity in our priors, the equation to obtain the posterion mean -can be derived analytically: +and covariance can be derived analytically: .. math:: - \mathbf{x} = \mathbf{x_0} + \mathbf{C}_x \mathbf{R}^T - (\mathbf{R} \mathbf{C}_x \mathbf{R}^T + \mathbf{C}_y)^{-1} (\mathbf{y} - + \mathbf{x} = \mathbf{x_0} + \mathbf{C}_{x_0} \mathbf{R}^T + (\mathbf{R} \mathbf{C}_{x_0} \mathbf{R}^T + \mathbf{C}_y)^{-1} (\mathbf{y} - \mathbf{R} \mathbf{x_0}) +and + +.. math:: + \mathbf{C}_x = \mathbf{C}_{x_0} - \mathbf{C}_{x_0} \mathbf{R}^T + (\mathbf{R} \mathbf{C}_x \mathbf{R}^T + \mathbf{C}_y)^{-1} + \mathbf{R} \mathbf{C}_{x_0} + """ import matplotlib.pyplot as plt @@ -80,14 +87,15 @@ def prior_realization(f0, a0, phi0, sigmaf, sigmaa, sigmaphi, dt, nt, nfft): sigmaa = [0.1, 0.5, 0.6] phi0 = [-90.0, 0.0, 0.0] sigmaphi = [0.1, 0.2, 0.4] -sigmad = 1e-2 +sigmad = 1 +scaling = 100 # Scale by a factor to allow noise std=1 # Prior models nt = 200 nfft = 2**11 dt = 0.004 t = np.arange(nt) * dt -xs = np.array( +xs = scaling * np.array( [ prior_realization(f0, a0, phi0, sigmaf, sigmaa, sigmaphi, dt, nt, nfft) for _ in range(nreals) @@ -95,7 +103,10 @@ def prior_realization(f0, a0, phi0, sigmaf, sigmaa, sigmaphi, dt, nt, nfft): ) # True model (taken as one possible realization) -x = prior_realization(f0, a0, phi0, [0, 0, 0], [0, 0, 0], [0, 0, 0], dt, nt, nfft) +x = scaling * prior_realization( + f0, a0, phi0, [0, 0, 0], [0, 0, 0], [0, 0, 0], dt, nt, nfft +) + ############################################################################### # We have now a set of prior models in time domain. We can easily use sample @@ -110,7 +121,7 @@ def prior_realization(f0, a0, phi0, sigmaf, sigmaa, sigmaphi, dt, nt, nfft): N = 30 # lenght of decorrelation diags = np.array([Cm[i, i - N : i + N + 1] for i in range(N, nt - N)]) diag_ave = np.average(diags, axis=0) -# add a taper at the end to avoid edge effects +# add a taper at the start and end to avoid edge effects diag_ave *= np.hamming(2 * N + 1) fig, ax = plt.subplots(1, 1, figsize=(12, 4)) @@ -157,65 +168,107 @@ def prior_realization(f0, a0, phi0, sigmaf, sigmaa, sigmaphi, dt, nt, nfft): ynmask = Rop.mask(x + n) ############################################################################### -# First we apply the Bayesian inversion equation -xbayes = x0 + Cm_op * Rop.H * ( +# First, since the problem is rather small, we construct the dense version of +# all our matrices and we compute the analytical posterior mean and covariance + +Cm = Cm_op.todense() +Cd = Cd_op.todense() +R = Rop.todense() + +# Bayesian analytical solution +xpost_ana = x0 + Cm @ R.T @ (np.linalg.solve(R @ Cm @ R.T + Cd, yn - R @ x0)) +Cmpost_ana = Cm - Cm @ R.T @ (np.linalg.solve(R @ Cm @ R.T + Cd, R @ Cm)) + +############################################################################### +# Next we solve the same Bayesian inversion equation iteratively. We will see +# that provided we use enough iterations we can retrieve the same values of +# the analytical posterior mean +xpost_iter = x0 + Cm_op * Rop.H * ( lsqr(Rop * Cm_op * Rop.H + Cd_op, yn - Rop * x0, iter_lim=400)[0] ) -# Visualize +############################################################################### +# But what is the problem did not allow creating dense matrices for both the +# operator and the input covariance matrices. In this case, we can resort to the +# Randomize-Then-Optimize algorithm of Bardsley et al., 2014, which simply solves +# the same problem that we solved to find the MAP solution repeatedly by adding +# random noise to the data. It can be shown that the sample mean and covariance +# of the solutions of the different perturbed problems provide a good +# approximation for the true posterior mean and covariance. + +# RTO number of solutions +nreals = 1000 + +xrto = [] +for ireal in range(nreals): + yreal = yn + Rop * np.random.normal(0, sigmad, nt) + xrto.append( + x0 + + Cm_op + * Rop.H + * (lsqr(Rop * Cm_op * Rop.H + Cd_op, yreal - Rop * x0, iter_lim=400))[0] + ) + +xrto = np.array(xrto) +xpost_rto = np.average(xrto, axis=0) +Cmpost_rto = ((xrto - xpost_rto).T @ (xrto - xpost_rto)) / nreals + +############################################################################### +# Finally we visualize the different results + +# Means fig, ax = plt.subplots(1, 1, figsize=(12, 5)) ax.plot(t, x, "k", lw=6, label="true") +ax.plot(t, xpost_ana, "r", lw=7, label="bayesian inverse (ana)") +ax.plot(t, xpost_iter, "g", lw=5, label="bayesian inverse (iter)") +ax.plot(t, xpost_rto, "b", lw=3, label="bayesian inverse (rto)") ax.plot(t, ymask, ".k", ms=25, label="available samples") ax.plot(t, ynmask, ".r", ms=25, label="available noisy samples") -ax.plot(t, xbayes, "r", lw=3, label="bayesian inverse") ax.legend() -ax.set_title("Signal") +ax.set_title("Mean reconstruction") ax.set_xlim(0, 0.8) -plt.tight_layout() -############################################################################### -# So far we have been able to estimate our posterion mean. What about its -# uncertainties (i.e., posterion covariance)? -# -# In real-life applications it is very difficult (if not impossible) -# to directly compute the posterior covariance matrix. It is much more -# useful to create a set of models that sample the posterion probability. -# We can do that by solving our problem several times using different prior -# realizations as starting guesses: - -xpost = [ - x0 - + Cm_op - * Rop.H - * (lsqr(Rop * Cm_op * Rop.H + Cd_op, yn - Rop * x0, iter_lim=400)[0]) - for x0 in xs[:30] -] -xpost = np.array(xpost) - -x0post = np.average(xpost, axis=0) -Cm_post = ((xpost - x0post).T @ (xpost - x0post)) / nreals - -# Visualize +# RTO realizations fig, ax = plt.subplots(1, 1, figsize=(12, 5)) ax.plot(t, x, "k", lw=6, label="true") -ax.plot(t, xpost.T, "--r", lw=1) -ax.plot(t, x0post, "r", lw=3, label="bayesian inverse") +ax.plot(t, xrto[::10].T, "--b", lw=0.5) +ax.plot(t, xpost_rto, "b", lw=3, label="bayesian inverse (rto)") ax.plot(t, ymask, ".k", ms=25, label="available samples") ax.plot(t, ynmask, ".r", ms=25, label="available noisy samples") ax.legend() -ax.set_title("Signal") +ax.set_title("RTO realizations") ax.set_xlim(0, 0.8) -fig, ax = plt.subplots(1, 1, figsize=(5, 4)) -im = ax.imshow( - Cm_post, interpolation="nearest", cmap="seismic", extent=(t[0], t[-1], t[-1], t[0]) +# Covariances +fig, axs = plt.subplots(1, 2, figsize=(12, 4)) +axs[0].imshow( + Cmpost_ana, + interpolation="nearest", + cmap="seismic", + vmin=-5e-1, + vmax=2, + extent=(t[0], t[-1], t[-1], t[0]), +) +axs[0].set_title(r"$\mathbf{C}_m^{post,ANA}$") +axs[0].axis("tight") + +axs[1].imshow( + Cmpost_rto, + interpolation="nearest", + cmap="seismic", + vmin=-5e-1, + vmax=2, + extent=(t[0], t[-1], t[-1], t[0]), ) -ax.set_title(r"$\mathbf{C}_m^{posterior}$") -ax.axis("tight") +axs[1].set_title(r"$\mathbf{C}_m^{post,RTO}$") +axs[1].axis("tight") plt.tight_layout() ############################################################################### # Note that here we have been able to compute a sample posterior covariance -# from its estimated samples. By displaying it we can see how both the overall +# from its estimated samples. By displaying it we can see how both the overall # variances and the correlation between different parameters have become -# narrower compared to their prior counterparts. +# narrower compared to their prior counterparts. Moreover, whilst the RTO +# covariance seems to be slightly under-estimated, this represents an appealing +# alternative to the closed-form solution for large-scale problems under +# Gaussian assumptions.