A Toy Example and Some Components

To help kick off the discussion, let's describe some general components of discovery, whether scientific or otherwise. These form the basis for the code we will run to compare, and is a lens thorugh which we can look at the work in some of the references.

This is very much a toy example. It exists in part as an introduction to some ideas, in part to introduce ideas in the code, and also in some way as a mini primer on exact inference. If you are comfortable with everything, feel free to skip ahead.

Some components of scientific discovery#

1. The World. The system our scientist interacts with. The scientist may have ideas about how it works (more below), but in principle all it can ever obtain are observations $o$ in response to carefully chosen actions or interventions $a$.
2. Hypotheses. Proposed models of the world. E.g. given a pendulum of length $l$ and mass $m$ in gravity $g$, what is its period $T$? One hypothesis — the correct one for a point mass at small angles — is $T = 2\pi\sqrt{l/g}$. In practice there are many more candidates, depending on the exact physical conditions. Ideas for what to test don't appear from a vacuum: in the approaches I look at here, an LLM's training data plus post-training reasoning ideally give it a good prior over models. The component that creates and updates the list of hypotheses is the Proposer, it has access to the Beliefs about the hypotheses.
3. Experiments. The actions or interventions the system can take to uncover the world's behaviour. Given everything observed so far — historical data $\mathcal{D}$ plus experiment/observation pairs $\{a_i, o_i\}$ — how do we choose the next one? This is the meaty bit. It can be handed directly to a reasoning system with no formalisation on top, or done with a strictly Bayesian approach and no LLM at all. Either way, the component is the Designer.
4. Observations. Fairly obvious; what is returned by the world when an experiment is performed.
5. Updating beliefs. The scientist gets an observation back and revises its belief about the true explanation. This can be done by formal machinery external to the scientist, left entirely to the scientist internally, or — as I'll also explore — something in between, where the agent is prompted to consider several different structures for handling this step. The component is the Updater.
6. We can also consider an Evaluator.

The final outcome is the scientist's belief about the correct hypothesis at the end of its turn budget, including its uncertainty.

Note that these sub-components are in part just conceptual. While we will make them concrete here to discuss the overlap between papers exploring these ideas, in principle, large parts could be handled entirely internally by a reasoning model. But it's useful to break it down.

Putting these together, one turn of the loop runs: the Proposer proposes hypotheses given the current Beliefs; the Designer chooses the next experiment $a$; the World returns an observation $o$; the Updater revises the Beliefs over the pool; and the Evaluator scores the result against held-out behaviour. This repeats until the turn budget is spent.

Part 1: A Toy Example#

In this post we're first going to play with a toy example to concretely ground what we mean by different components of the AI Scientist. Here, there is no agent involved in proposing hypotheses or designing experiments, but this is designed to allow us to "swap in" a reasoning model later. Sections 1.2-1.6 act as a small primer in exact inference and what we mean by likelihood, evidence, and posteriors, so skip straight over if you're already comfortable.

1.1 Setting up some machinery#

The simplest set up we can usefully imagine is a "World" consisting of a system of unknown polynomial equation. i.e. the hidden world has form $y=\sum_{k=0}^{d}\theta_k X^k+\varepsilon$ where $X$ represent the values at which observations are made (i.e. the set of $X$ is the set of experiments performed) and $y$ are the returned observations. $\theta$ represents the set of parameters of a given hypothesis. Here, a hypothesis itself is which degree $d$ which is unknown, and in some sense the "meta-hypothesis" is that we should be trying to fit a polynomial at all. Turtles all the way down.

To make this concrete, the hidden drop down cells here contain code pulling out some of these components. First, the Designs (experiments), and resulting Observations whihc together make a Dataset.

Formally, we are talking about $\mathcal D$, the dataset or history $[(a_1,o_1),(a_2,o_2),(a_3,o_3)]$. Each action $a_i$ sets input $x_i$, and observation $o_i$ contains measured value $y_i$.

"""First, the Design, Observation, and sequence of Observations that form the Dataset."""

from dataclasses import dataclass, field
from typing import Any, Optional, Protocol, Sequence, runtime_checkable

import numpy as np
from scipy.special import logsumexp

Array = np.ndarray
RNG = np.random.Generator

@dataclass(frozen=True)
class Design:
    """One executable experimental configuration (an intervention).

    `params` is a tuple of (name, value) pairs rather than a dict so Designs
    are hashable — useful for dedup and for logging which design was chosen
    each round.
    """

    params: tuple[tuple[str, float], ...]
    name: str = ""

    @staticmethod
    def from_dict(d: dict[str, float], name: str = "") -> "Design":
        return Design(params=tuple(sorted(d.items())), name=name)

    def as_dict(self) -> dict[str, float]:
        return dict(self.params)

    def __getitem__(self, key: str) -> float:
        return self.as_dict()[key]

