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

Using database version X4-2024-12-31 located in: /home/kyle/db/exfor/unpack_exfor-2024/X4-2024-12-31

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}

5. Unknown constant statistical noise, inferred alongside the model (option 2b):

   \begin{equation} \Sigma_{ij} = \delta_{ij} \epsilon_0^2 \end{equation}

6. An additive offset systematic, included as a rank-one mode with a constant basis (option 5):

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

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

# 2 (unknown statistical error: reported errors ignored, inferred instead)
obs_unknown_stat = rxmc.observation.Observation(x=x, y=y_exp)

# 3 (normalization systematic supplied as a covariance Term, scales with ym)
obs_sys_norm_correct = rxmc.observation.Observation(
    x=x,
    y=y_exp,
    y_stat_err=y_stat_err,
)

# 4 (the "wrong" fixed covariance, built from the DATA -> Peelle's Pertinent Puzzle)
obs_sys_norm_wrong = rxmc.observation.Observation(x=x, y=y_exp)
wrong_cov = 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 - the default Gaussian likelihood
likelihood = rxmc.likelihood_model.GaussianLikelihood()

# 2 - free noise-fraction parameter, sampled as a covariance nuisance
log_noise_fraction = rxmc.params.Parameter(
    "log noise fraction", float, latex_name=r"\log{\epsilon}", unit="dimensionless"
)

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

# 2
evidence_unknown_stat = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_unknown_stat],
            my_model,
            extra_terms=[rxmc.covariance.noise_fraction_term(log_noise_fraction)],
        )
    ]
)

# 3
evidence_sys_correct = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_correct],
            my_model,
            extra_terms=[
                rxmc.covariance.normalization_term(
                    magnitude=systematic_fractional_err,
                )
            ],
        )
    ]
)

# 4
evidence_sys_wrong = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_wrong],
            my_model,
            extra_terms=[rxmc.covariance.Term(wrong_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.72 s, sys: 25.3 ms, total: 3.75 s
Wall time: 3.72 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=[log_noise_fraction],
            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 8.51 s, sys: 70.6 ms, total: 8.58 s
Wall time: 8.49 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 2b: unknown constant noise

Option 2 inferred a fractional noise ($\Sigma_{ij} = \delta_{ij},\epsilon^2 y_m(x_j;\alpha)^2$, via noise_fraction_term). Its sibling noise_term infers a constant noise floor, \begin{equation} \Sigma_{ij} = \delta_{ij}, \epsilon_0^2 . \end{equation} Both are additive on top of any reported statistical diagonal — so to let the inferred noise replace the reported statistics, build the Observation without y_stat_err (as obs_unknown_stat is built), or pass include_statistical_term=False to the Constraint.

log_noise = rxmc.params.Parameter(
    "log noise", float, latex_name=r"\log{\epsilon_0}", unit="dimensionless"
)
evidence_unknown_const = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_unknown_stat],
            my_model,
            extra_terms=[rxmc.covariance.noise_term(log_noise)],
        )
    ]
)
# the true absolute noise scale is noise_fraction * y; center the prior there
const_noise_prior = stats.norm(loc=np.log(noise_fraction * np.mean(y_exp)), scale=1)

walker2b = 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_const,
    likelihood_samplers=[
        rxmc.param_sampling.MetropolisHastingsSampler(
            params=[log_noise],
            starting_location=const_noise_prior.mean(),
            proposal=proposal_distribution_noise,
            prior=const_noise_prior,
        )
    ],
    rng=rng,
)

%%time
walker2b.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.300
  Likelihood parameter acceptance fractions: [0.72]
Batch: 2/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.86]
Batch: 3/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.240
  Likelihood parameter acceptance fractions: [0.83]

Batch: 4/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.240
  Likelihood parameter acceptance fractions: [0.87]
Batch: 5/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.87]
Batch: 6/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.85]
Batch: 7/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.84]
Batch: 8/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.79]

