GetDist Jupyter Notebook Plot Gallery

Demonstrates the types of plot you can make with GetDist and how to make them. You can also run this notebook online.

# Show plots inline, and load main getdist plot module and samples class
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import os
import sys

sys.path.insert(0, os.path.realpath(os.path.join(os.getcwd(), "..")))
import getdist
import IPython
import matplotlib
import matplotlib.pyplot as plt
from getdist import MCSamples, plots

print("GetDist Version: %s, Matplotlib version: %s" % (getdist.__version__, matplotlib.__version__))

GetDist Version: 1.7.2, Matplotlib version: 3.10.3

# use this here *after* the above (instead of matplotlib inline) to use interactive plots
# %matplotlib notebook

# Get some random samples for demonstration:
# make random covariance, then independent samples from Gaussian
import numpy as np

ndim = 4
nsamp = 10000
random_state = np.random.default_rng(10)  # seed random generator
A = random_state.random((ndim, ndim))
cov = np.dot(A, A.T)
samps = random_state.multivariate_normal([0] * ndim, cov, size=nsamp)
A = random_state.random((ndim, ndim))
cov = np.dot(A, A.T)
samps2 = random_state.multivariate_normal([0] * ndim, cov, size=nsamp)

# Get the getdist MCSamples objects for the samples, specifying same parameter
# names and labels; if not specified weights are assumed to all be unity
names = ["x%s" % i for i in range(ndim)]
labels = ["x_%s" % i for i in range(ndim)]
samples = MCSamples(samples=samps, names=names, labels=labels)
samples2 = MCSamples(samples=samps2, names=names, labels=labels, label="Second set")

Removed no burn in
Removed no burn in

# Triangle plot (sometimes also called a corner plot)
g = plots.get_subplot_plotter()
g.triangle_plot([samples, samples2], filled=True)

# Here we are using inline plots, but if you wanted to export to file you'd just do e.g.
# g.export('output_file.pdf')

# 1D marginalized plot
g = plots.get_single_plotter(width_inch=4)
g.plot_1d(samples, "x2")

# 1D marginalized comparison plot
g = plots.get_single_plotter(width_inch=3)
g.plot_1d([samples, samples2], "x1")

# 1D normalized comparison plot
g = plots.get_single_plotter(width_inch=4)
g.plot_1d([samples, samples2], "x1", normalized=True)

# 2D line contour comparison plot with extra bands and markers
g = plots.get_single_plotter()
g.plot_2d([samples, samples2], "x1", "x2")
g.add_x_marker(0)
g.add_y_bands(0, 1)

# Filled 2D comparison plot with legend
g = plots.get_single_plotter(width_inch=4, ratio=1)
g.plot_2d([samples, samples2], "x1", "x2", filled=True)
g.add_legend(["sim 1", "sim 2"], colored_text=True);

# Shaded 2D comparison plot
g = plots.get_single_plotter(width_inch=4)
g.plot_2d([samples, samples2], "x1", "x2", shaded=True);

# Customized 2D filled comparison plot
g = plots.get_single_plotter(width_inch=6, ratio=3 / 5.0)
g.settings.legend_fontsize = 12
g.plot_2d([samples, samples2], "x1", "x2", filled=True, colors=["green", ("#F7BAA6", "#E03424")], lims=[-4, 7, -5, 5])
g.add_legend(["Sim ", "Sim 2"], legend_loc="upper right");

# Change the contours levels for marge stats and plots
# (note you need a lot of samples for 99% confidence contours to be accurate)
g = plots.get_single_plotter()
samples.updateSettings({"contours": [0.68, 0.95, 0.99]})
g.settings.num_plot_contours = 3
g.plot_2d(samples, "x1", "x2");

# 2D scatter (3D) plot
g = plots.get_single_plotter(width_inch=5)
g.plot_3d(samples, ["x1", "x2", "x3"])

# Multiple 1D subplots
g = plots.get_subplot_plotter(width_inch=5)
g.plots_1d(samples, ["x0", "x1", "x2", "x3"], nx=2);

# Multiple 2D subplots
g = plots.get_subplot_plotter(subplot_size=2.5)
g.settings.scaling = False  # prevent scaling down font sizes even though small subplots
g.plots_2d(samples, param_pairs=[["x0", "x1"], ["x2", "x3"]], nx=2, filled=True);

