Three different ways to do parameter samplingΒΆ

This code snipped shows how to use LensTools to perform cosmological parameter sampling with three diffenent methods: direct evaluation of the likelihood, Fisher Matrix and MCMC

from lenstools.statistics.ensemble import Series,Ensemble
from lenstools.statistics.constraints import Emulator
from lenstools.statistics.contours import ContourPlot

import numpy as np
import matplotlib.pyplot as plt


def lt_sample(emulator,test_data,covariance,p_value=0.68):

        #Check that the data types are correct
        assert isinstance(emulator,Emulator)
        assert isinstance(test_data,Series)
        assert isinstance(covariance,Ensemble)

        #Plot setup
        fig,ax = plt.subplots(figsize=(8,8))


        #Map the likelihood in the OmegaM-sigma8 plane, fix w to -1
        p = Ensemble.meshgrid({
                "Om":np.linspace(0.2,0.5,50),
                "sigma8":np.linspace(0.6,0.9,50)
        })

        p["w"] = -1.

        #Compute the chi squared scores of the test data on a variety of parameter points
        scores = emulator.score(p,test_data,covariance,correct=1000,method="chi2")
        scores["likelihood"] = np.exp(-0.5*scores[emulator.feature_names[0]])

        contour = ContourPlot.from_scores(scores,parameters=["Om","sigma8"],
                feature_names=["likelihood"],
                plot_labels=[r"$\Omega_m$",
                r"$\sigma_8$"],fig=fig,ax=ax)

        contour.getLikelihoodValues([p_value],precision=0.01)
        contour.plotContours(colors=["red"])
        contour.labels()

        #Approximate the emulator linearly around the maximum (Fisher matrix)
        fisher = emulator.approximate_linear(center=(0.26,-1.,0.8))

        #Consider (OmegaM,sigma8) only
        fisher.pop(("parameters","w"))
        fisher = fisher.iloc[[0,1,3]]

        #Fisher confidence ellipse
        ellipse = fisher.confidence_ellipse(covariance,
                correct=1000,
                observed_feature=test_data,
                parameters=["Om","sigma8"],
                p_value=p_value,
                fill=False,
                edgecolor="blue")

        ax.add_artist(ellipse)

        #MCMC sampling of (OmegaM,sigma8)
        samples = emulator.sample_posterior(test_data,
                features_covariance=covariance,
                correct=1000,
                pslice={"w":-1},
                sample="emcee")[emulator.feature_names[0]]

        ax.scatter(samples["Om"],samples["sigma8"],marker=".",color="black",s=1)
        ax.set_xlim(0.2,0.5)
        ax.set_ylim(0.6,0.9)

        #Save the figure
        fig.tight_layout()
        fig.savefig("parameter_sampling.png")

And this is the resulting figure:

../_images/parameter_sampling.png