Batch: 9/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.89]
Batch: 10/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.82]
Batch: 11/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.260
  Likelihood parameter acceptance fractions: [0.85]
Batch: 12/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.230
  Likelihood parameter acceptance fractions: [0.86]

Batch: 13/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.83]
Batch: 14/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.77]
Batch: 15/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.240
  Likelihood parameter acceptance fractions: [0.84]
Batch: 16/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.81]

Batch: 17/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.83]
Batch: 18/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.84]
Batch: 19/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.82]

Batch: 20/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.83]
Batch: 21/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.85]
Batch: 22/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.88]
Batch: 23/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.86]
Batch: 24/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.84]

Batch: 25/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.250
  Likelihood parameter acceptance fractions: [0.77]
Batch: 26/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.85]
Batch: 27/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.89]
Batch: 28/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.83]
Batch: 29/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.85]

Batch: 30/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.85]
Batch: 31/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.81]
Batch: 32/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.85]
Batch: 33/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.87]

Batch: 34/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.81]
Batch: 35/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.86]
Batch: 36/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.87]
Batch: 37/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.81]

Batch: 38/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.79]
Batch: 39/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.240
  Likelihood parameter acceptance fractions: [0.79]
Batch: 40/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.86]
Batch: 41/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.88]

Batch: 42/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.79]
Batch: 43/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
  Likelihood parameter acceptance fractions: [0.81]
Batch: 44/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.84]
Batch: 45/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.81]

Batch: 46/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.79]
Batch: 47/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.79]
Batch: 48/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.74]
Batch: 49/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.86]

Batch: 50/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.87]
Batch: 51/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.78]
Batch: 52/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.91]

Batch: 53/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.91]
Batch: 54/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.9]
Batch: 55/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.81]
Batch: 56/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.85]

Batch: 57/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.79]
Batch: 58/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.79]
Batch: 59/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.82]

Batch: 60/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.79]
Batch: 61/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.81]
Batch: 62/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.8]

Batch: 63/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.310
  Likelihood parameter acceptance fractions: [0.82]
Batch: 64/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.84]
Batch: 65/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.89]
Batch: 66/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.89]

Batch: 67/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.84]
Batch: 68/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.240
  Likelihood parameter acceptance fractions: [0.76]
Batch: 69/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.84]
Batch: 70/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.85]

Batch: 71/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.79]
Batch: 72/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.9]
Batch: 73/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.85]
Batch: 74/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.88]

Batch: 75/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.91]
Batch: 76/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.89]
Batch: 77/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.210
  Likelihood parameter acceptance fractions: [0.87]

Batch: 78/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.570
  Likelihood parameter acceptance fractions: [0.85]
Batch: 79/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.78]
Batch: 80/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.81]
Batch: 81/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.180
  Likelihood parameter acceptance fractions: [0.79]

Batch: 82/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.89]
Batch: 83/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.270
  Likelihood parameter acceptance fractions: [0.83]
Batch: 84/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.250
  Likelihood parameter acceptance fractions: [0.78]
Batch: 85/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.82]

Batch: 86/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.86]
Batch: 87/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.87]
Batch: 88/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.88]
Batch: 89/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.77]

Batch: 90/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.230
  Likelihood parameter acceptance fractions: [0.77]
Batch: 91/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.77]
Batch: 92/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.82]

Batch: 93/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.85]
Batch: 94/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.300
  Likelihood parameter acceptance fractions: [0.8]
Batch: 95/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
  Likelihood parameter acceptance fractions: [0.85]
Batch: 96/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.8]

Batch: 97/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.83]
Batch: 98/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.75]
Batch: 99/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.79]

Batch: 100/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.89]
CPU times: user 8.08 s, sys: 98.5 ms, total: 8.18 s
Wall time: 8.07 s