# Customized triangle plot
g = plots.get_subplot_plotter()
g.settings.figure_legend_frame = False
g.settings.alpha_filled_add = 0.4
g.settings.title_limit_fontsize = 14
g.triangle_plot(
    [samples, samples2],
    ["x0", "x1", "x2"],
    filled=True,
    legend_labels=["Simulation", "Simulation 2"],
    legend_loc="upper right",
    line_args=[{"ls": "--", "color": "green"}, {"lw": 2, "color": "darkblue"}],
    contour_colors=["green", "darkblue"],
    title_limit=1,  # first title limit (for 1D plots) is 68% by default
    markers={"x2": 0},
    marker_args={"lw": 1},
)

---

If you prefer to use one of the standard color series you can do that, using the matplotlib named colormaps

# tab10 is the standard discrete color table (good for color blindness etc.)
g = plots.get_subplot_plotter(subplot_size=2)

# Set line style default
g.settings.line_styles = "tab10"
g.plots_1d([samples, samples2], ["x1", "x2", "x3"], nx=3, legend_ncol=2, lws=[3, 2])

# or set explicitly
g.plots_1d([samples, samples2], ["x1", "x2", "x3"], nx=3, legend_ncol=2, colors="Set1", ls=["-", "--"])

# For filled contours set solid_colors setting (or set contour_colors argument as above)
g.settings.solid_colors = "tab10"
g.triangle_plot([samples, samples2], ["x1", "x2"], filled=True, contour_lws=2)

# 3D (scatter) triangle plot
g = plots.get_subplot_plotter(width_inch=6)
# you can adjust the scaling factor if font sizes are too small when
# making many subplots in a fixed size (default=2 would give smaller fonts)
g.settings.scaling_factor = 1
g.triangle_plot(
    [samples, samples2], ["x1", "x2", "x3"], plot_3d_with_param="x0", legend_labels=["Simulation", "Simulation 2"]
)

# You can reference g.subplots for manual tweaking,
# e.g. let's add a vertical axis line in the first column
# (this can also be done directly via the markers argument to triangle_plot)
for ax in g.subplots[:, 0]:
    ax.axvline(0, color="gray", ls="--")
IPython.display.display(g.fig)

# Rectangle 2D comparison plots
g = plots.get_subplot_plotter()
g.settings.figure_legend_frame = False
g.rectangle_plot(
    ["x0", "x1"],
    ["x2", "x3"],
    roots=[samples, samples2],
    filled=True,
    plot_texts=[["Test Label", None], ["Test 2", None]],
);

# Or join rows of 1D or 2D plots.
from matplotlib import cm

g.settings.linewidth = 2
g.plots_1d([samples, samples2], share_y=True, nx=2, legend_ncol=2, colors=cm.tab10)

g.rectangle_plot(["x1", "x2", "x3"], "x0", roots=[samples, samples2], colors=["k", "C2"]);

# Example of how to handle boundaries (samples are restricted to x0 >-0.5)
cut_samps = samps[samps[:, 0] > -0.5, :]
cut_samples = MCSamples(samples=cut_samps, names=names, labels=labels, ranges={"x0": (-0.5, None)}, label="Cut samples")
g = plots.get_subplot_plotter(subplot_size=2)
g.settings.title_limit_fontsize = 14  # reference size for 3.5 inch subplot
g.plots_1d(cut_samples, nx=4, title_limit=2)  # title by 95% limit
g = plots.get_single_plotter(width_inch=4, ratio=1)
g.plot_2d(cut_samples, "x0", "x1", filled=True);

Removed no burn in

# Periodic parameter: for samples within a perioidic range, you can also specify that by use a 3-tuple for the range
# This ensures the densities at the two sides agree, and reduces sampling noise near the boundaries 
# (but is not fully optimized for bandwidth selection)

def func(x): # perioid function
    return 0.8*np.sin(x+0.5)**2+0.2

rng = np.random.default_rng(10)
sample_array = rng.uniform(0, np.pi, 12500)
r = rng.uniform(0,1, len(sample_array))
p = []
for x, rand in zip(sample_array, r):
    if rand < func(x):
        p.append(x)

