Introduction to the likelihood

When we are fitting the parameters $\alpha$ in some model $\mathcal{M}(x;\alpha)$ to some data $\mathcal{D} = {x_i , y_i}$, how do we judge how well a given model prediction ${ x_i , y_m(x_i;\alpha) \equiv \mathcal{M}(x_i;\alpha) } $ describes our data $\mathcal{D}$?

We have to come up with a model for the likelihood. This is the conditional probability that our data $\mathcal{D}$ (also called evidence) is described by a given model prediction $y_m(x_i;\alpha)$. We can write it as $p(\mathcal{D} | y_m(x_i;\alpha))$ or $p(\mathcal{D} | \mathcal{M}, \alpha)$. Sometimes we use $\mathcal{L}$ instead of $p$ to write the likelihood, but it is just a conditional probability. It tells you the probability, given model $\mathcal{M}$, with a fixed set of parameters $\alpha$, that some data $\mathcal{D} = { x_i,y_i }$ is the result.

Every method for fitting or calibrating a model (whether it is explicitly stated or not) makes some assumption about the likelihood.

How do we come up with a defensible form for the likelihood for a given set of evidence and a given model?

A very simple likelihood model

What if there are no errors in our model or in ${ x_i, y_i}$ at all? Then our model must exactly reproduce $y_i$, and our lileihood is:

\begin{equation} p({ x_i,y_i } | y_m(x_i;\alpha) ) = \begin{cases} 1 & y_m(x_i;\alpha) = y_i \ 0 & \rm{otherwise} \end{cases} \end{equation}

Can you guess what issues this might have?

The general case

In real cases, our evidence $\mathcal{D}$ is more than just a set of numbers ${x_i , y_i}$. It is, in fact, a probability distribution itself: $\mathcal{D} \equiv p( \vec{x}, \vec{y} )$; a given ${x,y}$ pair is a measurement of a random variable. Based on the uncertainty information provided by the experimentalists, and the limitations in our model, we must come up with a reasonable model for what $p( \vec{x}, \vec{y} )$ is, and use it to construct a likelihood.

Then, when we do a frequentist model fit, we are searching through $\alpha$-space trying to find the $\alpha$ that corresponds to $\text{max} \left[ p(\mathcal{D}| \mathcal{M},\alpha) \right]$. This is called Maximum Likelihood Estimation (MLE).

As we shall see, minimizing the $\chi^2$ is just a special case of MLE, subject to certain assumptions (which may not always be appropriate!).

When we do a Bayesian calibration, we sample $\alpha$-space in such a way that the samples converge on a distribution that is related to the likelihood, called the posterior $p(\alpha | \mathcal{D}, \mathcal{M})$, which just modifies the likelihood to include prior belief about $\alpha$:

\begin{equation} p(\alpha | \mathcal{D}, \mathcal{M}) \propto p(\mathcal{D}| \mathcal{M},\alpha) p(\alpha). \end{equation}

This is just the result of Bayes theorem. In Bayesian calibration as well, one often sees the likelihood $p(\mathcal{D}| \mathcal{M},\alpha)$ modeled as $\exp(-\chi^2)$. Again, this may not always be appropriate!

Under what assumptions is a $\chi^2$-distribution appropriate for the likelihood?

Let us assume there is some ground truth ${x_i, y_i^\rm{true}}$, and futhermore assume that our model, for some unknown true parameters, is able to exactly replicate the truth:

\begin{equation} M(x_i;\alpha^\rm{true}) = y_i^\rm{true} \end{equation}

This may be a big assumption, and, in realistic scenarios, one may want to to instead use

\begin{equation} M(x_i;\alpha^\rm{true}) = y_i^\rm{true} + \epsilon_i \end{equation}

Where $\epsilon_i$ is another random variable describing the model error. We will ignore this for now.

Now we need more assumptions: the measured values $\vec{x}, \vec{y}$ are random vectors from some multivariate probability distribution $p(\vec{x},\vec{y})$.

Let’s assume:

• there is no error in $x$; $x_i = x_i^\rm{true}$

• each $y_i$ is really a number of counts in the bin $x_i = \left[ x_{i0} , x_{if} \right]$

• each count is independent of all the others (this is akin to saying there is no systematic error; later in the demo we will see a case for which this assumption is disastrous!)

If there are $N$ total counts, then $p_i = y_i^\rm{true}/N$ is the probability of a count being in bin $i$, and $ 1 - p_i$ is the probability of not being in bin $i$. Our model prediction is for the probability of a single count being in bin $i$ is just $y_m(x_i;\alpha) / N$. The probability of getting exactly $y$ counts in bin $i$ in $N$ trials is a binomial distribution:

\begin{equation} p(y_i | N, p_i ) = B(N,p_i) \equiv \binom{N}{y_i} p_i^{y_i} (1-p_i)^{N-y_i} \end{equation}

Plugging in our model prediction for $p_i$, we have the likelihood

\begin{equation} p(y_i | N, p_i = \mathcal{M}(x_i;\alpha)/N ) = B(N,\mathcal{M}(x_i;\alpha)/N) \equiv \binom{N}{y_i} \left( \frac{y_m(x_i;\alpha)}{N} \right)^{y_i} \left( 1- \frac{y_m(x_i;\alpha)}{ N} \right)^{N-y_i} \end{equation}

The total liklelihood is the product of this probability for every bin:

\begin{equation} p( \mathcal{D} | N, \mathcal{M} , \alpha)) = \prod_i B(N,\mathcal{M}(x_i;\alpha)/N) \equiv \prod_i \binom{N}{y_i} \left( \frac{y_m(x_i;\alpha)}{N} \right)^{y_i} \left( 1- \frac{y_m(x_i;\alpha)}{ N} \right)^{N-y_i} \end{equation}

As long as $N$ is known exactly (it’s own can of worms), then this is exactly what we are looking for! In principle, one could exactly use this formula in this situation.