@dataclass(frozen=True)
class Observation:
    """The outcome of running one design in the real world (noisy)."""

    design: Design
    y: Array  # observed outcome; shape is world-specific but fixed per world

Dataset = Sequence[Observation]

Next, we shall include the fact that our experiments may return observations with some noise. This is particularly relevant for any experimental observations. We set this up as noise that is added in to the observations of the World.

@runtime_checkable
class NoiseModel(Protocol):
    """Observation model linking noise-free hypothesis simulation output to data.

    Kept separate from Hypothesis so the same structure can be scored under
    different observation models (deterministic Gaussian now; Poisson spike
    counts or learned-summary-statistic likelihoods later, like NeuronBench setup).
    """

    def log_likelihood(self, y_obs: Array, y_pred: Array) -> float:
        """Log p(y_obs | y_pred)."""
        ...

    def sample(self, y_pred: Array, rng: RNG) -> Array:
        """Draw a noisy observation around the noise-free prediction."""
        ...


@dataclass(frozen=True)
class GaussianNoise(NoiseModel):
    """Isotropic Gaussian observation noise.

    Can compute the density of the observation given the prediction, 
    or be sample from the conditional posterior of observations given prediction.
    """

    sigma: float

    def log_likelihood(self, y_obs: Array, y_pred: Array) -> float:
        r = np.asarray(y_obs, float) - np.asarray(y_pred, float)
        n = r.size
        return float(
            -0.5 * np.sum(r**2) / self.sigma**2
            - n * np.log(self.sigma)
            - 0.5 * n * np.log(2.0 * np.pi)
        )

    def sample(self, y_pred: Array, rng: RNG) -> Array:
        return np.asarray(y_pred, float) + self.sigma * rng.standard_normal(
            np.shape(y_pred)
        )

Now, we can define an interface for a World, and give some very simple examples, i.e. our toy model, called LineWorld. We will later expand this to actually interesting worlds; for now we have defined an instance (our scientist is unable to read this initialisation, it can only see observations) with some semi-random initialization, drawn sample observations to make an initial Dataset and made a plot.

The components are the design_space which exposes which interventions can be run, and the run method itself.

import numpy as np

@runtime_checkable
class World(Protocol):
    """A ground-truth environment the agent experiments on.

    Note for an agent: The agent may call `run` at most `budget` times (AgentConfig.budget);
    `test_designs` are held-out interventions used only to measure interventional prediction.
    """

    name: str
    noise: NoiseModel

    def description(self) -> str:
        """Natural-language context handed to the Proposer, the description of the domain. We can contaminate by
        rewriting exactly this string (obfuscated names/units)."""
        ...

    def design_space(self) -> Sequence[Design]:
        """The discrete menu of interventions the agent may choose from."""
        ...

    def run(self, design: Design, rng: RNG) -> Observation:
        """Execute one experiment: ground-truth mechanism + noise."""
        ...

    def test_designs(self) -> Sequence[Design]:
        """Held-out interventions for evaluation. i.e. agent can't see these,
        they are to produce clean observations at test time.
        """
        ...

@dataclass
class LineWorld:
    """Hidden mechanism y = a*x + b, Gaussian noise, x chosen from a grid."""

    a: float = 2.0
    b: float = -1.0
    sigma: float = 0.1
    grid: tuple[float, ...] = (-2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 3.0)
    name: str = "line-world"

    def __post_init__(self) -> None:
        self.noise: NoiseModel = GaussianNoise(self.sigma)

    def description(self) -> str:
        return (
            "A scalar quantity y is measured after setting a single control "
            "variable x. The mechanism relating x to y is unknown. "
            "Measurements carry independent Gaussian noise."
        )

    def design_space(self) -> Sequence[Design]:
        return [Design.from_dict({"x": v}, name=f"x={v:g}") for v in self.grid]

    def _truth(self, design: Design) -> Array:
        # You ideally don't want to call this -- it's cheating!
        return np.array([self.a * design["x"] + self.b])

    def run(self, design: Design, rng: RNG) -> Observation:
        return Observation(design=design, y=self.noise.sample(self._truth(design), rng))

    def test_designs(self) -> Sequence[Design]:
        # Held-out interventions OUTSIDE the training grid on purpose:
        # interventional evaluation should include extrapolation.
        return [Design.from_dict({"x": v}, name=f"x={v:g}") for v in (-4.0, 5.0)]

SIGMA = 0.2
first_world = LineWorld(a=2.0, b=1.0, sigma=SIGMA, grid=(-2, -1, 0, 1))
first_rng = np.random.default_rng(7)
history: Dataset = tuple(first_world.run(d, first_rng)
                for d in first_world.design_space())

xs, ys = [], []
for observation in history:
    xs.append(observation.design["x"])
    ys.append(float(np.asarray(observation.y).reshape(-1)[0]))