from getdist import MCSamples
samps = MCSamples(samples = np.array(p), names=['x'], ranges = {'x':[0, np.pi]})
samps2 = MCSamples(samples = np.array(p), names=['x'], ranges = {'x':[0, np.pi, True]})

gplot = plots.get_single_plotter(width_inch=4)
gplot.plot_1d(samps, 'x')
gplot.plot_1d(samps2, 'x', colors=['r'])
x = np.arange(0, np.pi, 0.01)
plt.plot(x,func(x), ls = '--')
plt.xlim(0, np.pi)
plt.ylim(0, None);
plt.legend(['Non-periodic','Periodic', 'True distribution']);

Removed no burn in
Removed no burn in

# Add and plot a new derived parameter
# getParms gets p so that p.x0, p.x1.. are numpy vectors of sample values
# For individual parameters you can also just do samples['x0'] etc.
p = samples.getParams()
assert np.all(p.x1 == samples["x1"])
samples.addDerived((5 + p.x2) ** 2, name="z", label="z_d")
g = plots.get_subplot_plotter(subplot_size=4)
g.plots_2d(samples, "x1", ["x2", "z"], nx=2);

# Example of how to do importance sampling, modifying the samples by re-weighting by a new likelihood
# e.g. to modify samples to be from the original distribution times a Gaussian in x1
# (centered on 1, with sigma=1.2)
# Using samples['x'] retrieves the vector of sample values for the parameter named 'x'
new_samples = samples.copy()  # make a copy so don't change the original
new_loglike = (samples["x1"] - 1) ** 2 / 1.2**2 / 2
# code currently assumes existing loglikes are set, set to zero here
new_samples.loglikes = np.zeros(samples.numrows)
# re-weight to account for the new likelihood
new_samples.reweightAddingLogLikes(new_loglike)
g = plots.get_single_plotter(width_inch=4, ratio=1)
g.plot_2d([samples, new_samples], "x0", "x1", filled=True);

# You can also account for general 2D prior boundary cuts that are not aligned with the axes
# (Note that the auto-smoothing kernel size does not account for non-trivial mask
# so may want to manually adjust the kernel smoothing width [smooth_scale_2D])

# e.g. consider cut that is linear function of two parameters with these parameters
y0 = -0.7
x0 = -0.2
r = 0.3


def mask_function(minx, miny, stepx, stepy, mask):
    # define function to tell getdist which 2D points are excluded by the prior
    # Note this should not include min, max parameter range cuts aligned with axes, which are handled as above.

    x = np.arange(mask.shape[1]) * stepx + minx
    y = np.arange(mask.shape[0]) * stepy + miny
    # Create 2D coordinate grids
    X, Y = np.meshgrid(x, y)
    # Zero out the array where prior excluded
    mask[Y < y0 - r * (X - x0)] = 0


cut_samps = samples.copy()
p = cut_samps.getParams()
cut_samps.filter(p.x1 - y0 + (p.x0 - x0) * r > 0)  # make test samples with hard prior cut

g = plots.get_single_plotter(width_inch=4, ratio=0.9)
# mass in the mask function so getdist knows about the prior cut
g.plot_2d(cut_samps, "x0", "x1", filled=True, mask_function=mask_function)
g.add_2d_contours(cut_samps, "x0", "x1", filled=False, color="g", ls="--")
x = np.linspace(-5, 5, 100)
plt.plot(x, y0 - r * (x - x0), color="k", ls="--")
g.add_legend(["prior mask corrected", "uncorrected"])
g.export("z:\\boundaries2D.pdf")

# Many other things you can do besides plot, e.g. get latex
# Default limits are 1: 68%, 2: 95%, 3: 99% probability enclosed
# See  https://getdist.readthedocs.io/en/latest/analysis_settings.html
# and examples for below for changing analysis settings
# (e.g. 2hidh limits, and how they are defined)
print(cut_samples.getInlineLatex("x0", limit=2))
print(samples.getInlineLatex("x0", limit=2))

x_0 < 2.48
x_0 = 0.0^{+2.6}_{-2.6}

print(samples.getInlineLatex("x1", limit=1))

x_1 = 0.0\pm 1.2

print(samples.getTable().tableTex())

