Updating Beliefs with Observations

In Parts One and Two we set up the conceptual components of an AI scientist, defining the World with which it interacts, and the notion of Observations, experiment Designs and how experimental designs would be chosen with an eye to optimal information gain and model distinguishability. We also looked at how we might score a Hypothesis in a pool using the Bayes Factor and other metrics.

To integrate information from actual data into understanding which explanation is correct, we are always doing some form of inference. That can be entirely implicit (as humans we do that all the time), or we can make it explicit as inferring posterior distributions over explanations. However, so far we have assumed that we can exactly compute these values, and the resulting evidence for different models. Our toy model allowed us to do exact inference.

What happens when we can no longer do that? We may not be in a situation where our models can neatly provide log likelihoods and posteriors are non-trivial to obtain. This post addresses how the two main papers of interest(*reference to footnote with the references) address that problem, using the machinery we have built up so far. In addition, the code is written to allow for an LLM-driven model to do this implicitly.

This post focuses on parameter inference within each assumed model. Part 3b then compares the models and tests how inference affects predictions and experiment choices. We keep the evidence as a posterior normalizing constant here, but move model scoring and selection to that companion.

We'll keep using our simple LineWorld example and our pool of hypotheses -- no exciting Worlds or LLM-driven Proposers yet -- that's coming soon!

Summary of contents: how to infer parameters for a fixed model#

How do we make approximate methods recover each model's parameter posterior?

1. Define the observed data, likelihood, prior and parameter posterior.
2. Inspect prior importance sampling and tempered SMC.
3. Train a neural posterior estimator using labelled simulations.
4. Check parameter means, covariances and interval coverage against an exact reference.
5. Continue to Part 3b: comparing and selecting models.

The constant, line and quadratic hypotheses are treated separately here; we do not rank them. Part 3b contains both evidence-based and neural model comparison.

1. What is being inferred, and from which data?#

Fix a candidate model $h$ and let $\theta_h$ be its coefficient vector. For a line, these are its intercept and slope. We infer those coefficients conditional on this model, without deciding whether the model itself is the right explanation.

• $\mathcal D_0$ contains four measured outputs together with their known settings.
• $a_i$ and $o_i$ are the setting and measured output of additional experiment $i$.
• $\mathcal D_r=\mathcal D_0\cup\{(a_i,o_i)\}_{i=1}^r$ is the complete measured history after $r$ additional experiments. Repeats are separate records. Here we use $\mathcal D_0$; Part 3b adds experiments.
• $L_h(\theta_h;\mathcal D_r)=p(\mathcal D_r\mid\theta_h,h)$ is the likelihood: the density of the fixed measured outputs, conditional on their known settings, viewed as a function of $\theta_h$.
• $p(\theta_h\mid h)$ is the coefficient prior.

Bayes' rule gives

$Z_h$ normalizes this density over coefficients. It is also called the model evidence, but comparing evidences belongs to Part 3b. The integral uses the prior, not the posterior already fitted to these observations.

An approximate posterior representation need not supply an estimate of $Z_h$. MCMC, for example, can target the normalized posterior without evaluating this integral; NPE learns a normalized conditional density from simulations.

1.1 Approaches to inference#