Another component is the set of Hypothesis objects. Here we specify a super general protocol, but also a specific polynomial hypothesis. In practice, a scientist won't necessarily restrict to a very specific structure, but the methods of the Hypothesis object make clear what is needed. Here there needs to be a prior over parameters, a method to measure the log density, log_prior, and a method to sample from that prior. Most importantly, there must be a way to simulate the consequences of the hypothesis given the experiment design, hence the simulate method.

import sympy

@runtime_checkable
class Hypothesis(Protocol):
    """One candidate structure with free parameters theta.

    theta is always a flat float vector of length `n_params`; structure
    lives in code/symbolic form, parameters in theta, for now.

    This is meant to be a completely general way to represent any model of the world.

    In practice, this might be nice to extend to the concept of "meta-hypothesis",
    which is a language explanation with consequences for underlying model structures,
    without specifying precise the model. But that's a TODO for me :)
    """

    name: str

    @property
    def n_params(self) -> int: ...

    @property
    def param_names(self) -> Sequence[str]: ...

    def sample_prior(self, rng: RNG, n: int) -> Array:
        """(n, n_params) draws from the parameter prior."""
        ...

    def log_prior(self, theta: Array) -> float:
        """Log prior density of one theta vector."""
        ...

    def simulate(self, design: Design, theta: Array) -> Array:
        """Noise-free predicted outcome of `design` under parameters theta.

        Must return the same shape as Observations for the target world.
        This might be expensive and we should track the cost of these calls.

        The expense of course depends on the world being tested.
        """
        ...

    def sympy_form(self) -> Optional[Any]:
        """Symbolic form (sympy.Expr) if the model is a closed-form law,
        else None. Used by eval for symbolic-equivalence scoring
        (eg. ChemBench's sympy check); simulators (ODE worlds) return None.
        """
        ...

    def human_description(self) -> Optional[str]:
        """A human-readable description, or None when one is unavailable.

        Implementations must provide this method even when sympy_form exists.
        It can explain the model's variables and assumptions beyond its formula.
        """
        ...

@dataclass
class PolynomialHypothesis(Hypothesis):
    """Implements Hypothesis for y(x) = sum_k theta_k * x^k.

    The coefficients have independent N(0, prior_scale^2) priors.
    Explicit inheritance makes the implemented interface visible here.
    """

    degree: int
    prior_scale: float = 3.0
    name: str = ""

    def __post_init__(self) -> None:
        if not self.name:
            self.name = f"poly-deg{self.degree}"

    @property
    def n_params(self) -> int:
        return self.degree + 1

    @property
    def param_names(self) -> Sequence[str]:
        return [f"c{k}" for k in range(self.n_params)]

    def sample_prior(self, rng: RNG, n: int) -> Array:
        return self.prior_scale * rng.standard_normal((n, self.n_params))

    def log_prior(self, theta: Array) -> float:
        theta = np.asarray(theta, float)
        return float(
            -0.5 * np.sum((theta / self.prior_scale) ** 2)
            - self.n_params * np.log(self.prior_scale)
            - 0.5 * self.n_params * np.log(2.0 * np.pi)
        )

    def simulate(self, design: Design, theta: Array) -> Array:
        x = design["x"]
        powers = np.array([x**k for k in range(self.n_params)])
        return np.array([float(np.dot(np.asarray(theta, float), powers))])

    def sympy_form(self) -> Optional[Any]:
        x = sympy.Symbol("x")
        cs = sympy.symbols(f"c0:{self.n_params}")
        return sum(c * x**k for k, c in enumerate(cs))

    def human_description(self) -> str:
        return (
            f"A degree-{self.degree} polynomial in x: y = {self.sympy_form()}. "
            "The coefficients have independent Gaussian priors centered on zero "
            f"with standard deviation {self.prior_scale:g}. "
            "Predictions exclude observation noise."
        )

Here, we are going to choose a set of hypotheses explicitly, thus acting as the Proposer. We also have our dataset of current Observations.

constant = PolynomialHypothesis(0, prior_scale=3.0, name="constant")
linear = PolynomialHypothesis(1, prior_scale=3.0, name="linear")
quadratic = PolynomialHypothesis(2, prior_scale=3.0, name="quadratic")
models = [constant, linear, quadratic]
print("Known inputs x:", xs)
print("Measured outputs y:", ys)

Known inputs x: [-2, -1, 0, 1]
Measured outputs y: [-2.9997539693285034, -0.9402508924983061, 0.9451724289275565, 2.8218816322485454]

1.2 Identifying the best hypothesis in the exact inference case#