\begin{tabular} { l  c}

 Parameter &  95\% limits\\
\hline
{\boldmath$x_0            $} & $0.0^{+2.6}_{-2.6}         $\\

{\boldmath$x_1            $} & $0.0^{+2.4}_{-2.4}         $\\

{\boldmath$x_2            $} & $0.0^{+2.7}_{-2.7}         $\\

{\boldmath$x_3            $} & $0.0^{+3.1}_{-3.1}         $\\

$z_d                       $ & $27^{+30}_{-20}            $\\
\hline
\end{tabular}

# results from multiple chains
from getdist.types import ResultTable

print(
    ResultTable(
        ncol=1, results=[samples, new_samples], paramList=["x0", "x3"], limit=1, titles=["Samples", "Weighted samples"]
    ).tableTex()
)

\begin{tabular} { l  c c}
\noalign{\vskip 3pt}\hline\noalign{\vskip 1.5pt}\hline\noalign{\vskip 5pt}
 \multicolumn{1}{c}{\bf } &  \multicolumn{1}{c}{\bf Samples} &  \multicolumn{1}{c}{\bf Weighted samples}\\
\noalign{\vskip 3pt}\cline{2-3}\noalign{\vskip 3pt}

 Parameter &  68\% limits &  68\% limits\\
\hline
{\boldmath$x_0            $} & $0.0\pm 1.3                $ & $0.4\pm 1.1                $\\

{\boldmath$x_3            $} & $0.0\pm 1.6                $ & $0.6\pm 1.2                $\\
\hline
\end{tabular}

print(samples.PCA(["x1", "x2"]))

PCA for parameters:
         2 :x_1
         3 :x_2

Correlation matrix for reduced parameters
          x1 :  1.0000  0.7303
          x2 :  0.7303  1.0000

e-values of correlation matrix
PC 1:   0.2697
PC 2:   1.7303

e-vectors
  2: -0.7071  0.7071
  3:  0.7071  0.7071

Principal components
PC1 (e-value: 0.269710)
[1.229952]   (x_1-0.000555)/1.000000
[-1.229952]   (x_2-0.003484)/-1.137321
          = 0.000000 +- 0.903342

PC2 (e-value: 1.730290)
[1.229952]   (x_1-0.000555)/1.000000
[1.229952]   (x_2-0.003484)/1.137321
          = 0.000000 +- 2.288034

Correlations of principal components
       1       2
PC 1   1.000   0.000
PC 2   0.000   1.000
   1   0.029   0.834   (x_0)
   2   0.367   0.930   (x_1)
   3  -0.367   0.930   (x_2)
   4   0.018   0.990   (x_3)
   5  -0.363   0.912   (z_d)

stats = cut_samples.getMargeStats()
lims0 = stats.parWithName("x0").limits
lims1 = stats.parWithName("x1").limits
for conf, lim0, lim1 in zip(samples.contours, lims0, lims1):
    print("x0 %s%% lower: %.3f upper: %.3f (%s)" % (conf, lim0.lower, lim0.upper, lim0.limitType()))
    print("x1 %s%% lower: %.3f upper: %.3f (%s)" % (conf, lim1.lower, lim1.upper, lim1.limitType()))

x0 0.68% lower: -0.500 upper: 1.063 (one tail upper limit)
x1 0.68% lower: -0.493 upper: 1.491 (two tail)
x0 0.95% lower: -0.500 upper: 2.481 (one tail upper limit)
x1 0.95% lower: -1.382 upper: 2.532 (two tail)
x0 0.99% lower: -0.500 upper: 3.257 (one tail upper limit)
x1 0.99% lower: -1.792 upper: 3.371 (two tail)

# if samples have likelihood values, can also get best fit sample and extremal values of N-D confidence region
# Note in high dimensions best-fit sample is likely a long way from the best fit; N-D limits also often MC-noisy
print(new_samples.getLikeStats())
print(
    "x0 95% n-D confidence extrema:",
    new_samples.paramNames.parWithName("x0").ND_limit_bot[1],
    new_samples.paramNames.parWithName("x0").ND_limit_top[1],
)

Best fit sample -log(Like) = 0.000000
Ln(mean 1/like) = 0.529321
mean(-Ln(like)) = 0.337642
-Ln(mean like)  = 0.259510
2*Var(Ln(like)) = 0.438724