How do we get the $\chi^2$ form? Well, when $N \rightarrow \infty$ and $p_i^\rm{true} \rightarrow 0$ such that $N p_i^\rm{true} = y_i^\rm{true}$ stays constant, the binomial distribution limits to a Poission distribution, with mean $y_i^{\rm{true}}$ and standard deviation $\sigma_i = \sqrt{y_i^{\rm{true}}}$:

\begin{equation} p( \mathcal{D} | N, \mathcal{M} , \alpha) \rightarrow \prod_i \frac{ (y_i^{\rm{true}})^y }{ y! } e^{-y_i^{\rm{true}}} \end{equation}

Then using the Central Limit Theorem, we find in the limit as $N \rightarrow \infty$ and $0 < p_i < 1$ fixed, that our Poisson distribution becomes a normal distribution:

\begin{align} p( \mathcal{D} | N, \mathcal{M} , \alpha) &\rightarrow \prod_i \frac{1}{\sqrt{2 \pi y_i^{\rm{true}}}} \exp{ - \frac{\left( y_i - y_i^{\rm{true}}\right)^2}{ 2 y_i^{\rm{true}}} } \end{align}

Plugging in $y_m(x_i;\alpha)$ for $y_i^\rm{true}$ and $\sigma_i$ for the standard deviation, we have:

\begin{equation} p( \mathcal{D} | N, \mathcal{M} , \alpha) = \prod_i \frac{1}{\sqrt{2 \pi \sigma_i^2}} \exp{ \left( \frac{(y_i - y_m(x_i;\alpha))^2}{2 \sigma_i^2} \right)} \end{equation}

Take the log of this function, one of the terms will be of the $\chi^2$ form:

\begin{equation} \log p( \mathcal{D} | N, \mathcal{M} , \alpha) \propto -\frac{1}{2} \sum_i \frac{(y_i - y_m(x_i;\alpha))^2}{\sigma_i^2} + \dots \equiv \chi^2 + \dots \end{equation}

It is a useful excercise to take this log and see what the other terms are exactly. In some likelihood models, where the covariance is not known a priori but has parameters that are fit along side the physical model, these other terms will also play a role.

If you’ve taken a stats course, you’ve probably done these two limits (Binomial to Poission and Poission to Normal) at some point or another. A nice pedagogical source for the derivations is The Knolly Bible, Ch. 3. Also, I would like to mention that chatgpt was also useful in preparing this discussion.

too long; didn’t read:

The main point I want to get across is is that we had to make many assumptions to get our likeliood to look like a $\chi^2$! These assumptions are not valid in every case. Whether you are doing Bayesian calibration or frequentist model fitting, whether or not you’re even doing uncertainty quantification (which you should be doing), or you just care about finding the “best” parameters, it is your job to come up with a likelihood model. That is, a model for the function $p({x_i,y_i}| \mathcal{M},\alpha)$. You must carefully consider the data and its errors, as well as your physics model and its errors, to come up with such a model. You must also think carefully about the implications of any simplifying assumptions you make (e.g. ignoring systematic error), and transparently disclose them.

In this demo: the multivariate normal liklelihood model

In this demo, we will introduce some systematic error. That is, instead of measuring the number of counts in a bin, $y_i$ will instead be some derived quantity like a cross section, which includes some normalization. In particular, we will consider a scenario where this normalization is not known exactly.

This will break the assumption from above that each $y_i$ is indepent from each other. If the limiting conditions leading to a normal distribution are still valid, but there is some covariance between $y_i$’s for different points $x_i$, one can write the likelihood more generally as:

\begin{equation} p( \mathcal{D} | N, \mathcal{M} , \alpha) = \frac{1}{\sqrt{(2 \pi)^k |\mathbf{\Sigma}|}} \exp{\left( -\frac{1}{2} \vec{\Delta}^T \cdot \mathbf{\Sigma}^{-1} \cdot \vec{\Delta} \right)}, \end{equation}

where $\Delta \equiv (y_i - y_m(x_i;\alpha))$ and $k$ is the number of elements in the vector $\vec{\Delta}$. Then, one must find a model for the covariance matrix $\mathbf{\Sigma}$. If $\mathbf{\Sigma}$ is diagonal, this will resemble the $\exp{\left( -\frac{1}{2}\chi^2 \right)}$ form. The demo will begin with this assmumption, and test out different models for the covariance matrix $\mathbf{\Sigma}$.

In fact, one can show that this is a very good assumption for many cases due to the Central Limit Theorem, as any bounded, finite-variance distribution looks like a normal distribution near the mean. One can think of the generalized $\chi^2 \equiv \vec{\Delta}^T \cdot \mathbf{\Sigma}^{-1} \cdot \vec{\Delta}$ (which is also the squared Mahalanobis distance) as the truncation of the Taylor expansion of the log probability of an arbitrary distribution to 2nd order. See Bayesian Data Analysis by Gelman, Ch. 4. In other words, for cases in which there are many samples and the CLT applies, the above model for the likelihood function as a generalized multivariate normal is often adequate. The job of the modeler is then to come up with a model for $\mathbf{\Sigma}$, given their knowledge of the experimental evidence and the model. This is what we will be doing in the demo.

Some good reading is D’Agostini, 1994.

Comparing likelihood models: there is a right way, and many wrong ways!

from collections import OrderedDict

import corner
import numpy as np
from matplotlib import pyplot as plt
from scipy import stats

import rxmc

Plotting functions

def plot_chains(walker, model, true_params):
    fig, axes = plt.subplots(
        walker.model_sampler.chain.shape[1] + 1, 1, figsize=(8, 8), sharex=True
    )
    for i in range(walker.model_sampler.chain.shape[1]):
        axes[i].plot(walker.model_sampler.chain[:, i])
        axes[i].set_ylabel(f"${model.params[i].latex_name}$ [{model.params[i].unit}]")
        true_value = true_params[model.params[i].name]
        axes[i].hlines(
            true_value, 0, len(walker.model_sampler.chain), "r", linestyle="--"
        )

    axes[-1].plot(walker.model_sampler.logp_chain)
    axes[-1].set_ylabel(r"$\log{\mathcal{L}(\alpha_i | \mathcal{O})}$")
    axes[-1].set_xlabel(r"$i$")