plot_chains(walker=walker2b, model=my_model, true_params=true_params)
plt.figure()
plt.plot(walker2b.likelihood_samplers[0].chain)
plt.axhline(np.log(noise_fraction * np.mean(y_exp)), color="r", ls="--")
plt.ylabel(r"$\log{\epsilon_0}$")
plt.xlabel("$i$");

plot_predictive_post(
    walker=walker2b, model=my_model, x=x, y_exp=y_exp, y_err=y_stat_err, y_true=y_true
)
plt.title("option 2b: unknown constant noise, systematic ignored")

Text(0.5, 1.0, 'option 2b: unknown constant noise, 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.752
CPU times: user 4.25 s, sys: 7.86 ms, total: 4.25 s
Wall time: 4.25 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.735
CPU times: user 3.37 s, sys: 40.1 ms, total: 3.41 s
Wall time: 3.37 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')

Run option 5: an additive offset systematic

A different systematic: the whole dataset is shifted by one unknown additive offset (think background mis-subtraction), rather than rescaled. The reported offset uncertainty $\omega$ enters as a rank-one mode with a constant basis, \begin{equation} \Sigma_{ij} = \delta_{ij},\sigma_{stat,i}^2 + \omega^2 , \end{equation} via offset_term(magnitude=...). With parameter= instead of magnitude=, the magnitude becomes a free nuisance ($\omega = e^{\theta}$), exactly parallel to options 2 and 2b. As with the normalization bias above, we shift the data by one standard deviation of the reported offset error, and compare accounting for it against ignoring it.

offset_syst_err = 0.3
# one common additive shift, chosen 1 std deviation of the offset error above 0
delta = offset_syst_err
y_exp_off = y_true + rng.normal(scale=noise_fraction * y_true, size=len(x)) + delta
y_stat_err_off = noise_fraction * y_exp_off

obs_offset = rxmc.observation.Observation(x=x, y=y_exp_off, y_stat_err=y_stat_err_off)
obs_offset_ignored = rxmc.observation.Observation(
    x=x, y=y_exp_off, y_stat_err=y_stat_err_off
)

evidence_offset = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_offset],
            my_model,
            extra_terms=[rxmc.covariance.offset_term(magnitude=offset_syst_err)],
        )
    ]
)
evidence_offset_ignored = rxmc.evidence.Evidence(
    [rxmc.constraint.Constraint([obs_offset_ignored], my_model, likelihood)]
)

# the free-nuisance spelling of the same mode (constructed, not fit here)
log_omega = rxmc.params.Parameter("log omega", float, latex_name=r"\log{\omega}")
free_offset = rxmc.constraint.Constraint(
    [rxmc.observation.Observation(x=x, y=y_exp_off, y_stat_err=y_stat_err_off)],
    my_model,
    extra_terms=[rxmc.covariance.offset_term(parameter=log_omega)],
)
print(f"free-offset constraint has {free_offset.n_params} nuisance parameter(s)")

free-offset constraint has 1 nuisance parameter(s)

%%time
walker5 = rxmc.walker.Walker(
    rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence_offset,
    rng=rng,
)
walker5i = rxmc.walker.Walker(
    rxmc.param_sampling.MetropolisHastingsSampler(
        params=my_model.params,
        starting_location=prior_distribution.mean,
        proposal=proposal_distribution_model,
        prior=prior_distribution,
    ),
    evidence_offset_ignored,
    rng=rng,
)
walker5.walk(n_steps=10000, burnin=1000, batch_size=100)
walker5i.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.820
Batch: 2/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.760
Batch: 3/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740
Batch: 4/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.830

Batch: 5/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 6/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780
Batch: 7/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.820
Batch: 8/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.720
Batch: 9/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.830
Batch: 10/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.830
Batch: 11/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800

Batch: 12/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 13/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.730
Batch: 14/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.790
Batch: 15/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.760
Batch: 16/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.830
Batch: 17/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840
Batch: 18/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800
Batch: 19/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800
Batch: 20/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.690