parameter   bestfit        lower1         upper1         lower2         upper2
x0          9.0893191E-01 -2.1779183E+00  3.8679700E+00 -3.1333809E+00  4.1960094E+00   x_0
x1          1.0000762E+00  1.5285567E-02  1.9829577E+00 -9.3128975E-01  2.9312813E+00   x_1
x2          1.9388624E+00 -2.8283015E+00  4.7720226E+00 -3.2591460E+00  5.1498794E+00   x_2
x3          1.7908065E+00 -1.6800034E+00  4.1609904E+00 -2.9384466E+00  4.8321531E+00   x_3
z*          4.8147812E+01  4.7162744E+00  9.5492426E+01  3.0305726E+00  1.0302005E+02   z_d

x0 95% n-D confidence extrema: -3.133380940929877 4.196009420129642

# Save to file
import os
import tempfile

tempdir = os.path.join(tempfile.gettempdir(), "testchaindir")
if not os.path.exists(tempdir):
    os.makedirs(tempdir)
rootname = os.path.join(tempdir, "testchain")
samples.saveAsText(rootname)

# Load from file
from getdist import loadMCSamples

readsamps = loadMCSamples(rootname)

C:\Users\micro\AppData\Local\Temp\testchaindir\testchain.txt
Removed no burn in

# Make plots from chain files, loading automatically as needed by using root file name
g = plots.get_single_plotter(chain_dir=tempdir, width_inch=4)
g.plot_2d("testchain", "x1", "x2", shaded=True);

# Custom settings for all loaded chains can be set as follows;
# for example to use custom contours and remove the first 20% of each chain as burn in
g = plots.get_single_plotter(
    chain_dir=tempdir, analysis_settings={"ignore_rows": 0.2, "contours": [0.2, 0.4, 0.6, 0.8]}
)
g.settings.num_plot_contours = 4
g.plot_2d("testchain", "x1", "x2", filled=False);

C:\Users\micro\AppData\Local\Temp\testchaindir\testchain.txt
Removed 0.2 as burn in

# Silence messages about load
getdist.chains.print_load_details = False

# Chains can be loaded by searching in multiple directories by giving a list as chain_dir
# (note chain names must be unique)

# make second test chain in new temp dir
temp2 = os.path.join(tempdir, "chaindir2")
cut_samples.saveAsText(os.path.join(temp2, "testchain2"), make_dirs=True)
# Plot from chain files
g = plots.get_single_plotter(chain_dir=[tempdir, temp2])
g.plot_2d(["testchain", "testchain2"], "x1", "x2", filled=True);

# You can also load a samples object from the chain directories for further manipulation
read_samples = g.samples_for_root("testchain")
# e.g. add a new derived parameter
# Note our new variable is >0 by definition, so also need to define its finite range to
# get correct densities near zero
p = read_samples.getParams()
read_samples.addDerived(np.abs(p.x1), name="modx1", label="|x_1|", range=[0, None])
g.new_plot()
g.plot_2d(read_samples, "modx1", "x2", filled=True);

# cleanup test files
import shutil

shutil.rmtree(tempdir)

# The plotting scripts also let you plot Gaussian (or Gaussian mixture) contours
# This is useful for plotting smooth theory results, e.g. Fisher forecasts.
from getdist.gaussian_mixtures import GaussianND

covariance = [[0.001**2, 0.0006 * 0.05, 0], [0.0006 * 0.05, 0.05**2, 0.2**2], [0, 0.2**2, 2**2]]
mean = [0.02, 1, -2]
gauss = GaussianND(mean, covariance)
g = plots.get_subplot_plotter()
g.triangle_plot(gauss, filled=True)

# You can also explicitly name parameters so Gaussian mixtures can be plotted
# in combination with samples
from getdist.gaussian_mixtures import Mixture2D

cov1 = [[0.001**2, 0.0006 * 0.05], [0.0006 * 0.05, 0.05**2]]
cov2 = [[0.001**2, -0.0006 * 0.05], [-0.0006 * 0.05, 0.05**2]]
mean1 = [0.02, 0.2]
mean2 = [0.023, 0.09]
mixture = Mixture2D([mean1, mean2], [cov1, cov2], names=["zobs", "t"], labels=[r"z_{\rm obs}", "t"], label="Model")