def plot_posterior_corner(walker, true_params):
    fig = corner.corner(
        walker.model_sampler.chain,
        labels=[p.name for p in my_model.params],
        label="posterior",
        truths=[true_params["m"], true_params["b"]],
    )
    fig.suptitle("posterior")

def plot_predictive_post(walker, model, x, y_exp, y_err, y_true, x_true=None):
    if x_true is None:
        x_true = x
    n_posterior_samples = walker.model_sampler.chain.shape[0]
    y = np.zeros((n_posterior_samples, len(x)))
    for i in range(n_posterior_samples):
        sample = walker.model_sampler.chain[i, :]
        y[i, :] = model.y(x, *sample)

    upper, median, lower = np.percentile(y, [5, 50, 95], axis=0)
    plt.errorbar(
        x,
        y_exp,
        y_err,
        color="k",
        marker="o",
        linestyle="none",
        label="experiment",
    )

    plt.plot(x, y_true, "k--", label="truth")
    plt.plot(x, median, "m:", label="posterior median")
    plt.fill_between(
        x,
        lower,
        upper,
        alpha=0.5,
        zorder=2,
        label=r"posterior inner 90$^\text{th}$ pctl",
    )
    plt.legend()
    plt.xlabel("x")
    plt.ylabel("y")

make the model

class LinearModel(rxmc.physical_model.PhysicalModel):
    def __init__(self):
        params = [
            rxmc.params.Parameter("m", float, "no-units"),
            rxmc.params.Parameter("b", float, "y-units"),
        ]
        super().__init__(params)

    def y(self, x, m, b):
        return m * x + b

    def evaluate(self, observation, m, b):
        return self.y(observation.x, m, b)

my_model = LinearModel()

rng = np.random.default_rng(16)

set up prior

true_params = OrderedDict(
    [
        ("m", 0.6),
        ("b", 2),
    ]
)

prior_mean = OrderedDict(
    [
        ("m", 2),
        ("b", 5),
    ]
)
prior_std_dev = OrderedDict(
    [
        ("m", 2),
        ("b", 2),
    ]
)

covariance = np.diag(list(prior_std_dev.values())) ** 2
mean = np.array(list(prior_mean.values()))
prior_distribution = stats.multivariate_normal(mean, covariance)

systematic_fractional_err = 0.1
# choose a normalization 1 std deviation below the mean
N = 1 - systematic_fractional_err
noise_fraction = 0.05
x = np.linspace(0.01, 1.0, 15, dtype=float)
y_true = my_model.y(x, *list(true_params.values()))
y_exp = (y_true + rng.normal(scale=noise_fraction * y_true, size=len(x))) * N
y_stat_err = noise_fraction * y_exp * N

plt.errorbar(
    x,
    y_exp,
    y_stat_err,
    color="k",
    marker="o",
    linestyle="none",
    label="experiment with bias",
)
plt.plot(x, y_true, "k--", label="truth")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.title("experimental constraint with bias")

Text(0.5, 1.0, 'experimental constraint with bias')

Compare Likelihood Models

We will look at a few different cases:

1. Covariance is fixed to just statistical error (disregarding systematic error)

   \begin{equation} \Sigma_{ij} = \delta_{ij} \sigma_{stat,i}^2 \end{equation}

2. Covariance is just statistical error, but we fit the magnitude of the statistical noise $\eta$ (disregarding systematic error). This means the covariance is not fixed, but will be updated during calibration.

   \begin{equation} \Sigma_{ij} = \delta_{ij} \eta^2 y_m(x_j; \alpha)^2 \end{equation}

3. Systematic error is included properly in covariance, making the covariance a function of the model prediction. Again, this means the covariance is not fixed, but will be updated during calibration.

   \begin{equation} \Sigma_{ij}(\alpha) = \delta_{ij} \sigma_{stat,i}^2 + \sigma_N^2 y_m(x_i; \alpha) y_m(x_j; \alpha) \end{equation}

4. Systematic error is included improperly in covariance, using the experimental $y(x_i)$ instead of the model prediction $y_m(x_i;\alpha)$. In this case the covariance is again fixed.

   \begin{equation} \Sigma_{ij} = \delta_{ij} \sigma_{stat,i}^2 + \sigma_N^2 y(x_i) y(x_j) \end{equation}

# 1 and 2
obs_stat_only = rxmc.observation.Observation(
    x=x,
    y=y_exp,
    y_stat_err=y_stat_err,
)

# 3
obs_sys_norm_correct = rxmc.observation.Observation(
    x=x,
    y=y_exp,
    y_stat_err=y_stat_err,
    y_sys_err_normalization=systematic_fractional_err,
)

# 4
obs_sys_norm_wrong = rxmc.observation.FixedCovarianceObservation(
    x=x,
    y=y_exp,
    covariance=np.diag(y_stat_err**2)
    + systematic_fractional_err**2 * np.outer(y_exp, y_exp),
)

set up likelihood models and constraints

# 1 and 3
likelihood = rxmc.likelihood_model.LikelihoodModel()

# 4
likelihood_fixed_cov = rxmc.likelihood_model.FixedCovarianceLikelihood()


# 2 - a special likelihood model that takes in the noise fraction as a parameter
likelihood_unknown_stat = rxmc.likelihood_model.UnknownNoiseFractionErrorModel()

# 1
evidence_stat_only = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_stat_only],
            my_model,
            likelihood,
        )
    ]
)

# 2
evidence_unknown_stat = rxmc.evidence.Evidence(
    constraints=[],
    parametric_constraints=[
        rxmc.constraint.Constraint(
            [obs_stat_only],
            my_model,
            likelihood_unknown_stat,
        )
    ],
)

# 3
evidence_sys_correct = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_correct],
            my_model,
            likelihood,
        )
    ]
)