Batch: 21/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 22/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.880
Batch: 23/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.830
Batch: 24/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740
Batch: 25/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780
Batch: 26/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.860
Batch: 27/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.770

Batch: 28/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.760
Batch: 29/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 30/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 31/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840
Batch: 32/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840

Batch: 33/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.860
Batch: 34/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800
Batch: 35/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.790
Batch: 36/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780

Batch: 37/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.820
Batch: 38/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.770
Batch: 39/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740
Batch: 40/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.730
Batch: 41/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780
Batch: 42/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780
Batch: 43/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740

Batch: 44/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800
Batch: 45/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840

Batch: 46/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740
Batch: 47/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.830
Batch: 48/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780
Batch: 49/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780
Batch: 50/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 51/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.710
Batch: 52/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840

Batch: 53/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.830
Batch: 54/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.730

Batch: 55/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.730
Batch: 56/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740
Batch: 57/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.860
Batch: 58/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.830
Batch: 59/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 60/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840
Batch: 61/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840

Batch: 62/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.820
Batch: 63/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.790
Batch: 64/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780

Batch: 65/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.850
Batch: 66/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800
Batch: 67/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.760
Batch: 68/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840
Batch: 69/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800

Batch: 70/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.720
Batch: 71/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.770
Batch: 72/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780
Batch: 73/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 74/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780

Batch: 75/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.730
Batch: 76/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.820
Batch: 77/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740

Batch: 78/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 79/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740
Batch: 80/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.880
Batch: 81/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.710
Batch: 82/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.850
Batch: 83/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.810
Batch: 84/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840

Batch: 85/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.740
Batch: 86/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.790

Batch: 87/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.860
Batch: 88/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800
Batch: 89/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.760
Batch: 90/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.790
Batch: 91/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.750
Batch: 92/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.840
Batch: 93/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.780

Batch: 94/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.820
Batch: 95/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.820

Batch: 96/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.870
Batch: 97/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.820
Batch: 98/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.750
Batch: 99/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.770
Batch: 100/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.800
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.480
Batch: 2/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
Batch: 3/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 4/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
Batch: 5/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
Batch: 6/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440

Batch: 7/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
Batch: 8/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
Batch: 9/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
Batch: 10/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
Batch: 11/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
Batch: 12/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
Batch: 13/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
Batch: 14/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 15/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 16/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
Batch: 17/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.540
Batch: 18/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480

Batch: 19/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
Batch: 20/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
Batch: 21/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 22/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
Batch: 23/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 24/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
Batch: 25/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
Batch: 26/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430

Batch: 27/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440

Batch: 28/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
Batch: 29/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
Batch: 30/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
Batch: 31/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
Batch: 32/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
Batch: 33/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
Batch: 34/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410

Batch: 35/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
Batch: 36/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 37/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490

Batch: 38/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
Batch: 39/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
Batch: 40/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 41/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
Batch: 42/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
Batch: 43/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
Batch: 44/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460

Batch: 45/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
Batch: 46/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
Batch: 47/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.540

Batch: 48/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
Batch: 49/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.330
Batch: 50/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
Batch: 51/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
Batch: 52/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
Batch: 53/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
Batch: 54/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440

Batch: 55/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380

Batch: 56/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
Batch: 57/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
Batch: 58/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
Batch: 59/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
Batch: 60/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
Batch: 61/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
Batch: 62/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430

Batch: 63/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
Batch: 64/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480

Batch: 65/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
Batch: 66/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
Batch: 67/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 68/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
Batch: 69/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
Batch: 70/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410

Batch: 71/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
Batch: 72/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
Batch: 73/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500

Batch: 74/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
Batch: 75/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
Batch: 76/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
Batch: 77/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
Batch: 78/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490

Batch: 79/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
Batch: 80/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
Batch: 81/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
Batch: 82/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480

Batch: 83/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
Batch: 84/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
Batch: 85/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
Batch: 86/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
Batch: 87/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420

Batch: 88/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
Batch: 89/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
Batch: 90/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470

