Model Fitting E.g. 1: Cosine
Contents
import mmf_setup
mmf_setup.nbinit()
import logging
logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
This cell adds /home/docs/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/checkouts/fall-2021 to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:
- Choose "Trust Notebook" from the "File" menu.
- Re-execute this cell.
- Reload the notebook.
Model Fitting E.g. 1: Cosine¶
Here we go through the step-by-step procedure for fitting data to the following model:
We shall assume that the errors are independently distributed, considering both the case of gaussian errors, and non-gaussian:
We start with the gaussian case, and assume that the errors are sufficiently small that the standard analysis techniques based on \(\chi^2(\vect{a})\) apply as discussed in chapter 15 of [Press et al., 2007]. We then test these using Monte Carlo (MC) and then use an MCMC approach to solve the general case.
Linear Gaussian Approximation¶
We start with an “exact” solution, then generate some data to analyze at \(N_t\) equally spaced points for \(t \in [0, t_\max]\).
from IPython.display import Latex
from collections import namedtuple
from uncertainties import correlated_values, unumpy as unp
from scipy.optimize import least_squares
from scipy.stats import chi2
Nt = 12
#Nt = 100
t_max = 10.0
# Exact parameter values
Params = namedtuple("Params", ["w", "c", "A", "phi"])
a_exact = Params(w=2 * np.pi / 5, c=2.1, A=3.4, phi=5.6)
w, c, A, phi = a_exact
def f(t, w, c, A, phi, np=np):
"""Model function."""
return c + A * np.cos(w * t + phi)
# Exact data and experimental errors
t = np.linspace(0, t_max, Nt)
y = f(t, *a_exact)
sigmas = 0.5 * np.ones(Nt)
# Randomly generated data... An "Experiment"
rng = np.random.default_rng(seed=2)
ydata = rng.normal(loc=y, scale=sigmas)
ts = np.linspace(0, t_max) # Many points for a smooth curve.
fig, ax = plt.subplots(figsize=(10, 5))
ax.errorbar(t, ydata, yerr=sigmas, fmt="C0.", ecolor="C1", label="data")
ax.plot(ts, f(ts, *a_exact), "-C2", label="exact")
ax.set(xlabel="$t$", ylabel="$f(t)$")
ax.legend();
Chi Square Fit¶
Now we do the fit. We use scipy.optimize.least_squares() which requires a
function that returns the list of weighted residuals:
This minimizes the following cost function, for which the Hessian is the inverse covariance matrix:
def get_residuals(p, xdata=t, ydata=ydata, sigmas=sigmas):
"""Return residuals for least_squares."""
return (ydata - f(xdata, *p))/sigmas
# By using simple numpy functions, we can compute the jacobian
# with the complex-step method.
kw = dict(jac='cs')
res = least_squares(fun=get_residuals, x0=a_exact, **kw)
p = Params(*res.x)
C = np.linalg.inv(res.jac.T @ res.jac)
r = get_residuals(p)
nu = len(r) - len(p)
chi2_r = np.sum(abs(r)**2) / nu
Q = 1 - chi2.cdf(chi2_r*nu, df=nu)
Latex(rf"$\chi^2_r = {chi2_r:.2g}, \qquad Q = {Q:.2g}$")
Here we have used the Chi-Squared Distribution \(P_{\nu,\chi^2}(\chi^2)\) to calculate the chi-square probability \(Q\), sometimes called the tail distribution:
This is the probability that \(\chi^2\) exceeds the value found, and is the complement of the cumulative distribution function (CDF). This is sometimes called a one-sided \(p\)-value and is used as a test statistic to assess whether or not the model is reasonable. As discussed in section 15.1 of [Press et al., 2007], values of \(Q> 0.001\) are not too bad:
Truly wrong models have values will often be rejected with vastly smaller values of \(Q\), \(10^{-18}\), say.
Here we plot the distribution of the reduced chi squared for our dataset with \(\nu = N - M\) degrees of freedom:
from myst_nb import glue
# Generate actual distribution of chi2_r using MC
from functools import partial
def get_chi2_r(ydata=None, random=rng.normal):
"""Return chi2_r from a sample experiment and fit."""
global y, sigmas, a_exact
if ydata is None:
ydata = random(loc=y, scale=sigmas)
fun = partial(get_residuals, ydata=ydata)
res = least_squares(fun=fun, x0=a_exact, jac='cs')
p = Params(*res.x)
chi2_r = sum(abs(fun(p))**2) / nu
return chi2_r
Ns = 2000 # number of samples
chi2_rs = [get_chi2_r(random=rng.normal) for n in range(Ns)]
_chi2_r = np.linspace(0, 3, 500)
_i = np.where(_chi2_r >= chi2_r)[0][0]
fig, ax = plt.subplots()
ax.plot(_chi2_r, nu*chi2.pdf(nu*_chi2_r, df=nu), "-C0", label=rf"PDF ($\nu={nu}$)")
ax.plot(_chi2_r, nu*chi2.pdf(nu*_chi2_r, df=nu-1), ":C0", label=rf"$\nu={nu}\pm 1$")
ax.plot(_chi2_r, nu*chi2.pdf(nu*_chi2_r, df=nu+1), ":C0")
ax.plot(_chi2_r, chi2.cdf(nu*_chi2_r, df=nu), "--C1", label="CDF")
kw = dict(bins=50, histtype='step', alpha=0.8, density=True)
ax.hist(chi2_rs, ec="C2", label=f"{Ns} samples", **kw)
ax.fill_between(_chi2_r[_i:], nu*chi2.pdf(nu*_chi2_r[_i:], df=nu), fc="C0")
ax.axvline([chi2_r], ls=":", c="y")
ax.axhline([1-Q], ls=":", c="y")
ax.set(xlim=(0, 4),
xlabel=r"$\chi^2_r=\chi^2/\nu$",
title=fr"Distribution of $\chi^2_r$ with $\nu={nu} = {Nt}-{len(p)}$ degrees of freedom")
ax.legend();
Parameter Estimates and Uncertainties¶
The results of the fit are characterized by the best-fit parameter values \(\bar{\vect{a}}\) and the covariance matrix \(\mat{C}\) such that, to lowest order, \(\chi^2(\vect{a})\) is quadratic:
The confidence region corresponding to confidence level \(p\) is given by the region within the contour
where the contour \(\chi^2_{p}\) is carefully chosen so that this region contains fraction \(p\) of the samples. This can bed computed from the inverse of the cumulative distribution function (CDF) for the chi squared distribution with \(\nu\) degrees of freedom:
if:
the fit is good,
the posterior distribution is well approximated by a multivariate normal distribution (i.e. if the experimental errors are gaussian and the model is linear in the region of these errors), and
you use the correct value of \(\nu\) corresponding to the number of free parameters you wish to consider. If you want to consider a marginal distribution of \(\tilde{\nu} < \nu\) parameters, then just extract the corresponding \(\tilde{nu}\) rows and columns from \(\mat{C}\) (not \(\mat{C}^{-1}\)).
In our case, the errors are gaussian as we saw above, so we can proceed with the simplified approach. What we do here is use the uncertainties package to perform a forward error propagation of the parameter covariance matrix \(\mat{C}\) through to the actual function values for the model \(y = f(x, \vect{a})\). We then treat each \(y\) as a single parameter (\(\tilde{\nu}=1\)) and use the \(n\sigma\) values where \(\sigma\) is the standard deviation computed by the uncertainties package to produce a \(1\sigma\) and \(2\sigma\) uncertainty band corresponding to our parameter values:
# Here we use the uncertainties package to do linear
# error propagation of are parameter covariances
p_ = Params(*correlated_values(p, C))
pr_ = Params(*correlated_values(p, C*chi2_r))
def plot_band(t, y_, sigma=1, ax=None, **kw):
"""Plot a band using correlated errors in y_."""
### To Do: These bands are valid in the Gaussian case, but should
### be properly scaled using the proper confidence analysis.
if ax is None:
ax = plt.gca()
y = unp.nominal_values(y_)
dy = unp.std_devs(y_)
ax.fill_between(t, y-dy*sigma, y+dy*sigma, **kw)
return ax
ts = np.linspace(0, t_max) # Many points for a smooth curve.
fig, axs = plt.subplots(2, 1, figsize=(10, 10))
for _p_, ax in zip((p_, pr_), axs):
ax.errorbar(t, ydata, yerr=sigmas, fmt="C0.", ecolor="C1", label="data")
ax.plot(ts, f(ts, *a_exact), "-C2", label="exact")
ax.plot(ts, f(ts, *p), "--C3", label="best fit")
plot_band(t=ts, y_=f(ts, *_p_, np=unp), ax=ax, fc="C3", alpha=0.5,
sigma=1, label=r"$1\sigma$ band")
plot_band(t=ts, y_=f(ts, *_p_, np=unp), ax=ax, fc="C3", alpha=0.2,
sigma=2, label=r"$2\sigma$ band")
ax.set(xlabel="$t$", ylabel="$f(t)$")
ax.legend();
axs[0].set(title=fr"$\chi^2_r={chi2_r:.2f}$, $Q={Q:.2g}$");
axs[1].set(title=r"Rescaled covariance so $\chi^2_r=1$");
In the lower plot, we have rescaled the covariance matrix \(\mat{C} \rightarrow \chi^2_r \mat{C}\) as if we did not trust our error estimates, and interpreted the large \(\chi^2_r\) value as a sign that we underestimated the errors.
We can now look at some of the marginal distributions. One can use a set of corner plots to summarize this information, but we first do this explicitly.
Single Parameter Distributions¶
We start with the individual parameter marginal distributions. In this case \(\tilde{\nu} = 1\). Since our posterior is well approximated by a gaussian, we can simply look at the standard deviations given as the square root of the corresponding diagonal entry of the covariance matrix: \(\sigma_i = \sqrt{C_{ii}}\). This is nicely summarized by printing the parameters using the uncertainties package as we show below:
from scipy.stats import norm
fig, ax = plt.subplots()
_x = np.linspace(0, 7, 200)
for _name, _p, _Cii in zip(Params._fields, p_, np.diag(C)):
print(f"{_name}: {_p:.2uS}, (sqrt(C_ii) = {np.sqrt(_Cii):.2g})")
ax.plot(_x, norm.pdf(_x, loc=_p.n, scale=_p.s), label=_name)
ax.legend();
w: 1.237(21), (sqrt(C_ii) = 0.021)
c: 2.14(15), (sqrt(C_ii) = 0.15)
A: 3.47(20), (sqrt(C_ii) = 0.2)
phi: 5.79(12), (sqrt(C_ii) = 0.12)
Pairwise Distributions¶
Supposed we are most interested in the frequency \(\omega\) and phase \(\phi\). We can consider the marginal distribution of these two parameters by extracting the appropriate columns and rows. Note that we must now choose our contours using \(\tilde{\nu} = 2\). We display these as contours of \(\chi^2\):
from scipy.stats import norm, chi2
i, j = map(Params._fields.index, ['w', 'phi'])
w_, phi_ = p_.w, p_.phi
# Extract rows and colums. Note: C[[i, j], [i, j]] does not work.
Cij = C[[i, j], :][:, [i, j]]
# We can also do this with the uncertainties package
from uncertainties import covariance_matrix
assert np.allclose(Cij, covariance_matrix([w_, phi_]))
_sigma = 5
ws, phis = np.meshgrid(
np.linspace(w_.n-_sigma*w_.s, w_.n+_sigma*w_.s, 200),
np.linspace(phi_.n-_sigma*phi_.s, phi_.n+_sigma*phi_.s, 200))
# Deviation matrix to calculate chi2
dws_phis = np.array([ws-w_.n, phis-phi_.n])
# This is just the matrix product, but we want to do this
# over a bunch of parameter values, so we use einsum.
dchi2 = np.einsum('ij,ixy,jxy->xy',
np.linalg.inv(Cij), dws_phis, dws_phis)
# Use a normal distribution to compute the confidence levels
sigmas_ = np.array([1, 2, 3, 4])
ps = norm.cdf(sigmas_) - norm.cdf(-sigmas_)
# Now use chi2 distribution to get the corresponding contours
levels = chi2.ppf(ps, df=2)
fig, ax = plt.subplots(figsize=(10, 6))
_cs = ax.contour(ws, phis, dchi2, levels=levels, cmap="winter")
fmt = dict([(_l, fr"${_n}\sigma$")
for _l, _n, _p in zip(levels, sigmas_, ps)])
# Draw single-parameter confidence limits to show that
# they are smaller.
for _n, _sigma in enumerate(sigmas_):
kw = dict(c=_cs.collections[_n].get_edgecolor(),
alpha=0.5, zorder=-100)
for _s in [1, -1]:
ax.axhline(phi_.n + _s*_sigma*phi_.s, ls=":", **kw)
ax.axvline(w_.n + _s*_sigma*w_.s, ls="--", **kw)
ax.clabel(_cs, _cs.levels, inline=True, fmt=fmt, fontsize=10)
ax.set(xlabel="$\omega$", ylabel="$\phi$");
We show the \(1\sigma\), \(2\sigma\), \(3\sigma\), and \(4\sigma\) \(\tilde{\nu} = 2\)-parameter confidence regions demonstrating the correlations between \(\omega\) and \(\phi\). For comparison, we show the corresponding 1-parameter regions as dotted lines. E.g. for the \(1\sigma\) confidence regions, the same amount of posterior probability (68.27%) lies within the center blue ellipse as lies between the central vertical blue dashed lines, or between the central horizontal blue dotted lines.
In other words, the single-parameter confidence regions include points outside of the 2-parameter ellipse, hence must be smaller as shown to keep the same probability.
Here is the full corner plot using phys_581_2021.plotting.corner_plot():
from phys_581_2021.plotting import corner_plot
fig, axs = corner_plot(
p_, labels=[r"$\omega$", r"$c$", r"$A$", r"$\phi$"]);
Change of Variables¶
To compute a change of variables, we need to keep track of the Jacobian of the transformation so we can properly transform the covariance matrix. The uncertainties package does this automatically, computing the derivatives with automatic differentiation. This is only valid if the errors are small enough that the model is approximately linear, but is very convenient. Here we replace \(\phi\) with the location of the peak \(t_0\):
# make sure t0_ is positive and small by making phi0_ between -2pi and 0
phi0_ = phi_ % (2*np.pi) - 2*np.pi
t0_ = - phi0_ / w_
Cij = covariance_matrix([w_, t0_])
_sigma = 5
ws, t0s = np.meshgrid(
np.linspace(w_.n-_sigma*w_.s, w_.n+_sigma*w_.s, 200),
np.linspace(t0_.n-_sigma*t0_.s, t0_.n+_sigma*t0_.s, 200))
# Deviation matrix to calculate chi2
dws_t0s = np.array([ws-w_.n, t0s-t0_.n])
dchi2 = np.einsum('ij,ixy,jxy->xy',
np.linalg.inv(Cij), dws_t0s, dws_t0s)
# Now use chi2 distribution to get the corresponding contours
fig, ax = plt.subplots(figsize=(10, 6))
_cs = ax.contour(ws, t0s, dchi2, levels=levels, cmap="winter")
fmt = dict([(_l, fr"${_n}\sigma$")
for _l, _n, _p in zip(levels, sigmas_, ps)])
# Draw single-parameter confidence limits to show that
# they are smaller.
for _n, _sigma in enumerate(sigmas_):
kw = dict(c=_cs.collections[_n].get_edgecolor(),
alpha=0.5, zorder=-100)
for _s in [1, -1]:
ax.axhline(t0_.n + _s*_sigma*t0_.s, ls=":", **kw)
ax.axvline(w_.n + _s*_sigma*w_.s, ls="--", **kw)
ax.clabel(_cs, _cs.levels, inline=True, fmt=fmt, fontsize=10)
ax.set(xlabel="$\omega$", ylabel="$t_0$");
Non-Gaussian Errors¶
We now repeat the analysis, but with non-gaussian errors. Instead, we use the Gumbel distribution:
If we consider the normalized errors \(\tilde{e}_n\), they are distributed as:
The log-likelihood will be a sum of
We can use this to generate the various simulated \(P_{\nu}(\chi^2)\) distributions for Gumbel distributed errors.
en = np.linspace(-3, 10, 100)
plt.plot(en, en + np.exp(-en))
[<matplotlib.lines.Line2D at 0x7fa87ce28670>]
from scipy.interpolate import InterpolatedUnivariateSpline, UnivariateSpline
def get_cdf_ppf(samples=None):
"""Return interpolated cdf and ppf functions.
Results
-------
cdf : function
Cumulative distribution function `q=cdf(chi2_r)`.
ppf : function
Percent point function (inverse of `cdf`).
"""
Ns = len(samples)
ps = np.linspace(0, 1, Ns)
x = sorted(samples)
cdf = InterpolatedUnivariateSpline(x, ps, k=1, ext='const')
ppf = InterpolatedUnivariateSpline(ps, x, k=1, ext='const')
return cdf, ppf
def get_chi2_cdf_ppf(nu=1, random=rng.normal, Ns=10000):
"""Return `(cdf, pdf, samples)` for the chi2 distribution.
Arguments
---------
nu : int
Degrees of freedom.
random : function
Random number generator for errors.
Ns : int
Number of samples.
"""
en = random(size=(Ns, nu))
chi2 = (en**2).sum(axis=-1)
cdf, ppf = get_cdf_ppf(chi2)
return cdf, ppf, chi2
import scipy.stats
sp = scipy
chi2r_ = np.linspace(0, 4, 100)
hist_kw = dict(bins=np.linspace(0, 4, 100), histtype='step',
alpha=0.8, density=True)
fig, axs = plt.subplots(1, 2, figsize=(10, 4))
for ax, label, random in zip(axs,
["normal", "gumbel"],
[rng.normal, rng.gumbel]):
for _n, nu in enumerate([1, 2, 3, 10, 50]):
cdf, ppf, chi2r_samples = get_chi2_cdf_ppf(nu=nu, random=random)
chi2r_samples /= nu
c = f"C{_n}"
l, = ax.plot(chi2r_, cdf(nu*chi2r_),
"--", lw=1, c=c, label=fr"CDF $\nu={nu}$")
if label == "normal":
ax.plot(chi2r_, nu*sp.stats.chi2.pdf(nu*chi2r_, df=nu),
ls=":", c=c)
ax.plot(chi2r_, sp.stats.chi2.cdf(nu*chi2r_, df=nu),
ls=":", c=c)
ax.set(xlim=(0, 2))
else:
ax.set(xlim=(0, 4))
ax.hist(chi2r_samples, ec=c, **hist_kw)
ax.set(title=rf"$\chi^2_r$ distribution for {label} errors")
ax.set(ylim=(0, 2), xlabel=r"$\chi^2_r$")
ax.legend(loc='upper right')
plt.tight_layout()
On the left we show the CDF and PDF histogram for \(\chi^2_r\) generated from a sample of random numbers for a standard normal distribution, comparing with the analytic forms (dotted lines). On the right we use the same method to compute the CDF for the standard Gumbel distribution.
We now generate experimental data and proceed with an analysis:
# Randomly generated data...
# An "Experiment" with non-gaussian errors.
rng = np.random.default_rng(seed=2)
ydata = rng.gumbel(loc=y, scale=sigmas)
fig, ax = plt.subplots(figsize=(10, 5))
ax.errorbar(t, ydata, yerr=sigmas, fmt="C0.", ecolor="C1", label="data")
ax.plot(ts, f(ts, *a_exact), "-C2", label="exact")
ax.set(xlabel="$t$", ylabel="$f(t)$")
ax.legend();
Chi Square Fit¶
res = least_squares(fun=partial(get_residuals, ydata=ydata),
x0=a_exact, jac='cs')
p = Params(*res.x)
C = np.linalg.inv(res.jac.T @ res.jac)
r = get_residuals(p)
nu = len(r) - len(p)
chi2_r = np.sum(abs(r)**2) / nu
Q_wrong = 1 - sp.stats.chi2.cdf(chi2_r*nu, df=nu)
cdf_gumbel, pdf_gumbel, chi2s_gumbel = get_chi2_cdf_ppf(nu=nu, random=rng.gumbel)
Q = 1 - cdf_gumbel(chi2_r*nu)
Latex(r", \qquad ".join([
rf"$\chi^2_r = {chi2_r:.2g}",
rf"Q = {Q:.2g}",
rf"Q_{{\mathrm{{wrong}}}} = {Q_wrong:.2g}$"]))
With the non-gaussian errors, estimating the \(Q\) value with the chi squared distribution gives a very wrong interpretation about the quality of the fit. Instead, we must use the corresponding CDF for our Gumbel-distributed errors. To obtain the confidence region with confidence level \(p\), we must find the value of \(\chi^2_p\) where:
This is just the inverse CDF. We can check our method (slowly) by actually fitting the data repeatedly to generate the distribution of \(\chi^2\):
# Generate actual distribution of chi2_r using MC
from functools import partial
Ns = 2000 # number of samples
chi2s_data = np.array([nu*get_chi2_r(random=rng.gumbel) for n in range(Ns)])
chi2s_normal = np.array([nu*get_chi2_r(random=rng.normal) for n in range(Ns)])
cdf_data, ppf_data = get_cdf_ppf(chi2s_data)
_chi2_r = np.linspace(0, np.max(chi2s_gumbel)/nu, 500)
fig, ax = plt.subplots()
ax.plot(_chi2_r, cdf_data(nu*_chi2_r), '-C0', label="data")
ax.plot(_chi2_r, cdf_gumbel(nu*_chi2_r), ':C1', label="gumbel")
ax.plot(_chi2_r, sp.stats.chi2.cdf(nu*_chi2_r, df=nu), '--C2', label="normal")
ax.set(xlabel="$\chi^2_r$", ylabel="CDF",
title=fr"CDF for $P(\chi^2_r)$ for $\nu={nu}$ Gumbel-distributed errors.")
kw = dict(bins=100, histtype='step', alpha=0.8, density=True)
ax.hist(chi2s_data/nu, ec="C0", **kw)
ax.hist(chi2s_gumbel/nu, ec="C1", **kw)
ax.legend();
fig, ax = plt.subplots()
random = rng.gumbel
for Ns in [500, 1000]:
chi2s_data = np.array([nu*get_chi2_r(random=random) for n in range(Ns)])
cdf_data, ppf_data = get_cdf_ppf(chi2s_data)
cdf_gumbel, pdf_gumbel, chi2s_gumbel = get_chi2_cdf_ppf(nu=nu, random=random)
ax.plot(_chi2_r, cdf_data(nu*_chi2_r), '-', label=f"data {Ns}")
ax.plot(_chi2_r, cdf_gumbel(nu*_chi2_r), ':', label="gumbel")
ax.legend()
<matplotlib.legend.Legend at 0x7fa876d26fa0>
_chi2_r = np.linspace(0, 5, 500)
_i = np.where(_chi2_r >= chi2_r)[0][0]
fig, ax = plt.subplots()
ax.plot(_chi2_r, nu*sp.stats.chi2.pdf(nu*_chi2_r, df=nu), "-C1", label=rf"$\nu={nu}$")
ax.plot(_chi2_r, nu*sp.stats.chi2.pdf(nu*_chi2_r, df=nu-1), ":C1", label=rf"$\nu={nu}\pm 1$")
ax.plot(_chi2_r, nu*sp.stats.chi2.pdf(nu*_chi2_r, df=nu+1), ":C1")
ax.plot(_chi2_r, cdf_gumbel(nu*_chi2_r), "--C0", label=rf"CDF (Gumbel)")
ax.plot(_chi2_r, sp.stats.chi2.cdf(nu*_chi2_r, df=nu), "--C1", label=rf"CDF $\nu={nu}$")
# Make sure bins end on chi2_r
bins = np.linspace(0, chi2_r, 20)
_d = np.diff(bins).mean()
bins = np.concatenate([bins, np.arange(chi2_r+_d, 4+_d, _d)])
kw = dict(bins=bins, histtype='step', alpha=0.8, density=True)
res = ax.hist(chi2s_gumbel/nu, ec="C0", label=f"{Ns} samples (Gumbel)", **kw)
# Add filled region
xy = res[2][0].xy[res[2][0].xy[:, 0] >= chi2_r, :]
xy[0, 1] = 0
ax.add_patch(plt.Polygon(xy, fc="C0"))
ax.hist(chi2s_normal/nu, ls=":", ec="C1", label=fr"{Ns} samples ($\chi^2_{{\nu=3}}$)", **kw)
ax.axvline([chi2_r], ls=":", c="y")
ax.axhline([1-Q], ls=":", c="y")
ax.set(xlim=(0, 5),
xlabel=r"$\chi^2_r=\chi^2/\nu$",
title=fr"Distribution of $\chi^2_r$ for $\nu={nu}$ gumbel errors; Q={Q:.2g}")
ax.legend();
References¶
Samples, samples, everywhere…: A nice overview of different MCMC software accessible with python.