Now, we have three measurements and two candidate explanations (we'll ignore the other for now). There are two unknowns we need to discover:

1. Which stucture? Is constant or linear good enough or do we need quadratic? This is the choice of hypothesis, $h$, from the set $H$.
2. Which coefficients $\theta$, given $h$? For the simplest case, what intercept and slope of line, and their uncertainty.

For this toy example, a candidate $h$ predicts

Here $d_h$ is the degree chosen by the hypothesis, $\theta_k$ is its coefficient of $x^k$, and $\varepsilon_i$ is independent measurement noise. The hypothesis also specifies a coefficient prior. We have chosen independent zero-mean Gaussian priors with standard deviation $\tau=3$. Measurement-noise standard deviation is known here: $\sigma=0.2$.

Importantly, we are in effect doing inference on two levels here: 1) to find the distribution that best represents each hypothesis' parameters (if we assume that is the correct model) 2) to find the distribution over the actual set of hypotheses themselves.

I will use lowercase $x_i$ for a known input setting, and uppercase $X_h$ for the design matrix built from all the observed inputs for hypothesis $h$. Its row $i$ is $[1,x_i,\ldots,x_i^{d_h}]$.

The complete available dataset is still $\mathcal D=[(a_i,o_i)]_{i=1}^3$. No new measurements are made while we fit or compare these candidates. Everything below uses the objects defined above; no project modules or saved datasets are imported.

Let's inspect the design matrices for each of the hypothesis models

def observed_arrays(data: Dataset) -> tuple[Array, Array]:
    """Read known scalar inputs and measured outputs, retaining record order."""
    inputs, outputs = [], []
    for obs in data:
        if set(obs.design.as_dict()) != {"x"} or np.size(obs.y) != 1:
            raise ValueError("This example requires a known x and one measured y per record.")
        inputs.append(obs.design["x"])
        outputs.append(float(np.asarray(obs.y).reshape(-1)[0]))
    x_values, y_values = np.asarray(inputs, float), np.asarray(outputs, float)
    if not np.all(np.isfinite(x_values)) or not np.all(np.isfinite(y_values)):
        raise ValueError("Inputs and measurements must be finite.")
    return x_values, y_values

def design_matrix(hypothesis: PolynomialHypothesis, inputs: Array) -> Array:
    """Return one row per input and one column per polynomial coefficient."""
    return np.vander(np.asarray(inputs, float),
                     N=hypothesis.n_params, increasing=True)


x_observed, y_observed = observed_arrays(history)
np.set_printoptions(precision=5, suppress=True)
for h in models:
    X_h = design_matrix(h, x_observed)
    print(f"{h.name}: {h.sympy_form()}, design matrix shape={X_h.shape}")
    print(X_h)

constant: c0, design matrix shape=(4, 1)
[[1.]
 [1.]
 [1.]
 [1.]]
linear: c0 + c1*x, design matrix shape=(4, 2)
[[ 1. -2.]
 [ 1. -1.]
 [ 1.  0.]
 [ 1.  1.]]
quadratic: c0 + c1*x + c2*x**2, design matrix shape=(4, 3)
[[ 1. -2.  4.]
 [ 1. -1.  1.]
 [ 1.  0.  0.]
 [ 1.  1.  1.]]

1.3 Computing the likelihood of the observations for a single hypothesis#

The noise model allows is to define the distribution of observations given the hypothesis, with parameters $\theta$:

This is the likelihood, where we have of course already assumed the mean of observations follows the relationship $f$. For $n$ independent measurements, log likelihoods add:

$\mathbf y$ is the vector of the $n$ measured outputs. The final sum is the squared residual error. All terms, not just the residuals, belong to the density.

A small terminology clarification for the earlier noise class: its sample method draws from this observation distribution conditional on a mean. It does not by itself draw from a posterior over coefficients. We'll get to that posterior next.

These are probability densities, not probabilities of exact points.

Let's compute the log likelihood explicitly.

def dataset_log_likelihood(
    hypothesis: Hypothesis, theta: Array, data: Dataset, noise: NoiseModel
) -> float:
    """Score fixed observations under the means predicted by one coefficient vector."""
    return float(sum(
        noise.log_likelihood(obs.y, hypothesis.simulate(obs.design, theta))
        for obs in data
    ))


observation_noise = GaussianNoise(SIGMA)
for theta in [np.array([1.0, 0.0]), np.array([1.0, 2.0])]:
    predictions = design_matrix(linear, x_observed) @ theta # this is a linear prediction now, as the design matrix makes that possible
    score = dataset_log_likelihood(linear, theta, history, observation_noise)
    direct = (-len(history)/2*np.log(2*np.pi) - len(history)*np.log(SIGMA)
              - np.sum((y_observed-predictions)**2)/(2*SIGMA**2))
    np.testing.assert_allclose(score, direct)
    print("Coefficients [intercept, slope]:", theta)
    print("Predictions:", predictions, "; log likelihood:", round(score, 5))

Coefficients [intercept, slope]: [1. 0.]
Predictions: [1. 1. 1. 1.] ; log likelihood: -285.7988
Coefficients [intercept, slope]: [1. 2.]
Predictions: [-3. -1.  1.  3.] ; log likelihood: 2.28322