# 4
evidence_sys_wrong = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_wrong],
            my_model,
            likelihood_fixed_cov,
        )
    ]
)

def proposal_distribution_model(x, rng):
    return stats.multivariate_normal.rvs(
        mean=x, cov=prior_distribution.cov / 1000, random_state=rng
    )

Run option 1: fixed covariance, statistical error only

walker1 = rxmc.walker.Walker(
    rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence_stat_only,
    rng=rng,
)

%%time
walker1.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.336
CPU times: user 3.33 s, sys: 124 ms, total: 3.45 s
Wall time: 3.46 s

plot_chains(walker=walker1, model=my_model, true_params=true_params)

plot_posterior_corner(walker=walker1, true_params=true_params)

plot_predictive_post(
    walker=walker1, model=my_model, x=x, y_exp=y_exp, y_err=y_stat_err, y_true=y_true
)
plt.title("option 1: fixed statistical error, systematic ignored")

Text(0.5, 1.0, 'option 1: fixed statistical error, systematic ignored')

Run option 2: unknown statistical error

We need to come up with a prior for the noise for option 2. We will keep it fairly wide and centered about the reported value.

noise_prior = stats.norm(loc=np.log(noise_fraction), scale=1)

def proposal_distribution_noise(x, rng):
    return np.atleast_1d(stats.norm.rvs(loc=x, scale=0.1, random_state=rng))

walker2 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_unknown_stat,
    likelihood_samplers=[
        rxmc.param_sampling.MetropolisHastingsSampler(
            params=likelihood_unknown_stat.params,
            starting_location=noise_prior.mean(),
            proposal=proposal_distribution_noise,
            prior=noise_prior,
        )
    ],
    rng=rng,
)

%%time
walker2.walk(
    n_steps=10000,
    burnin=1000,
    batch_size=100,
)

Burn-in batch 1/10 completed, 100 steps.
Burn-in batch 2/10 completed, 100 steps.
Burn-in batch 3/10 completed, 100 steps.
Burn-in batch 4/10 completed, 100 steps.
Burn-in batch 5/10 completed, 100 steps.
Burn-in batch 6/10 completed, 100 steps.
Burn-in batch 7/10 completed, 100 steps.
Burn-in batch 8/10 completed, 100 steps.
Burn-in batch 9/10 completed, 100 steps.
Burn-in batch 10/10 completed, 100 steps.
Batch: 1/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.8]
Batch: 2/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.78]
Batch: 3/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.86]
Batch: 4/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.260
  Likelihood parameter acceptance fractions: [0.83]
Batch: 5/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.82]
Batch: 6/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.82]
Batch: 7/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.81]
Batch: 8/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.79]
Batch: 9/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.81]
Batch: 10/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.73]
Batch: 11/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.85]
Batch: 12/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.88]
Batch: 13/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.91]
Batch: 14/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.82]
Batch: 15/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.79]
Batch: 16/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.85]
Batch: 17/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.86]
Batch: 18/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.84]
Batch: 19/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.540
  Likelihood parameter acceptance fractions: [0.89]
Batch: 20/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.79]
Batch: 21/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.79]
Batch: 22/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.9]
Batch: 23/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.82]
Batch: 24/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.8]
Batch: 25/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.78]
Batch: 26/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.89]
Batch: 27/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.82]
Batch: 28/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.83]
Batch: 29/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.200
  Likelihood parameter acceptance fractions: [0.86]
Batch: 30/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.86]
Batch: 31/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.85]
Batch: 32/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.220
  Likelihood parameter acceptance fractions: [0.85]
Batch: 33/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.85]
Batch: 34/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.8]
Batch: 35/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.8]
Batch: 36/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.82]
Batch: 37/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.83]
Batch: 38/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.220
  Likelihood parameter acceptance fractions: [0.9]
Batch: 39/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.76]
Batch: 40/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.86]
Batch: 41/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.88]
Batch: 42/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.83]
Batch: 43/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.81]
Batch: 44/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.82]
Batch: 45/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.83]
Batch: 46/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.84]
Batch: 47/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.84]
Batch: 48/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.86]
Batch: 49/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.86]
Batch: 50/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.84]
Batch: 51/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.79]
Batch: 52/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.89]
Batch: 53/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.220
  Likelihood parameter acceptance fractions: [0.83]
Batch: 54/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.250
  Likelihood parameter acceptance fractions: [0.85]
Batch: 55/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.74]
Batch: 56/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.660
  Likelihood parameter acceptance fractions: [0.82]
Batch: 57/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.85]
Batch: 58/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.85]
Batch: 59/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.87]
Batch: 60/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.250
  Likelihood parameter acceptance fractions: [0.84]
Batch: 61/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.85]
Batch: 62/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.270
  Likelihood parameter acceptance fractions: [0.84]
Batch: 63/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.84]
Batch: 64/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.82]
Batch: 65/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.78]
Batch: 66/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.72]
Batch: 67/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.81]
Batch: 68/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.250
  Likelihood parameter acceptance fractions: [0.85]
Batch: 69/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.8]
Batch: 70/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.270
  Likelihood parameter acceptance fractions: [0.78]
Batch: 71/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.230
  Likelihood parameter acceptance fractions: [0.75]
Batch: 72/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.79]
Batch: 73/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.79]
Batch: 74/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.92]
Batch: 75/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.220
  Likelihood parameter acceptance fractions: [0.83]
Batch: 76/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.86]
Batch: 77/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.83]
Batch: 78/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.240
  Likelihood parameter acceptance fractions: [0.8]
Batch: 79/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.82]
Batch: 80/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.77]
Batch: 81/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.84]
Batch: 82/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.82]
Batch: 83/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.86]
Batch: 84/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.83]
Batch: 85/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.81]
Batch: 86/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.84]
Batch: 87/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.81]
Batch: 88/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.82]
Batch: 89/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.84]
Batch: 90/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
  Likelihood parameter acceptance fractions: [0.79]