Batch: 91/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
Batch: 92/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
Batch: 93/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
Batch: 94/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
Batch: 95/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
Batch: 96/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430

Batch: 97/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
Batch: 98/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440

Batch: 99/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
Batch: 100/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
CPU times: user 6.78 s, sys: 83.6 ms, total: 6.87 s
Wall time: 6.76 s

for name, w in [("offset accounted", walker5), ("offset ignored", walker5i)]:
    ch = w.model_sampler.chain
    print(
        f"{name:18s} m = {ch[:, 0].mean():.3f} ± {ch[:, 0].std():.3f}   "
        f"b = {ch[:, 1].mean():.3f} ± {ch[:, 1].std():.3f}"
    )

fig = corner.corner(
    walker5.model_sampler.chain,
    labels=["m", "b"],
    truths=[true_params["m"], true_params["b"]],
    truth_color="k",
    color="tab:blue",
)
corner.corner(walker5i.model_sampler.chain, fig=fig, color="tab:red")
plt.plot([], [], color="tab:blue", label="offset accounted")
plt.plot([], [], color="tab:red", label="offset ignored")
fig.legend(loc="upper right");

offset accounted   m = 0.473 ± 0.108   b = 2.357 ± 0.299
offset ignored     m = 0.472 ± 0.099   b = 2.332 ± 0.056

plot_predictive_post(
    walker=walker5,
    model=my_model,
    x=x,
    y_exp=y_exp_off,
    y_err=y_stat_err_off,
    y_true=y_true,
)
plt.title("option 5: additive offset accounted via offset_term")

Text(0.5, 1.0, 'option 5: additive offset accounted via offset_term')

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)
obs_unknown_stat2 = rxmc.observation.Observation(x=x2, y=y_exp2)

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

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

set up likelihood models and constraints

N1 = obs_stat_only.n_data_pts
N2 = obs_stat_only2.n_data_pts
s1 = np.arange(N1)
s2 = np.arange(N1, N1 + N2)

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

# 2 (shared noise-fraction parameter across both datasets -> case B)
evidence_unknown_stat = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_unknown_stat, obs_unknown_stat2],
            my_model,
            extra_terms=[
                rxmc.covariance.noise_fraction_term(log_noise_fraction, support=s1),
                rxmc.covariance.noise_fraction_term(log_noise_fraction, support=s2),
            ],
        )
    ]
)

# 3
evidence_sys_correct = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_correct, obs_sys_norm_correct2],
            my_model,
            extra_terms=[
                rxmc.covariance.normalization_term(
                    magnitude=systematic_fractional_err,
                    support=s1,
                ),
                rxmc.covariance.normalization_term(
                    magnitude=systematic_fractional_err2,
                    support=s2,
                ),
            ],
        )
    ]
)

# 4
evidence_sys_wrong = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_wrong, obs_sys_norm_wrong2],
            my_model,
            extra_terms=[
                rxmc.covariance.Term(wrong_cov, support=s1),
                rxmc.covariance.Term(wrong_cov2, support=s2),
            ],
        )
    ]
)

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 3.72 s, sys: 18.9 ms, total: 3.74 s
Wall time: 3.72 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=[log_noise_fraction],
            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.530
  Likelihood parameter acceptance fractions: [0.67]
Batch: 2/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.75]

Batch: 3/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.75]

Batch: 4/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.67]
Batch: 5/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.81]

Batch: 6/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.69]

Batch: 7/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.75]
Batch: 8/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.69]

Batch: 9/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.75]

Batch: 10/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.66]
Batch: 11/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.74]

Batch: 12/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.62]

Batch: 13/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.77]
Batch: 14/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.73]

Batch: 15/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.78]

Batch: 16/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.73]
Batch: 17/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.77]

Batch: 18/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.77]

Batch: 19/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.7]
Batch: 20/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
  Likelihood parameter acceptance fractions: [0.77]