Great, so for our linear hypothesis $h_{linear}$, we can compute the likelihood in the case of different possible parameter sets, given the data. In practice, choosing those parameters is also part of hypothesis creation but here we are treating the hypothesis selection over the model degree, $d$.

1.4 Obtaining the posterior over coefficients for a single hypothesis#

So the next thing we need is actually to get the coefficient posterior, where all we have are these likelihoods and the prior.

Here we are claiming that our coefficients are Gaussian-distributed, we now want the distribution over all plausible vectors after seeing out Dataset of Observations, $\mathcal{D}$.

For a fixed hypothesis with $p=d_h+1$ coefficients, our assumptions are

$I_p$ and $I_n$ are identity matrices of sizes $p$ and $n$.

Bayes' rule says posterior is proportional to likelihood times prior. Keeping the terms that depend on $\theta$, its log density is

Expand the square and collect terms. Define the matrix $A$ and vector $q$:

The terms involving $\theta$ are $\theta^\top A\theta-2\theta^\top q$. Choose $b$ so that $Ab=q$. Completing the square gives

The last term does not depend on $\theta$. We can therefore recognise the posterior as

$b$ is the entire length-$p$ coefficient-mean vector, not just the intercept or the world's b setting. $V$ is the $p\times p$ coefficient covariance.

Here, we create an object to store the posterior for a single fixed hypothesis $h$, and compute the mean and covariance in the linear hypothesis case.

The code solves $Ab=q$ and $AV=I_p$, rather than explicitly calculating a matrix inverse.

This calculation is exact for the polynomial/Gaussian assumptions. It does not make those assumptions universally correct.

@dataclass
class ParameterPosterior:
    """Gaussian coefficient uncertainty conditional on one fixed hypothesis."""
    hypothesis: PolynomialHypothesis
    mean: Array
    covariance: Array

    def predict(self, inputs: Array) -> tuple[Array, Array]:
        """Return noise-free output predicted mean and variance at each requested input."""
        X = design_matrix(self.hypothesis, inputs)
        mean = X @ self.mean
        variance = np.sum((X @ self.covariance) * X, axis=1)
        return mean, np.maximum(variance, 0.0)


def fit_parameters(
    hypothesis: PolynomialHypothesis, data: Dataset, sigma: float
) -> ParameterPosterior:
    """Refit from this model's prior using the full current observed dataset."""
    if not np.isfinite(sigma) or sigma <= 0:
        raise ValueError("Noise SD must be positive and finite.")
    tau = hypothesis.prior_scale
    if not np.isfinite(tau) or tau <= 0:
        raise ValueError("Coefficient-prior SD must be positive and finite.")
    inputs, outputs = observed_arrays(data)
    X = design_matrix(hypothesis, inputs)
    identity = np.eye(hypothesis.n_params)
    A = identity/tau**2 + X.T @ X/sigma**2
    q = X.T @ outputs/sigma**2
    covariance = np.linalg.solve(A, identity)
    mean = np.linalg.solve(A, q)
    return ParameterPosterior(hypothesis, mean, covariance)


line_posterior = fit_parameters(linear, history, SIGMA)
print("Line coefficient order:", linear.param_names)
print("Posterior mean:", line_posterior.mean)
print("Posterior covariance:\n", line_posterior.covariance)
print("Posterior SDs:", np.sqrt(np.diag(line_posterior.covariance)))

Line coefficient order: ['c0', 'c1']
Posterior mean: [0.92219 1.93291]
Posterior covariance:
 [[0.01198 0.00399]
 [0.00399 0.00799]]
Posterior SDs: [0.10946 0.08939]

The posterior mean is close to intercept 1 and slope 2. We also have uncertainty around those values. For this symmetric set of inputs the line's coefficient covariance is diagonal; it will not generally be diagonal for other inputs or models.

There is a useful distinction here: ordinary least squares finds the coefficient vector that minimises training error. This Bayesian calculation gives a distribution and incorporates the prior. With a finite prior width, its mean need not be the least-squares estimate.

Neither result yet tells us whether the line is more plausible than the quadratic. Both assume a particular structure. To compare structures, we need another calculation.

This type of calculation is what would be used by an Updater to choose how to weight the set of running hypotheses, and choose to add another if required, so it is what we must consider next.

1.5 Compare hypothesis structures by integrating out their coefficients, the model evidence#

A model's evidence asks how much density its prior assigned to the measured outputs, after accounting for observation noise. Or put another way it is the likelihood of the data integrated over the prior distribution of model parameters:

We integrate over $\theta$, keeping the observed $\mathbf y$ and known input settings fixed. Importantly, we average over the coefficient prior, not the fitted posterior. A model is not scored only at its best-fitting coefficients, or put another way, we really shouldn't be fitting to the data twice!

For this example we can do the integral exactly. Before observing $\mathbf y$, it is the sum $X_h\theta+\varepsilon$ of two independent Gaussian quantities. Its mean is zero and its covariance is