We can map out a few different approaches to inference, (note we've already used exact inference because we could):

• Exact inference / numerical quadrature: evaluate posterior integrals analytically when possible, or approximate them by weighted evaluations on a grid. A grid becomes expensive as the number of parameters grows. Our toy has an analytic reference.
• MCMC: Explore a posterior with a chain of samples, avoiding the need to calculate an integral over parameters. But if you actually want to calculate the evidence you need to do extra work. We do not build a standalone MCMC sampler here, but the Metropolis moves inside Murphy's SMC are MCMC steps.
• Importance sampling / SMC: Approximate distributions with weighted parameter samples, we explore below in Sections 2 and 3 as this was used extensively in Murphy's parameter inference.
• Variational inference: Fit a chosen distribution to approximate the posterior, we don't explore here, for now.
• Neural Posterior Estimation (NPE): Learn posteriors for parameters conditional on simulated observations as in program-synthesis SBI parameter inference. See Section 4..
• Neural likelihood / ratio estimation: learn observation densities or density ratios from simulations. The model-ratio route is in Part 3b.

Distinguish the likelihood cannot be evaluated from a likelihood can be evaluated but its parameter integral is hard. We start in the second case, using an exactly solvable toy as a reference. The neural route deliberately learns from simulations without calling that analytic reference during training.

1.2 Representing the posterior#

We have the ParameterPosterior object from before, where we specifically looked at the GaussianPosterior. We are now going to add ParticlePosterior objects for sampling methods and GaussianNeuralPosterior for learned methods.

These store and expose parameter moments (eg. mean and covariance in the Gaussian case), have a sample method to allow sampling from the distribution as well as predict to allow prediction of values at a given input or design, for instance predicting noise free means. They do not all promise a smooth density or a Gaussian shape. The neural version here happens to be a conditional Gaussian, it doesn't have to be.

for cls in (GaussianPosterior, ParticlePosterior, GaussianNeuralPosterior):
    print(cls.__name__, "inherits", cls.__bases__[0].__name__)
    assert issubclass(cls, ParameterPosterior)
print(inspect.getsource(ParameterPosterior))

SIGMA = 0.2
INITIAL_X = (-2.0, -1.0, 0.0, 1.0)
MENU_X = (-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0)
world = LineWorld(a=2, b=1, sigma=SIGMA, grid=MENU_X)
pool = tuple(PolynomialHypothesis(degree, prior_scale=3, name=name)
             for degree, name in enumerate(("constant", "line", "quadratic")))
measurement_rng = np.random.default_rng(7)
data = tuple(world.run(Design.from_dict({"x": x}), measurement_rng) for x in INITIAL_X)
exact_parameters = {h.name: fit_parameters(h, data, SIGMA) for h in pool}
table([{"x": obs.design["x"], "observed y": obs.y[0]} for obs in data])
assert len(data) == 4

GaussianPosterior inherits ParameterPosterior
ParticlePosterior inherits ParameterPosterior
GaussianNeuralPosterior inherits ParameterPosterior
class ParameterPosterior(ABC):
    """Coefficient uncertainty for the current polynomial examples.

    Predictions are noise-free moments, not necessarily Gaussian densities.
    Samples have shape (n, hypothesis.n_params). This deliberately does not
    require density evaluation: a weighted particle measure has no smooth PDF.
    """
    hypothesis: PolynomialHypothesis
    mean: Array
    covariance: Array

    @abstractmethod
    def sample(self, rng: np.random.Generator, n: int) -> Array: ...

    @abstractmethod
    def predict(self, inputs: Array) -> tuple[Array, Array]: ...

2. Prior importance sampling#

We have a fixed model $h$ and fixed data $\mathcal D_0$, and want a posterior over its coefficients $\theta$.

Draw $N$ parameter vectors from the prior, $\theta^{(j)}\sim p(\theta\mid h)$, and evaluate the observed-data likelihood at each. Here $j=1,\ldots,N$ indexes parameter particles, not observations. Define their log likelihoods $\lambda_j$ and normalized weights $W_j$:

Because the proposal is the prior, the prior/proposal density ratio cancels. The weighted particles approximate the parameter posterior, with mean $\sum_j W_j\theta^{(j)}$.

Weight effective sample size (ESS) is $1/\sum_jW_j^2$. It ranges from 1 (all weight on one particle) to $N$ (equal weights). It measures weight concentration, not how well the particles cover parameter space. High ESS does not establish that we have found the posterior.

We can inspect the particles and weights returned by polynomial_importance_fit. Its ParticleFit record also contains a normalizing-constant estimate, used in Part 3b.

# Set up world and initial measurements
SIGMA = 0.2
INITIAL_X = (-2.0, -1.0, 0.0, 1.0)
MENU_X = (-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0)
world = LineWorld(a=2, b=1, sigma=SIGMA, grid=MENU_X)
measurement_rng = np.random.default_rng(7)
data = tuple(world.run(Design.from_dict({"x": x}), measurement_rng) for x in INITIAL_X)

# Propose a hypothesis
line = PolynomialHypothesis(1, prior_scale=3, name="line")

# Use the importance_fitting function that PolynomialParticleUpdater would use to update Beliefs.
prior_fit: ParticleFit = polynomial_importance_fit(line, data, SIGMA, n_particles=1500,
                           rng=np.random.default_rng(31))
trace = prior_fit.steps[0]
largest = np.argsort(trace.weights)[-8:][::-1]

# Fit each model's exact parameter posterior for comparison.
exact_parameters = {h.name: fit_parameters(h, data, SIGMA) for h in pool}

table([{"particle": int(j), "intercept": trace.particles_before[j, 0],
        "slope": trace.particles_before[j, 1],
        "log likelihood": trace.log_likelihoods[j], "weight": trace.weights[j]}
       for j in largest])
print("Parameter weight ESS:", trace.ess, "of", len(trace.weights))
print("Estimated / exact coefficient mean:", prior_fit.posterior.mean, exact_parameters["line"].mean)

Parameter weight ESS: 4.109384153487241 of 1500
Estimated / exact coefficient mean: [0.83149 1.90979] [0.92219 1.93291]

Clearly, after an importance sampling step we haven't managed to exactly find the posterior of the parameters. The challenge with IS is that it is reliant on our prior proposal distribution being a good choice, so that our particles are actually exploring useful areas of parameter space.

This motivates a method that permits movement within the sampling space, rather than being strongly influenced by the initial draw from a prior. Hence, the use of Sequential Monte Carlo (SMC).

3. Sequential Monte Carlo: move to the posterior by introducing the observed likelihood gradually#

SMC reduces our reliance on the initial prior draws by repeatedly reweighting, resampling and moving the particles. The prior remains part of the model. Movement helps reach regions missed by the initial draws, but finite runs can still miss posterior modes, and zero prior support cannot be overcome.

We keep model $h$ and observed data $\mathcal D_0$ fixed. Write $L(\theta)=p(\mathcal D_0\mid\theta,h)$ for the likelihood of those observations.

3.1 Move from prior to posterior#

We gradually increase $\beta$ from 0 to 1:

At $\beta=0$ this is the prior; at $\beta=1$ it is the posterior. Intermediate values introduce the likelihood gradually, without collecting new data.

3.2 Reweight the particles#

When we increase $\beta$ by $\Delta\beta$, we multiply each incoming particle weight by the additional likelihood factor, then normalize:

Here $\theta_j$ is a current parameter particle, $W_j$ its incoming normalized weight, and $\widetilde W_j$ its new weight before resampling. We use only the increase in likelihood power: the previous target already included the earlier power; this reweighting is taking in to account the size of the step and the likelihood at the new position of the particle to reweight, so the higher likelihood particles are upweighted.

3.3 Resample, then move#

Resample particles according to their new weights, then reset the weights to $1/N$, where $N$ is the population size. This copies promising particles, but doesn't create new parameter values.

Next, use Metropolis–Hastings moves to explore nearby values while preserving the current target $\pi_\beta$. These moves account for both the prior and the tempered likelihood. They help prevent the population from becoming just a collection of duplicates.

3.4 Put the steps together#

Our implementation resamples at every temperature, so incoming weights are always $1/N$. The reweighting calculation is:

increment = (next_beta - beta) * ell
weights = np.exp(increment - logsumexp(increment))

Here ell contains current log likelihoods, and logsumexp evaluates the log of the sum of exponentials stably.

Repeat increase $\beta$ → reweight → resample → move until $\beta=1$. Below we inspect a fixed schedule. The code also tracks normalizing constants; their use for evidence-based model comparison is explained in Part 3b. Murphy's parameter-SMC algorithm motivates this component; model search and latent-state inference remain separate.

smc_fit: ParticleFit = polynomial_tempered_smc(
    line, data, SIGMA, n_particles=1500, rng=np.random.default_rng(31),
    temperatures=[0, 0.001, 0.005, 0.02, 0.08, 0.25, 0.6, 1],
    moves=5,
)
table([{"rung": k+1, "beta": step.beta, "ESS before resampling": step.ess,
        "MH acceptance": step.acceptance_rate}
       for k, step in enumerate(smc_fit.steps)])
print("Full-dataset likelihood evaluations:", smc_fit.likelihood_evaluations)
assert smc_fit.steps[-1].beta == 1
print("SMC / exact coefficient mean:", smc_fit.posterior.mean, exact_parameters["line"].mean)

Full-dataset likelihood evaluations: 54000
SMC / exact coefficient mean: [0.91614 1.92836] [0.92219 1.93291]

The next plot shows how the parameter population changes. Compare its final location and spread with the exact posterior, rather than judging it only by the particle weights.

4. Learn the posterior instead#

The idea here is to learn the posterior: given noisy outputs, predict a distribution over the coefficients that could have produced them.

For each candidate model $h$, create a separate collection of $S$ simulated datasets:

• $\theta^{(s)}$ is a coefficient vector drawn from model $h$'s prior; for our simple line example, it contains an intercept and a slope.
• $\mathbf x=(-2,-1,0,1)$ is the ordered list of known settings. $\widetilde{\mathbf y}^{(s)}$ contains one simulated noisy output at each setting, all generated using the same $\theta^{(s)}$.
• $s$ indexes simulated datasets, not individual experiments. $\mathcal S_h$ is this synthetic collection, not the real observed dataset $\mathcal D_0$.

Our estimator $q_\psi(\theta\mid\mathbf y,h,\mathbf x)$ is a normalized density over coefficients $\theta$, given outputs $\mathbf y$. Its trainable parameters $\psi$ are the network's weights, biases, and covariance parameters—not the scientific coefficients $\theta$.

Use the first $S_{\rm train}=\lfloor0.8S\rfloor$ pairs for training and hold out the rest for validation. Below, $S=3000$ per model: 2400 training pairs and 600 validation pairs. Fit $\psi$ by minimizing

This is the mean negative log density, labelled NLL in the code. We reward the estimator for assigning high density to the known generating coefficients. This is not the observation likelihood: we score coefficients given outputs, not outputs given coefficients. No exact posterior labels are needed.

After training, we can supply $\mathbf y_{\rm obs}$, the actual outputs from $\mathcal D_0$ in the same setting order. The resulting $q_\psi(\theta\mid\mathbf y_{\rm obs},h,\mathbf x)$ approximates our parameter posterior.

Here PyTorch learns a Gaussian with a linear mean and a full covariance matrix. This family can represent the true toy posterior; we could have instead used a flow, which will be implemented in the code for our final fair comparison. Each estimator is fitted for one model and one ordered setting list. We must retrain when that list changes.

Note something important: the neural posterior estimator is trained to capture the behaviour of a hypothesis, and then effectively "transferred" by inputting the true observations. If the hypothesis is really far from the truth, the observations will be far from the training data for this parameter set, so we might expect high error.

Here we call train_npe and condition_npe directly, so no model-comparison network is trained. Part 3b combines parameter inference with model weights through NeuralUpdater.

print(inspect.getsource(ConditionalGaussian))
npe_rng = np.random.default_rng(40)
started = perf_counter()
npe_fits = {
    h.name: train_npe(h, INITIAL_X, SIGMA, n_simulations=3000, max_iter=220,
                      seed=int(npe_rng.integers(0, 2**31-1)))
    for h in pool
}
observed_y = np.array([obs.y[0] for obs in data])
neural_parameters = {h.name: condition_npe(npe_fits[h.name], h, observed_y) for h in pool}
table([{"model": name, "simulated datasets for NPE": fit.simulations,
        "initial train NLL": fit.train_losses[0],
        "final train NLL": fit.train_losses[-1], "validation NLL": fit.validation_loss}
       for name, fit in npe_fits.items()])
print("NPE fitting time (no model comparison):", round(perf_counter()-started, 3), "seconds")
assert all(isinstance(p, GaussianNeuralPosterior) for p in neural_parameters.values())

class ConditionalGaussian(nn.Module):
    """Inspectable NPE: linear conditional mean and a shared learned Cholesky factor.

    A neural density estimator need not start as a deep network. For this toy,
    the exact conditional mean is linear and covariance does not depend on y.
    A richer family is needed for nonlinear means or multimodal posteriors.
    Outputs are in prior-standardized coefficient coordinates.
    """
    def __init__(self, n_observations, n_parameters):
        super().__init__()
        self.mean_head = nn.Linear(n_observations, n_parameters, dtype=torch.float64)
        nn.init.zeros_(self.mean_head.weight)
        nn.init.zeros_(self.mean_head.bias)
        self.raw_cholesky = nn.Parameter(torch.zeros(n_parameters, n_parameters, dtype=torch.float64))

    def scale_tril(self):
        """Exponentiate the diagonal to make covariance positive definite."""
        raw = self.raw_cholesky
        return torch.tril(raw, diagonal=-1) + torch.diag(torch.exp(torch.diag(raw)))

    def forward(self, standardized_y):
        """Return q(theta_standardized | y), not a distribution over observations."""
        return torch.distributions.MultivariateNormal(
            self.mean_head(standardized_y), scale_tril=self.scale_tril())

NPE fitting time (no model comparison): 0.499 seconds

A low training loss is not enough.

We make two different checks below. Neither uses new real-world observations.

Coverage: does the uncertainty behave as expected? Generate 1000 fresh simulated datasets from the line model, using coefficients drawn from its prior. These datasets were not used for training or validation. For each one:

1. Give its noisy outputs to the trained NPE estimator.
2. Form a 95% posterior interval for each coefficient: posterior mean $\pm\,1.96$ posterior standard deviations. This formula applies because our estimator returns a Gaussian.
3. Check whether the interval contains the coefficient that actually generated that dataset.

Coverage is the fraction of these 1000 intervals that contain their generating coefficient. For example, 950 successes means 95% coverage. The table and left plot report this separately for the intercept and slope; these are marginal intervals, not a joint region for both coefficients. Values near 95% are the target, with some variation from the finite test sample. Lower coverage means the intervals miss too often; higher coverage means they contain the generating values more often than the stated level.

This checks performance averaged over our prior and simulator. It does not establish accuracy for our particular real dataset, or for a world that differs from the simulator. Even returning the prior every time can achieve 95% average coverage while ignoring the observations. We therefore also compare posterior means and covariances with the exact solution later in the post.

Training loss: did fitting improve the estimator? The right plot tracks the average negative log density assigned to the known training coefficients; lower is better. Each curve is a separate model's estimator, not a score for choosing that model. Coefficients are scaled by their prior standard deviation during training, so these losses are not model evidences.

The horizontal axis counts optimizer trial evaluations, including rejected trials—not epochs or new datasets. Temporary spikes can occur. The symmetric-log vertical scale compresses large magnitudes while retaining zero and negative losses. Negative values are possible because a probability density can exceed one.

5. Check the inferred coefficients without choosing a model#

We can compare each approximate posterior with the exact posterior under the same assumed model and data. The constant model also has a well-defined posterior even though it cannot represent the world's slope.

The error in the mean is the L2 norm of the difference between coefficient means. The error in the covariance is the Frobenius norm: the square root of the sum of squared entrywise differences. Zero means agreement with the exact reference.

We can see that different inference methods produce different results, and that is also dependent on the Hypothesis model under test. But we should remember that numerical error is not the same as scientific uncertainty, which remains even under exact inference.

is_rng, smc_rng = np.random.default_rng(31), np.random.default_rng(31)
parameter_methods = {
    "Exact": exact_parameters,
    "Prior IS": {h.name: polynomial_importance_fit(h, data, SIGMA, n_particles=1500, rng=is_rng).posterior for h in pool},
    "Tempered SMC": {h.name: polynomial_tempered_smc(h, data, SIGMA, n_particles=1500, rng=smc_rng).posterior for h in pool},
    "NPE": neural_parameters,
}
table([{"method": method, "assumed model": h.name,
        "mean error L2": np.linalg.norm(fits[h.name].mean-exact_parameters[h.name].mean),
        "covariance error Frobenius": np.linalg.norm(fits[h.name].covariance-exact_parameters[h.name].covariance)}
       for method, fits in parameter_methods.items() for h in pool])
assert len(data) == 4

Next: which model should receive our trust?#

We now have parameter posteriors conditional on each model, but they don't tell us how much probability to assign to the models themselves, which we need to perform model selection from the Hypothesis pool, or indeed know whether ot expand the pool.

Continue to Part 3b: comparing and selecting models.