Batch: 21/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.520
  Likelihood parameter acceptance fractions: [0.75]

Batch: 22/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.71]
Batch: 23/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.74]

Batch: 24/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.280
  Likelihood parameter acceptance fractions: [0.83]

Batch: 25/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.78]
Batch: 26/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.71]

Batch: 27/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.570
  Likelihood parameter acceptance fractions: [0.73]

Batch: 28/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.69]
Batch: 29/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.69]

Batch: 30/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.230
  Likelihood parameter acceptance fractions: [0.79]

Batch: 31/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.81]
Batch: 32/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.560
  Likelihood parameter acceptance fractions: [0.75]

Batch: 33/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.74]

Batch: 34/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.560
  Likelihood parameter acceptance fractions: [0.71]
Batch: 35/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.320
  Likelihood parameter acceptance fractions: [0.74]

Batch: 36/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.65]

Batch: 37/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.69]
Batch: 38/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.620
  Likelihood parameter acceptance fractions: [0.8]

Batch: 39/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.570
  Likelihood parameter acceptance fractions: [0.76]

Batch: 40/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.490
  Likelihood parameter acceptance fractions: [0.64]
Batch: 41/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.75]

Batch: 42/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.74]

Batch: 43/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.81]
Batch: 44/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.75]

Batch: 45/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.8]

Batch: 46/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.7]
Batch: 47/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.74]

Batch: 48/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.78]

Batch: 49/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.73]
Batch: 50/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.77]

Batch: 51/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.480
  Likelihood parameter acceptance fractions: [0.72]

Batch: 52/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.77]
Batch: 53/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.450
  Likelihood parameter acceptance fractions: [0.78]

Batch: 54/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.550
  Likelihood parameter acceptance fractions: [0.75]

Batch: 55/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.7]
Batch: 56/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.71]

Batch: 57/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.73]

Batch: 58/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.77]
Batch: 59/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.65]

Batch: 60/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.75]

Batch: 61/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.340
  Likelihood parameter acceptance fractions: [0.71]
Batch: 62/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.75]

Batch: 63/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.290
  Likelihood parameter acceptance fractions: [0.68]

Batch: 64/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.71]
Batch: 65/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.67]

Batch: 66/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.75]

Batch: 67/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.72]
Batch: 68/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.560
  Likelihood parameter acceptance fractions: [0.69]

Batch: 69/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.72]

Batch: 70/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.540
  Likelihood parameter acceptance fractions: [0.8]
Batch: 71/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.68]

Batch: 72/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.360
  Likelihood parameter acceptance fractions: [0.77]

Batch: 73/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.74]
Batch: 74/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.74]

Batch: 75/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.71]

Batch: 76/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.72]
Batch: 77/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.410
  Likelihood parameter acceptance fractions: [0.67]

Batch: 78/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.420
  Likelihood parameter acceptance fractions: [0.67]

Batch: 79/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.440
  Likelihood parameter acceptance fractions: [0.69]
Batch: 80/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.77]

Batch: 81/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.72]

Batch: 82/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.470
  Likelihood parameter acceptance fractions: [0.78]
Batch: 83/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.530
  Likelihood parameter acceptance fractions: [0.76]

Batch: 84/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.460
  Likelihood parameter acceptance fractions: [0.79]

Batch: 85/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.580
  Likelihood parameter acceptance fractions: [0.66]
Batch: 86/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.74]

Batch: 87/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.370
  Likelihood parameter acceptance fractions: [0.68]

Batch: 88/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.73]
Batch: 89/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.75]

Batch: 90/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.540
  Likelihood parameter acceptance fractions: [0.72]

Batch: 91/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.500
  Likelihood parameter acceptance fractions: [0.72]
Batch: 92/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.510
  Likelihood parameter acceptance fractions: [0.71]

Batch: 93/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.390
  Likelihood parameter acceptance fractions: [0.75]

Batch: 94/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.400
  Likelihood parameter acceptance fractions: [0.73]