This $n\times n$ matrix describes possible measured output vectors given the prior. It is not the $p\times p$ posterior coefficient covariance $V$ we just calculated.

Evaluating the Gaussian density at our actual measured vector and taking its logarithm gives

Now the normalising terms must be retained: they change with the model and its prior. The code below uses a log determinant and a linear solve to evaluate this expression.

def log_model_evidence(
    hypothesis: PolynomialHypothesis, data: Dataset, sigma: float
) -> float:
    """Integrate likelihood over the coefficient prior, not the fitted posterior."""
    tau = hypothesis.prior_scale
    if not np.all(np.isfinite([sigma, tau])) or sigma <= 0 or tau <= 0:
        raise ValueError("Noise and coefficient-prior SDs must be positive and finite.")
    inputs, outputs = observed_arrays(data)
    n = len(outputs)
    if n == 0:
        return 0.0  # With no observations, evidence is 1 and log evidence is 0.
    X = design_matrix(hypothesis, inputs)
    C = sigma**2*np.eye(n) + tau**2*(X @ X.T)
    sign, log_determinant = np.linalg.slogdet(C)
    if sign <= 0:
        raise ValueError("The prior-predictive covariance must be positive definite.")
    quadratic_term = outputs @ np.linalg.solve(C, outputs)
    return float(-0.5*(n*np.log(2*np.pi) + log_determinant + quadratic_term))


for h in models:
    print(h.name, "log evidence:", round(log_model_evidence(h, history, SIGMA), 6))

constant log evidence: -234.783415
linear log evidence: -4.529732
quadratic log evidence: -7.823344

1.6 Scoring among models: model weights in the model pool#

A coefficient prior $p(\theta\mid h)$ and a model prior $p(h)$ are different things. We want to know how much weight each hypothesis in the pool deserves when we've integrated over the coefficient prior

For prior model weights $\pi_h=p(h)$, Bayesian model probabilities are $\pi_h * p(\mathcal{D}|h)/ \sum_j \pi_j * p(\mathcal{D}|h_j)$

Here we give each candidate hypothesis $h$ equal prior weight. We calculate in log space using logsumexp, which is a stable way to take the logarithm of a sum of exponentials.

Our Updater will return a Beliefs record. Hypotheses remain model definitions; their fitted parameters and within-pool probabilities live in this separate record. Numerical fields can be absent, so a later text-only updater does not have to invent a Bayesian posterior just to satisfy the interface.

In this particular case we are demonstrating a fixed ExactGaussianUpdater that performs exact inference. In practice, this inference may be much more complex, both the parameter posterior inference, and the computation of model evidence. We leave that for later.

In addition, we have only used one definition of model weights. There are of course better ways to do this as discussed later in this series.

@dataclass
class Beliefs:
    """An assessment of a pool, with optional numerical support."""
    method: str
    assessment: str
    parameters: Optional[dict[str, ParameterPosterior]] = None
    model_probabilities: Optional[dict[str, float]] = None
    log_evidences: Optional[dict[str, float]] = None
    support_kind: Optional[str] = None


class Updater(Protocol):
    """Assess a supplied candidate pool using the complete observed history."""
    def update(self, data: Dataset, hypotheses: Sequence[PolynomialHypothesis]) -> Beliefs:
        ...


@dataclass
class ExactGaussianUpdater(Updater):
    """Exact polynomial inference with known Gaussian noise and equal model priors."""
    sigma: float

    def update(self, data: Dataset, hypotheses: Sequence[PolynomialHypothesis]) -> Beliefs:
        """Fit each model and normalize its evidence within this candidate pool."""
        names = [h.name for h in hypotheses]
        if not names or len(set(names)) != len(names):
            raise ValueError("Supply a non-empty pool with unique candidate names.")
        # This performs the exact inference to get parameter posteriors
        fits = {h.name: fit_parameters(h, data, self.sigma) for h in hypotheses}
        # This computes the log evidence of each hypothesis in the set
        scores = {h.name: log_model_evidence(h, data, self.sigma) for h in hypotheses}
        if not np.all(np.isfinite(list(scores.values()))):
            raise ValueError("Log evidences must be finite.")
        normalizer = logsumexp(list(scores.values()))
        probabilities = {name: float(np.exp(score-normalizer))
                         for name, score in scores.items()}
        return Beliefs(
            method="exact_gaussian",
            assessment="Relative support in this pool, not proof that any candidate is adequate.",
            parameters=fits,
            model_probabilities=probabilities,
            log_evidences=scores,
            support_kind="bayesian_model_probability",
        )


updater = ExactGaussianUpdater(SIGMA)
beliefs = updater.update(history, models)
print(f"{'model':10} {'LS squared error':>18} {'log evidence':>14} {'probability':>13}")
for h in models:
    X = design_matrix(h, x_observed)
    theta_ols = np.linalg.lstsq(X, y_observed, rcond=None)[0]
    best_fit_error = np.sum((y_observed-X @ theta_ols)**2)
    print(f"{h.name:10} {best_fit_error:18.6f} "
          f"{beliefs.log_evidences[h.name]:14.6f} "
          f"{beliefs.model_probabilities[h.name]:13.6f}")