Batch: 91/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.82]
Batch: 92/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.81]
Batch: 93/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.81]
Batch: 94/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.89]
Batch: 95/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.83]
Batch: 96/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.8]
Batch: 97/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.77]
Batch: 98/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.88]
Batch: 99/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.260
  Likelihood parameter acceptance fractions: [0.81]
Batch: 100/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.89]
CPU times: user 7.14 s, sys: 353 ms, total: 7.49 s
Wall time: 7.23 s

def plot_chains_with_err(walker, model, true_params):
    fig, axes = plt.subplots(
        walker.model_sampler.chain.shape[1] + 3, 1, figsize=(8, 8), sharex=True
    )
    for i in range(walker.model_sampler.chain.shape[1]):
        axes[i].plot(walker.model_sampler.chain[:, i])
        axes[i].set_ylabel(f"${model.params[i].latex_name}$ [{model.params[i].unit}]")
        true_value = true_params[model.params[i].name]
        axes[i].hlines(
            true_value, 0, len(walker.model_sampler.chain), "r", linestyle="--"
        )

    axes[-3].plot(walker.model_sampler.logp_chain)
    axes[-3].set_ylabel(r"$\log{\mathcal{L}(\alpha_i | \mathcal{O})}$")

    lmp = walker.likelihood_samplers[0].params[0]
    axes[-2].plot(walker.likelihood_samplers[0].chain)
    axes[-2].set_ylabel(f"${lmp.latex_name}$ [{lmp.unit}]")
    axes[-2].hlines(
        np.log(noise_fraction),
        0,
        len(walker.likelihood_samplers[0].chain),
        "r",
        linestyle="--",
    )

    axes[-1].plot(walker.likelihood_samplers[0].logp_chain)
    axes[-1].set_ylabel(r"$\log{\mathcal{L}(\alpha_i | \mathcal{O})}$")

    axes[-1].set_xlabel(r"$i$")
    # plt.legend(title="chains", ncol=3,)

plot_chains_with_err(walker2, my_model, true_params)

fig = corner.corner(
    np.hstack([walker2.model_sampler.chain, walker2.likelihood_samplers[0].chain]),
    labels=[p.name for p in my_model.params]
    + [walker2.likelihood_samplers[0].params[0].name],
    label="posterior",
    truths=[true_params["m"], true_params["b"], np.log(noise_fraction)],
)
fig.suptitle("posterior")

Text(0.5, 0.98, 'posterior')

plot_predictive_post(
    walker=walker2, model=my_model, x=x, y_exp=y_exp, y_err=y_stat_err, y_true=y_true
)
plt.title("option 2: unknown statistical error, systematic ignored")

Text(0.5, 1.0, 'option 2: unknown statistical error, systematic ignored')

Run option 3: correct formulation of the systematic error

walker3 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_sys_correct,
    rng=rng,
)

%%time
walker3.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.750
CPU times: user 3.45 s, sys: 22.2 ms, total: 3.47 s
Wall time: 4.39 s

plot_chains(walker=walker3, model=my_model, true_params=true_params)

plot_posterior_corner(walker=walker3, true_params=true_params)

plot_predictive_post(
    walker=walker3, model=my_model, x=x, y_exp=y_exp, y_err=y_stat_err, y_true=y_true
)
plt.title("option 3: systematic included correctly")

Text(0.5, 1.0, 'option 3: systematic included correctly')

Run option 4: incorrect formulation of the systematic error

walker4 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_sys_wrong,
    rng=rng,
)

%%time
walker4.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.733
CPU times: user 2.67 s, sys: 72.3 ms, total: 2.74 s
Wall time: 2.78 s

plot_chains(walker=walker4, model=my_model, true_params=true_params)

plot_posterior_corner(walker=walker4, true_params=true_params)

plot_predictive_post(
    walker=walker4, model=my_model, x=x, y_exp=y_exp, y_err=y_stat_err, y_true=y_true
)
plt.title("option 4: systematic included incorrectly")

Text(0.5, 1.0, 'option 4: systematic included incorrectly')

Multiple constraints

Let’s choose a second constraint, with the same normalization bias in the opposite direction and the same coverage over the $x$-domain.

systematic_fractional_err2 = 0.1
# choose a normalization 1 std deviation above the mean this time
N2 = 1 + systematic_fractional_err2
noise_fraction2 = 0.025
x2 = np.linspace(0.01, 0.8, 27, dtype=float)
y_true2 = my_model.y(x2, *list(true_params.values()))
y_exp2 = (y_true2 + rng.normal(scale=noise_fraction2 * y_true2, size=len(x2))) * N2
y_stat_err2 = noise_fraction2 * y_exp2 * N2

x_full = np.linspace(-1, 2, 100)
y_true_full = my_model.y(x_full, *list(true_params.values()))

plt.errorbar(
    x,
    y_exp,
    y_stat_err,
    marker="o",
    linestyle="none",
    label="experiment 1",
    color="tab:purple",
)
plt.errorbar(
    x2,
    y_exp2,
    y_stat_err2,
    marker="o",
    linestyle="none",
    label="experiment 2",
    color="tab:cyan",
)

plt.plot(x_full, y_true_full, "k--", label="truth")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.title("multiple experimental constraint with opposite bias")

Text(0.5, 1.0, 'multiple experimental constraint with opposite bias')

# 1 and 2
obs_stat_only2 = rxmc.observation.Observation(
    x=x2,
    y=y_exp2,
    y_stat_err=y_stat_err2,
)

# 3
obs_sys_norm_correct2 = rxmc.observation.Observation(
    x=x2,
    y=y_exp2,
    y_stat_err=y_stat_err2,
    y_sys_err_normalization=systematic_fractional_err2,
)

# 4
obs_sys_norm_wrong2 = rxmc.observation.FixedCovarianceObservation(
    x=x2,
    y=y_exp2,
    covariance=np.diag(y_stat_err2**2)
    + systematic_fractional_err2**2 * np.outer(y_exp2, y_exp2),
)

set up likelihood models and constraints

# 1
evidence_stat_only = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_stat_only, obs_stat_only2],
            my_model,
            likelihood,
        )
    ]
)