Batch: 95/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.380
  Likelihood parameter acceptance fractions: [0.76]

Batch: 96/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.550
  Likelihood parameter acceptance fractions: [0.72]

Batch: 97/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.430
  Likelihood parameter acceptance fractions: [0.72]
Batch: 98/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.350
  Likelihood parameter acceptance fractions: [0.61]

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.610
  Likelihood parameter acceptance fractions: [0.73]
CPU times: user 9.78 s, sys: 10.8 ms, total: 9.79 s
Wall time: 9.77 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.592
CPU times: user 5.4 s, sys: 841 μs, total: 5.4 s
Wall time: 5.4 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 3.63 s, sys: 30.1 ms, total: 3.66 s
Wall time: 3.63 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)
obs_unknown_stat2 = rxmc.observation.Observation(x=x2, y=y_exp2)

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

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

set up likelihood models and constraints

N1 = obs_stat_only.n_data_pts
N2 = obs_stat_only2.n_data_pts
s1 = np.arange(N1)
s2 = np.arange(N1, N1 + N2)

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

# 2 (shared noise-fraction parameter across both datasets -> case B)
evidence_unknown_stat = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_unknown_stat, obs_unknown_stat2],
            my_model,
            extra_terms=[
                rxmc.covariance.noise_fraction_term(log_noise_fraction, support=s1),
                rxmc.covariance.noise_fraction_term(log_noise_fraction, support=s2),
            ],
        )
    ]
)

# 3
evidence_sys_correct = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_correct, obs_sys_norm_correct2],
            my_model,
            extra_terms=[
                rxmc.covariance.normalization_term(
                    magnitude=systematic_fractional_err,
                    support=s1,
                ),
                rxmc.covariance.normalization_term(
                    magnitude=systematic_fractional_err2,
                    support=s2,
                ),
            ],
        )
    ]
)

# 4
evidence_sys_wrong = rxmc.evidence.Evidence(
    [
        rxmc.constraint.Constraint(
            [obs_sys_norm_wrong, obs_sys_norm_wrong2],
            my_model,
            extra_terms=[
                rxmc.covariance.Term(wrong_cov, support=s1),
                rxmc.covariance.Term(wrong_cov2, support=s2),
            ],
        )
    ]
)

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.331
CPU times: user 3.58 s, sys: 16 ms, total: 3.6 s
Wall time: 3.58 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=[log_noise_fraction],
            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: 1.000
  Likelihood parameter acceptance fractions: [0.83]
Batch: 2/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.940
  Likelihood parameter acceptance fractions: [0.76]

Batch: 3/100 completed, 100 steps. 
  Model parameter acceptance fraction: 1.000
  Likelihood parameter acceptance fractions: [0.71]

Batch: 4/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.69]
Batch: 5/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.79]

Batch: 6/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.8]

Batch: 7/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.75]
Batch: 8/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.73]

Batch: 9/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.74]

Batch: 10/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.910
  Likelihood parameter acceptance fractions: [0.7]
Batch: 11/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.73]

Batch: 12/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.72]

Batch: 13/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.67]
Batch: 14/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.71]

Batch: 15/100 completed, 100 steps. 
  Model parameter acceptance fraction: 1.000
  Likelihood parameter acceptance fractions: [0.72]

Batch: 16/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.66]
Batch: 17/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.68]

Batch: 18/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.930
  Likelihood parameter acceptance fractions: [0.71]

Batch: 19/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.71]
Batch: 20/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.73]

Batch: 21/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.75]

Batch: 22/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.75]
Batch: 23/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.940
  Likelihood parameter acceptance fractions: [0.67]

Batch: 24/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.940
  Likelihood parameter acceptance fractions: [0.74]

Batch: 25/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.7]
Batch: 26/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.71]

Batch: 27/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.75]

Batch: 28/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.76]
Batch: 29/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.67]

Batch: 30/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.76]