np.testing.assert_allclose(sum(beliefs.model_probabilities.values()), 1.0)

model        LS squared error   log evidence   probability
constant            18.731484    -234.783415      0.000000
linear               0.009721      -4.529732      0.964209
quadratic            0.001367      -7.823344      0.035791

The quadratic can interpolate these three observations, but the line receives about 96% of the model probability. The quadratic receives about 3.6%; the constant receives almost none.

That is not a contradiction. Training error asks whether a model contains a close-fitting coefficient vector. Evidence asks whether the model's prior, averaged over coefficients, predicts data like these. More flexibility can improve the best fit while, but the prior is less predictive. When integrating over the coefficients, the quadratic case gives lower evidence of the actual observations. In other words, a classic example of the bias-variance tradeoff.

These numbers depend on the candidate pool, coefficient priors, model priors, and noise assumption. They are not probabilities that we have found the final scientific explanation. In particular, a pool containing only the constant would give it probability one, despite its visibly poor predictions.

SO, the actual residual errors ultimately also count

constant_only = updater.update(history, [constant])
constant_mean, constant_variance = constant_only.parameters["constant"].predict(x_observed)
print("Constant's probability if it is the only hypothesis:", constant_only.model_probabilities["constant"])
print("Constant's coefficient SD:",
      np.sqrt(constant_only.parameters["constant"].covariance[0, 0]))
print("Residuals at positions of measurements:", (y_observed-constant_mean)/SIGMA)
assert constant_only.model_probabilities["constant"] == 1.0
assert np.max(np.abs((y_observed-constant_mean)/SIGMA)) > 5

Constant's probability if it is the only hypothesis: 1.0
Constant's coefficient SD: 0.09994449069791544
Residuals at positions of measurements: [-14.78282  -4.48531   4.94181  14.32536]

1.7 Make an episode: proposing hypotheses and updating posteriors and evidences in a one round of the loop#

So far, we have acted as the proposer by selecting the three candidates ourselves. To make the order explicit, we'll start with only the constant hypothesis, assess it, add our scripted alternatives, and reassess the expanded pool.

The flow is:

start with an observed history → assess constant hypothesis → propose line and quadratic as extra hypotheses → assess expanded pool

The runner receives observations, hypotheses, a proposer, and an updater. It does not receive the world or held-out answers. Each update refits from the original priors on the full available history; it does not take the old posterior and count the same data again.

There is no adaptive experiment design yet. The three input settings were chosen in the setup. This is one proposal round, with zero new measurements, not one new experimental round.

Later, we'll add lots of things, like experiment choice, and hypothesis proposal that doesn't just come from us, as well as dealing with cases where inference is hard.

def propose_alternatives(data, current_hypotheses, current_beliefs) -> list[Hypothesis]:
    """Hand-chosen additions for this example; no LLM or automated discovery,
    but we can imagine the beliefs and hypotheses being used here."""
    return [
        PolynomialHypothesis(1, prior_scale=3.0, name="linear"),
        PolynomialHypothesis(2, prior_scale=3.0, name="quadratic"),
    ]


def run_episode(data, initial_hypotheses, proposer, updater: Updater) -> tuple[Beliefs]:
    """Assess a pool, add proposals, and reassess on exactly the same data."""
    pool = list(initial_hypotheses)
    before: Beliefs = updater.update(data, pool)
    additions: list[Hypothesis] = proposer(data, tuple(pool), before)
    pool.extend(additions)
    after: Beliefs = updater.update(data, pool)
    return before, after


history_snapshot = [(obs.design, obs.y.copy()) for obs in history]
before, after = run_episode(history, [constant], propose_alternatives, updater)
print("Before:", before.model_probabilities)
print("After: ", after.model_probabilities)
print("Observed records:", len(history), "; new measurements in the episode: 0")

# Candidate expansion changes pool probabilities, not the constant's own evidence.
np.testing.assert_allclose(before.log_evidences["constant"], after.log_evidences["constant"])
np.testing.assert_allclose(before.parameters["constant"].mean, after.parameters["constant"].mean)
np.testing.assert_allclose(
    list(after.model_probabilities.values()), list(beliefs.model_probabilities.values()))
for obs, (design, measured_y) in zip(history, history_snapshot):
    assert obs.design == design
    np.testing.assert_array_equal(obs.y, measured_y)

Before: {'constant': 1.0}
After:  {'constant': 9.688740668414482e-101, 'linear': 0.9642090122802123, 'quadratic': 0.03579098771978728}
Observed records: 4 ; new measurements in the episode: 0

1.8 Evaluation of the current pool of models#