# 2
evidence_unknown_stat = rxmc.evidence.Evidence(
    constraints=[],
    parametric_constraints=[
        rxmc.constraint.Constraint(
            [obs_stat_only, obs_stat_only2],
            my_model,
            likelihood_unknown_stat,
        )
    ],
)

# 3
evidence_sys_correct = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_correct, obs_sys_norm_correct2],
            my_model,
            likelihood,
        )
    ]
)

# 4
evidence_sys_wrong = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_wrong, obs_sys_norm_wrong2],
            my_model,
            likelihood_fixed_cov,
        )
    ]
)

walker1 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_stat_only,
    rng=rng,
)

%%time
walker1.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.151
CPU times: user 4.27 s, sys: 148 ms, total: 4.42 s
Wall time: 4.79 s

upper1, med1, lower1 = np.percentile(
    [my_model.y(x_full, *p) for p in walker1.model_sampler.chain],
    [5, 50, 95],
    axis=0,
)

walker2 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_unknown_stat,
    likelihood_samplers=[
        rxmc.param_sampling.MetropolisHastingsSampler(
            params=likelihood_unknown_stat.params,
            starting_location=noise_prior.mean(),
            proposal=proposal_distribution_noise,
            prior=noise_prior,
        )
    ],
    rng=rng,
)

%%time
walker2.walk(
    n_steps=10000,
    burnin=1000,
    batch_size=100,
)

Burn-in batch 1/10 completed, 100 steps.
Burn-in batch 2/10 completed, 100 steps.
Burn-in batch 3/10 completed, 100 steps.
Burn-in batch 4/10 completed, 100 steps.
Burn-in batch 5/10 completed, 100 steps.
Burn-in batch 6/10 completed, 100 steps.
Burn-in batch 7/10 completed, 100 steps.
Burn-in batch 8/10 completed, 100 steps.
Burn-in batch 9/10 completed, 100 steps.
Burn-in batch 10/10 completed, 100 steps.
Batch: 1/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.72]
Batch: 2/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.73]
Batch: 3/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.67]
Batch: 4/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.77]
Batch: 5/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.76]
Batch: 6/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.68]
Batch: 7/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.550
  Likelihood parameter acceptance fractions: [0.68]
Batch: 8/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.77]
Batch: 9/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.7]
Batch: 10/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.77]
Batch: 11/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.67]
Batch: 12/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.73]
Batch: 13/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.74]
Batch: 14/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.82]
Batch: 15/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.77]
Batch: 16/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.82]
Batch: 17/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.75]
Batch: 18/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.78]
Batch: 19/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.64]
Batch: 20/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.8]
Batch: 21/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.71]
Batch: 22/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.69]
Batch: 23/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.68]
Batch: 24/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.64]
Batch: 25/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.72]
Batch: 26/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.73]
Batch: 27/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.78]
Batch: 28/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.69]
Batch: 29/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.72]
Batch: 30/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.78]
Batch: 31/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.78]
Batch: 32/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.65]
Batch: 33/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.580
  Likelihood parameter acceptance fractions: [0.75]
Batch: 34/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.560
  Likelihood parameter acceptance fractions: [0.76]
Batch: 35/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.72]
Batch: 36/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.73]
Batch: 37/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
  Likelihood parameter acceptance fractions: [0.75]
Batch: 38/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.73]
Batch: 39/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.73]
Batch: 40/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.69]
Batch: 41/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.75]
Batch: 42/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
  Likelihood parameter acceptance fractions: [0.63]
Batch: 43/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.75]
Batch: 44/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.540
  Likelihood parameter acceptance fractions: [0.73]
Batch: 45/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.79]
Batch: 46/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.73]
Batch: 47/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.77]
Batch: 48/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.71]
Batch: 49/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.65]
Batch: 50/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
  Likelihood parameter acceptance fractions: [0.72]
Batch: 51/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.74]
Batch: 52/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.74]
Batch: 53/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.78]
Batch: 54/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.73]
Batch: 55/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.75]
Batch: 56/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.71]
Batch: 57/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.77]
Batch: 58/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.7]
Batch: 59/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.78]
Batch: 60/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.73]
Batch: 61/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.6]
Batch: 62/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.71]
Batch: 63/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.69]
Batch: 64/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.74]
Batch: 65/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.72]
Batch: 66/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.69]
Batch: 67/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.76]
Batch: 68/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.72]
Batch: 69/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.77]
Batch: 70/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.74]
Batch: 71/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.71]
Batch: 72/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.82]
Batch: 73/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.73]
Batch: 74/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.82]
Batch: 75/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.75]
Batch: 76/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.540
  Likelihood parameter acceptance fractions: [0.77]
Batch: 77/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.76]
Batch: 78/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.68]
Batch: 79/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.79]
Batch: 80/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.72]
Batch: 81/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.65]
Batch: 82/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.69]
Batch: 83/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.7]
Batch: 84/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.65]
Batch: 85/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.71]
Batch: 86/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.74]
Batch: 87/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.69]
Batch: 88/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.71]
Batch: 89/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.7]
Batch: 90/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.75]
Batch: 91/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.75]
Batch: 92/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.74]
Batch: 93/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.71]
Batch: 94/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.75]
Batch: 95/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.67]
Batch: 96/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.65]
Batch: 97/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.68]
Batch: 98/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.74]
Batch: 99/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.72]
Batch: 100/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.76]
CPU times: user 8.81 s, sys: 15.8 ms, total: 8.82 s
Wall time: 8.82 s

upper2, med2, lower2 = np.percentile(
    [my_model.y(x_full, *p) for p in walker2.model_sampler.chain],
    [5, 50, 95],
    axis=0,
)

walker3 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_sys_correct,
    rng=rng,
)

%%time
walker3.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.590
CPU times: user 4.21 s, sys: 7.66 ms, total: 4.21 s
Wall time: 4.22 s