# Generate samples from the mixture as simple example
mix_samples = mixture.MCSamples(3000, label="Samples")

g = plots.get_subplot_plotter()
# compare the analytic mixture to the sample density
g.triangle_plot([mix_samples, mixture], filled=False)

# For long tick labels GetDist will space them intelligently so they don't overlap.
# You can also rotate tick labels
p = mix_samples.getParams()
mix_samples.addDerived(10 * p.zobs - p.t / 10, name="xt", label="x_t")

g = plots.get_subplot_plotter(subplot_size=3)
g.settings.axes_fontsize = 12
g.settings.axis_tick_x_rotation = 45
g.settings.axis_tick_y_rotation = 90
g.settings.colorbar_tick_rotation = 90
g.triangle_plot(mix_samples, ["t", "zobs"], plot_3d_with_param="xt", filled=True)

# Double triangle plot
g = plots.get_subplot_plotter()
g.triangle_plot(
    [samples, new_samples],
    ["x0", "x1", "x2"],
    filled=True,
    upper_roots=[samples2],
    upper_kwargs={"contour_colors": ["green"]},
    legend_labels=["Samples", "Reweighted samples", "Second samples"],
);

# Variants
from matplotlib import cm

g = plots.get_subplot_plotter()

upper_kwargs = {
    "contour_colors": cm.tab10.colors[4:],
    "contour_ls": ["-", "--"],
    "filled": [True, False],
    "show_1d": [True, True],
    "contour_lws": [1, 2],
}
g.settings.solid_contour_palefactor = 0.9
g.settings.alpha_filled_add = 0.6
g.triangle_plot(
    [samples, new_samples],
    ["x0", "x1", "x2"],
    filled=True,
    contour_colors=["C3", "green"],
    markers={"x0": -0.5},
    upper_roots=[samples2, cut_samples],
    upper_kwargs=upper_kwargs,
    upper_label_right=True,
    legend_labels=["Samples", "Reweighted", "Second", "Cut"],
);

# Sets of 3D plots
g = plots.get_subplot_plotter(width_inch=8)
g.plots_3d(mix_samples, [["xt", "t", "zobs"], ["t", "zobs", "xt"]], nx=2)
# and this is how to manually add contours only to a specific plot
g.add_2d_contours(mixture, "t", "zobs", ax=[0, 1]);

# 4D x-y-z-color scatter plots:
# Use ""%matplotlib notebook" to interactively rotate etc.

g = plots.get_single_plotter()
g.plot_4d(
    [samples, samples2],
    ["x0", "x1", "x2", "x3"],
    cmap="viridis",
    color_bar=False,
    azim=75,
    alpha=[0.3, 0.1],
    shadow_color=False,
    compare_colors=["k"],
)

# with projection onto axes, colorbar, custom limits, colors, etc.
g = plots.get_single_plotter()
g.plot_4d(
    samples,
    ["x0", "x1", "x2", "x3"],
    cmap="jet",
    alpha=0.4,
    shadow_alpha=0.05,
    shadow_color=True,
    max_scatter_points=6000,
    lims={"x2": (-3, 3), "x3": (-3, 3)},
    colorbar_args={"shrink": 0.6},
)

# use animate=True and mp4_filename option to export an rotation animation
# this shows example output

from IPython.display import Video

Video("https://cdn.cosmologist.info/antony/sample_rotation.mp4", html_attributes="controls loop", width=600)

---

Using styles to change settings for consistent plot scripts

If you want to change several default settings, you can make a module containing a new plotter class, which defines its own settings and behaviour, and then use it as a style. A couple of sample styles are included: getdist.styles.tab10 and getdist.styles.planck. Using the same style in each script will then give consistent outputs without changing settings each time.

# use the 'tab10' style, which uses default matplotlib colourmaps
from getdist.styles.tab10 import style_name

plots.set_active_style(style_name)

g = plots.get_subplot_plotter(width_inch=6)
g.triangle_plot(
    [samples, samples2], ["x1", "x2", "x3"], plot_3d_with_param="x0", legend_labels=["Simulation", "Simulation 2"]
)