We should judge explanations by their predictions, not just their pool probabilities or model evidence. First distinguish a noise-free world value $f_*$ at a requested input $x_*$ from a future noisy measurement $y_*$.

For candidate $h$, let $v_{*,h}=[1,x_*,\ldots,x_*^{d_h}]^\top$. Its coefficient posterior has mean $b_h$ and covariance $V_h$, so

Averaging across models with weights $w_h$ gives

This is the law of total variance: uncertainty remains both about the coefficients within each explanation and about which explanation to use. For a future measurement with independent noise, $\operatorname{Var}(y_*\mid\mathcal D)=\operatorname{Var}(f_*\mid\mathcal D)+\sigma^2$.

These are predictive quantities, not yet a value-of-information score. A design method must also specify what we want a new observation to teach us.

def predictive_moments(state: Beliefs, inputs: Array) -> dict[str, Array]:
    """Combine noise-free predictions, keeping the two variance sources separate."""
    if (state.support_kind != "bayesian_model_probability"
            or state.parameters is None or state.model_probabilities is None):
        raise ValueError("Numerical predictions require model probabilities and coefficient posteriors.")
    names = list(state.model_probabilities)
    weights = np.array([state.model_probabilities[name] for name in names])
    moments = [state.parameters[name].predict(inputs) for name in names]
    means = np.array([item[0] for item in moments])
    variances = np.array([item[1] for item in moments])
    mean = weights @ means
    within = weights @ variances
    between = weights @ (means-mean)**2
    return {"mean": mean, "within": within, "between": between, "variance": within+between}


future_x_vals = predictive_moments(after, np.array([2.0]))
for label, value in future_x_vals.items():
    print(label, "at x=2:", value[0])
print("Future measured-output variance:", future_x_vals["variance"][0] + SIGMA**2)

mean at x=2: 4.779649065825904
within at x=2: 0.06882927900245274
between at x=2: 0.0018788230212134992
variance at x=2: 0.07070810202366624
Future measured-output variance: 0.11070810202366625

1.9 Keep evaluation separate from fitting#

Because we built this world, we can check forecasts against its true noise-free output at held-out settings. This is an evaluator's privilege, not an input to the proposer or updater.

Only the following evaluation cell reads the hidden mechanism through _truth. We do not append these answers to $\mathcal D$ or use them to update the pool. The calculation is mean squared prediction error against the true mean, not error against another noisy sample.

# Evaluator-only inputs and answers. Never passed to run_episode or update.
evaluation_designs = first_world.test_designs()
evaluation_inputs = np.array([d["x"] for d in evaluation_designs])
evaluation_truth = np.array([first_world._truth(d)[0] for d in evaluation_designs])
old_forecast = predictive_moments(before, evaluation_inputs)["mean"]
new_forecast = predictive_moments(after, evaluation_inputs)["mean"]
error_before = np.mean((old_forecast-evaluation_truth)**2)
error_after = np.mean((new_forecast-evaluation_truth)**2)
print(f"{'held-out x':>12} {'true mean':>12} {'before':>12} {'after':>12}")
for x, truth, old, new in zip(evaluation_inputs, evaluation_truth, old_forecast, new_forecast):
    print(f"{x:12.1f} {truth:12.5f} {old:12.5f} {new:12.5f}")
print(f"Mean squared error: {error_before:.6f} -> {error_after:.6f}")
assert error_after < error_before
assert len(history) == 4

  held-out x    true mean       before        after
        -4.0     -7.00000      1.00024     -6.97142
         5.0     11.00000      1.00024     10.68310
Mean squared error: 81.999519 -> 0.050621

The expanded pool predicts this particular world's held-out means much better. This is not a general result about scientific discovery: we deliberately supplied a candidate containing the true mechanism where we could perform easy exact inference.

A small check#

Suppose we widen only the quadratic hypothesis coefficient prior distribution, so the SD goes from 3 to 30, keeping our measured history the same. If we called the updater, how does this change our beliefs over the posteriors, and the relative model probalities?

wide_quadratic = PolynomialHypothesis(2, prior_scale=30.0, name="quadratic")
changed_models = [constant, linear, wide_quadratic]
changed_beliefs = updater.update(history, changed_models)
print("Quadratic coefficients with original prior:", after.parameters["quadratic"].mean)
print("Quadratic coefficients with wider prior:   ", changed_beliefs.parameters["quadratic"].mean)
print("Original quadratic probability:", after.model_probabilities["quadratic"])
print("Wider-prior quadratic probability:", changed_beliefs.model_probabilities["quadratic"])

Quadratic coefficients with original prior: [ 0.96881  1.88626 -0.04667]
Quadratic coefficients with wider prior:    [ 0.96997  1.8893  -0.04571]
Original quadratic probability: 0.03579098771978728
Wider-prior quadratic probability: 4.768646447397203e-05

Up Next#

1. Choosing which experiment -- the job of the Designer
2. Infering the best hypothesis when the inference is not exact, a better Updater.