upper3, med3, lower3 = np.percentile(
    [my_model.y(x_full, *p) for p in walker3.model_sampler.chain],
    [5, 50, 95],
    axis=0,
)

walker4 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_sys_wrong,
    rng=rng,
)

%%time
walker4.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.578
CPU times: user 2.72 s, sys: 40.3 ms, total: 2.76 s
Wall time: 2.74 s

upper4, med4, lower4 = np.percentile(
    [my_model.y(x_full, *p) for p in walker4.model_sampler.chain],
    [5, 50, 95],
    axis=0,
)

plt.errorbar(
    x,
    y_exp,
    y_stat_err,
    marker=".",
    linestyle="none",
    label="experiment 1",
    color="tab:purple",
    zorder=999,
)
plt.errorbar(
    x2,
    y_exp2,
    y_stat_err2,
    marker=".",
    linestyle="none",
    label="experiment 2",
    color="tab:cyan",
    zorder=999,
)
p = plt.fill_between(x_full, lower1, upper1, label="stat only", alpha=0.5, zorder=99)
plt.plot(x_full, med1, "--", color=p.get_facecolor(), alpha=1, zorder=100)


p = plt.fill_between(x_full, lower2, upper2, label="stat fit", alpha=0.5, zorder=89)
plt.plot(x_full, med2, "--", color=p.get_facecolor(), alpha=1, zorder=90)


p = plt.fill_between(x_full, lower3, upper3, label="full covariance", alpha=0.5)
plt.plot(x_full, med3, "--", color=p.get_facecolor(), alpha=1, zorder=89)

p = plt.fill_between(x_full, lower4, upper4, label="full covariance, wrong", alpha=0.25)
plt.plot(x_full, med4, "--", color=p.get_facecolor(), alpha=0.5)


plt.plot(x_full, y_true_full, "k--", label="truth", zorder=999)
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.title("multiple constraints with systematic normalization error")

Text(0.5, 1.0, 'multiple constraints with systematic normalization error')

Multiple constraints with offset domain

Let’s choose a second constraint, with the same normalization bias in the opposite direction and slightly different coverage opver $x$.

x2 = np.linspace(0.6, 1.4, 27, dtype=float)
y_true2 = my_model.y(x2, *list(true_params.values()))
y_exp2 = (y_true2 + rng.normal(scale=noise_fraction2 * y_true2, size=len(x2))) * N2
y_stat_err2 = noise_fraction2 * y_exp2 * N2

plt.errorbar(
    x,
    y_exp,
    y_stat_err,
    marker="o",
    linestyle="none",
    label="experiment 1",
    color="tab:purple",
)
plt.errorbar(
    x2,
    y_exp2,
    y_stat_err2,
    marker="o",
    linestyle="none",
    label="experiment 2",
    color="tab:cyan",
)

plt.plot(x_full, y_true_full, "k--", label="truth")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.title("$x$-offset experimental constraint with opposite bias")

Text(0.5, 1.0, '$x$-offset experimental constraint with opposite bias')

# 1 and 2
obs_stat_only2 = rxmc.observation.Observation(
    x=x2,
    y=y_exp2,
    y_stat_err=y_stat_err2,
)

# 3
obs_sys_norm_correct2 = rxmc.observation.Observation(
    x=x2,
    y=y_exp2,
    y_stat_err=y_stat_err2,
    y_sys_err_normalization=systematic_fractional_err2,
)

# 4
obs_sys_norm_wrong2 = rxmc.observation.FixedCovarianceObservation(
    x=x2,
    y=y_exp2,
    covariance=np.diag(y_stat_err2**2)
    + systematic_fractional_err2**2 * np.outer(y_exp2, y_exp2),
)

set up likelihood models and constraints

# 1
evidence_stat_only = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_stat_only, obs_stat_only2],
            my_model,
            likelihood,
        )
    ]
)

# 2
evidence_unknown_stat = rxmc.evidence.Evidence(
    constraints=[],
    parametric_constraints=[
        rxmc.constraint.Constraint(
            [obs_stat_only, obs_stat_only2],
            my_model,
            likelihood_unknown_stat,
        )
    ],
)

# 3
evidence_sys_correct = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_correct, obs_sys_norm_correct2],
            my_model,
            likelihood,
        )
    ]
)

# 4
evidence_sys_wrong = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_wrong, obs_sys_norm_wrong2],
            my_model,
            likelihood_fixed_cov,
        )
    ]
)

walker1 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_stat_only,
    rng=rng,
)

%%time
walker1.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.144
CPU times: user 4.18 s, sys: 3.95 ms, total: 4.19 s
Wall time: 4.19 s

upper1, med1, lower1 = np.percentile(
    [my_model.y(x_full, *p) for p in walker1.model_sampler.chain],
    [5, 50, 95],
    axis=0,
)

walker2 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_unknown_stat,
    likelihood_samplers=[
        rxmc.param_sampling.MetropolisHastingsSampler(
            params=likelihood_unknown_stat.params,
            starting_location=noise_prior.mean(),
            proposal=proposal_distribution_noise,
            prior=noise_prior,
        )
    ],
    rng=rng,
)

%%time
walker2.walk(
    n_steps=10000,
    burnin=1000,
    batch_size=100,
)

Burn-in batch 1/10 completed, 100 steps.
Burn-in batch 2/10 completed, 100 steps.
Burn-in batch 3/10 completed, 100 steps.
Burn-in batch 4/10 completed, 100 steps.
Burn-in batch 5/10 completed, 100 steps.
Burn-in batch 6/10 completed, 100 steps.
Burn-in batch 7/10 completed, 100 steps.
Burn-in batch 8/10 completed, 100 steps.
Burn-in batch 9/10 completed, 100 steps.
Burn-in batch 10/10 completed, 100 steps.
Batch: 1/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.72]
Batch: 2/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.230
  Likelihood parameter acceptance fractions: [0.68]
Batch: 3/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.73]
Batch: 4/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.550
  Likelihood parameter acceptance fractions: [0.78]