Batch: 31/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.69]
Batch: 32/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.67]

Batch: 33/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.72]

Batch: 34/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.74]
Batch: 35/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.66]

Batch: 36/100 completed, 100 steps. 
  Model parameter acceptance fraction: 1.000
  Likelihood parameter acceptance fractions: [0.76]

Batch: 37/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.66]
Batch: 38/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.72]

Batch: 39/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.68]

Batch: 40/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.6]
Batch: 41/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.920
  Likelihood parameter acceptance fractions: [0.7]

Batch: 42/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.910
  Likelihood parameter acceptance fractions: [0.75]

Batch: 43/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.72]
Batch: 44/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.7]

Batch: 45/100 completed, 100 steps. 
  Model parameter acceptance fraction: 1.000
  Likelihood parameter acceptance fractions: [0.74]

Batch: 46/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.71]
Batch: 47/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.76]

Batch: 48/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.75]

Batch: 49/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.920
  Likelihood parameter acceptance fractions: [0.73]
Batch: 50/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.69]

Batch: 51/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.75]

Batch: 52/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.71]
Batch: 53/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.67]

Batch: 54/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.74]

Batch: 55/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.72]
Batch: 56/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.7]

Batch: 57/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.66]

Batch: 58/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.940
  Likelihood parameter acceptance fractions: [0.69]
Batch: 59/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.69]

Batch: 60/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.75]

Batch: 61/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.7]
Batch: 62/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.76]

Batch: 63/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.67]

Batch: 64/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.75]
Batch: 65/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.68]

Batch: 66/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.75]

Batch: 67/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.73]
Batch: 68/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.7]

Batch: 69/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.71]

Batch: 70/100 completed, 100 steps. 
  Model parameter acceptance fraction: 1.000
  Likelihood parameter acceptance fractions: [0.7]
Batch: 71/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.65]

Batch: 72/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.66]

Batch: 73/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.71]
Batch: 74/100 completed, 100 steps. 
  Model parameter acceptance fraction: 1.000
  Likelihood parameter acceptance fractions: [0.73]

Batch: 75/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.7]

Batch: 76/100 completed, 100 steps. 
  Model parameter acceptance fraction: 1.000
  Likelihood parameter acceptance fractions: [0.69]
Batch: 77/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.68]

Batch: 78/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.72]

Batch: 79/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.78]
Batch: 80/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.66]

Batch: 81/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.65]

Batch: 82/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.82]
Batch: 83/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.950
  Likelihood parameter acceptance fractions: [0.69]

Batch: 84/100 completed, 100 steps. 
  Model parameter acceptance fraction: 1.000
  Likelihood parameter acceptance fractions: [0.57]

Batch: 85/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.71]
Batch: 86/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.73]

Batch: 87/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.74]

Batch: 88/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.76]
Batch: 89/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.77]

Batch: 90/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.72]

Batch: 91/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.69]
Batch: 92/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.940
  Likelihood parameter acceptance fractions: [0.77]

Batch: 93/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.72]

Batch: 94/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.71]
Batch: 95/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.74]

Batch: 96/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.990
  Likelihood parameter acceptance fractions: [0.73]

Batch: 97/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.980
  Likelihood parameter acceptance fractions: [0.73]
Batch: 98/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.72]

Batch: 99/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.970
  Likelihood parameter acceptance fractions: [0.77]

Batch: 100/100 completed, 100 steps. 
  Model parameter acceptance fraction: 0.960
  Likelihood parameter acceptance fractions: [0.76]
CPU times: user 9.66 s, sys: 10.9 ms, total: 9.67 s
Wall time: 9.65 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.748
CPU times: user 5.42 s, sys: 897 μs, total: 5.42 s
Wall time: 5.42 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.735
CPU times: user 3.65 s, sys: 10.9 ms, total: 3.66 s
Wall time: 3.65 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")

Text(0.5, 1.0, '$x$-offset constraints with systematic normalization error')