Batch: 5/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.73]
Batch: 6/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.73]
Batch: 7/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.76]
Batch: 8/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.74]
Batch: 9/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.69]
Batch: 10/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.69]
Batch: 11/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.76]
Batch: 12/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.69]
Batch: 13/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.74]
Batch: 14/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.73]
Batch: 15/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.75]
Batch: 16/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.76]
Batch: 17/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.73]
Batch: 18/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.75]
Batch: 19/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.66]
Batch: 20/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.550
  Likelihood parameter acceptance fractions: [0.79]
Batch: 21/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.72]
Batch: 22/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.77]
Batch: 23/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.76]
Batch: 24/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.74]
Batch: 25/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.62]
Batch: 26/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.72]
Batch: 27/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.270
  Likelihood parameter acceptance fractions: [0.68]
Batch: 28/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.69]
Batch: 29/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.72]
Batch: 30/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.74]
Batch: 31/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.67]
Batch: 32/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.200
  Likelihood parameter acceptance fractions: [0.69]
Batch: 33/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.75]
Batch: 34/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.77]
Batch: 35/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.67]
Batch: 36/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.73]
Batch: 37/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.240
  Likelihood parameter acceptance fractions: [0.69]
Batch: 38/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.67]
Batch: 39/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.76]
Batch: 40/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.73]
Batch: 41/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.7]
Batch: 42/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.7]
Batch: 43/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.78]
Batch: 44/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.7]
Batch: 45/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.65]
Batch: 46/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.79]
Batch: 47/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.71]
Batch: 48/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.63]
Batch: 49/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.72]
Batch: 50/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.73]
Batch: 51/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.66]
Batch: 52/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.67]
Batch: 53/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.65]
Batch: 54/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.71]
Batch: 55/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.73]
Batch: 56/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.240
  Likelihood parameter acceptance fractions: [0.75]
Batch: 57/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.77]
Batch: 58/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.79]
Batch: 59/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.72]
Batch: 60/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.74]
Batch: 61/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.76]
Batch: 62/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.220
  Likelihood parameter acceptance fractions: [0.75]
Batch: 63/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.71]
Batch: 64/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.75]
Batch: 65/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.77]
Batch: 66/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.65]
Batch: 67/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.72]
Batch: 68/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.7]
Batch: 69/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.77]
Batch: 70/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.270
  Likelihood parameter acceptance fractions: [0.75]
Batch: 71/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.72]
Batch: 72/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.75]
Batch: 73/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.65]
Batch: 74/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.69]
Batch: 75/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.77]
Batch: 76/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.71]
Batch: 77/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.69]
Batch: 78/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.77]
Batch: 79/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.73]
Batch: 80/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.84]
Batch: 81/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.74]
Batch: 82/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.7]
Batch: 83/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.63]
Batch: 84/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.74]
Batch: 85/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.66]
Batch: 86/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.73]
Batch: 87/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.66]
Batch: 88/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.76]
Batch: 89/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.73]
Batch: 90/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.73]
Batch: 91/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.72]
Batch: 92/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.73]
Batch: 93/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.73]
Batch: 94/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.7]
Batch: 95/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.270
  Likelihood parameter acceptance fractions: [0.81]
Batch: 96/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.76]
Batch: 97/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.78]
Batch: 98/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.76]
Batch: 99/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.270
  Likelihood parameter acceptance fractions: [0.79]
Batch: 100/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.81]
CPU times: user 8.96 s, sys: 223 ms, total: 9.19 s
Wall time: 8.92 s

upper2, med2, lower2 = np.percentile(
    [my_model.y(x_full, *p) for p in walker2.model_sampler.chain],
    [5, 50, 95],
    axis=0,
)

walker3 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_sys_correct,
    rng=rng,
)

%%time
walker3.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.626
CPU times: user 4.17 s, sys: 3.99 ms, total: 4.17 s
Wall time: 4.17 s

upper3, med3, lower3 = np.percentile(
    [my_model.y(x_full, *p) for p in walker3.model_sampler.chain],
    [5, 50, 95],
    axis=0,
)

walker4 = rxmc.walker.Walker(
    model_sampler=rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence=evidence_sys_wrong,
    rng=rng,
)

%%time
walker4.walk(n_steps=10000, burnin=1000)

Burn-in batch 1/1 completed, 1000 steps.
Batch: 1/1 completed, 10000 steps. 
  Model parameter acceptance fraction: 0.618
CPU times: user 2.77 s, sys: 113 ms, total: 2.88 s
Wall time: 2.78 s

upper4, med4, lower4 = np.percentile(
    [my_model.y(x_full, *p) for p in walker4.model_sampler.chain],
    [5, 50, 95],
    axis=0,
)

plt.errorbar(
    x,
    y_exp,
    y_stat_err,
    marker=".",
    linestyle="none",
    label="experiment 1",
    color="tab:purple",
    zorder=999,
)
plt.errorbar(
    x2,
    y_exp2,
    y_stat_err2,
    marker=".",
    linestyle="none",
    label="experiment 2",
    color="tab:cyan",
    zorder=999,
)
p = plt.fill_between(x_full, lower1, upper1, label="stat only", alpha=0.5, zorder=99)
plt.plot(x_full, med1, "--", color=p.get_facecolor(), alpha=1, zorder=100)


p = plt.fill_between(x_full, lower2, upper2, label="stat fit", alpha=0.5, zorder=89)
plt.plot(x_full, med2, "--", color=p.get_facecolor(), alpha=1, zorder=90)


p = plt.fill_between(x_full, lower3, upper3, label="full covariance", alpha=0.5)
plt.plot(x_full, med3, "--", color=p.get_facecolor(), alpha=1, zorder=89)

p = plt.fill_between(x_full, lower4, upper4, label="full covariance, wrong", alpha=0.25)
plt.plot(x_full, med4, "--", color=p.get_facecolor(), alpha=0.5)


plt.plot(x_full, y_true_full, "k--", label="truth", zorder=999)
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.title("$x$-offset constraints with systematic normalization error")
plt.savefig("systematic_err.pdf")