Model comparison

Model

We again consider the sine model with gaussian measurement errors.

\[y = A_1 \sin\left(2 \pi \left(\frac{t}{P_1} + t_1\right)\right) + B + \epsilon\]

where \(\epsilon \sim \mathrm{Normal}(0, \sigma)\)

We want to test if this is preferred over pure noise.

[1]:
import numpy as np
from numpy import pi, sin

def sine_model1(t, B, A1, P1, t1):
    return A1 * sin((t / P1 + t1) * 2 * pi) + B

def sine_model0(t, B):
    return B + t*0

The model has four unknown parameters per component:

  • the signal offset \(B\)

  • the amplitude \(A\)

  • the period \(P\)

  • the time offset \(t_0\)

Generating data

Lets generate some data following this model:

[2]:
np.random.seed(42)

n_data = 50

# time of observations
t = np.random.uniform(0, 5, size=n_data)
# measurement values
yerr = 1.0
y = np.random.normal(sine_model1(t, B=1.0, A1=0.9, P1=3, t1=0), yerr)

Visualise the data

Lets plot the data first to see what is going on:

[3]:
%matplotlib inline
import matplotlib.pyplot as plt

plt.figure()
plt.xlabel('x')
plt.ylabel('y')
plt.errorbar(x=t, y=y, yerr=yerr,
             marker='o', ls=' ', color='orange')
t_range = np.linspace(0, 5, 1000)

_images/example-sine-modelcomparison_6_0.svg

A beautiful noisy data set, with some hints of a modulation.

Now the question is: what model parameters are allowed under these data?

First, we need to define the parameter ranges through a prior:

[4]:
parameters1 = ['B', 'A1', 'P1', 't1']

def prior_transform1(cube):
    # the argument, cube, consists of values from 0 to 1
    # we have to convert them to physical scales

    params = cube.copy()
    # let background level go from -10 to +10
    params[0] = cube[0] * 20 - 10
    # let amplitude go from 0.1 to 100
    params[1] = 10**(cube[1] * 3 - 1)
    # let period go from 1 to 100
    params[2] = 10**(cube[2] * 2)
    # let time go from 0 to 1
    params[3] = cube[3]
    return params

parameters0 = ['B']

def prior_transform0(cube):
    # the argument, cube, consists of values from 0 to 1
    # we have to convert them to physical scales

    params = cube.copy()
    # let background level go from -10 to +10
    params[0] = cube[0] * 20 - 10
    return params

Define the likelihood, which measures how far the data are from the model predictions. More precisely, how often the parameters would arise under the given parameters. We assume gaussian measurement errors of known size (yerr).

\[\chi^2 = \sum\left(\frac{m_i-y_i}{\sigma}\right)^2\]
\[\log \cal{L} = -\chi^2 / 2\]

where the model is the sine_model function from above at time \(t_i\).

[5]:
def log_likelihood1(params):
    # unpack the current parameters:
    B, A1, P1, t1 = params

    # compute for each x point, where it should lie in y
    y_model = sine_model1(t, B=B, A1=A1, P1=P1, t1=t1)
    # compute likelihood
    loglike = -0.5 * (((y_model - y) / yerr)**2).sum()

    return loglike

def log_likelihood0(params):
    B, = params

    y_model = sine_model0(t, B=B)
    # compute likelihood
    loglike = -0.5 * (((y_model - y) / yerr)**2).sum()

    return loglike

Solve the problem:

[6]:
import ultranest

sampler1 = ultranest.ReactiveNestedSampler(parameters1, log_likelihood1, prior_transform1)

sampler0 = ultranest.ReactiveNestedSampler(parameters0, log_likelihood0, prior_transform0)

[7]:
result1 = sampler1.run(min_num_live_points=400)
sampler1.print_results()
[ultranest] Sampling 400 live points from prior ...
[ultranest] Explored until L=-2e+01
[ultranest] Likelihood function evaluations: 207256
[ultranest]   logZ = -32.7 +- 0.1341
[ultranest] Effective samples strategy satisfied (ESS = 2447.0, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.04 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.13, need <0.5)
[ultranest]   logZ error budget: single: 0.16 bs:0.13 tail:0.01 total:0.13 required:<0.50
[ultranest] done iterating.

logZ = -32.731 +- 0.227
  single instance: logZ = -32.731 +- 0.156
  bootstrapped   : logZ = -32.697 +- 0.227
  tail           : logZ = +- 0.010
insert order U test : converged: True correlation: inf iterations

    B                   : -2.65 │ ▁▁   ▁ ▁ ▁▁▁▇▇▁▁▁▁▁▁▁▁ ▁   ▁        ▁ │7.51      1.02 +- 0.20
    A1                  : 0.10  │▁▁▂▅▇▆▃▁▁▁▁▁▁▁ ▁▁▁ ▁  ▁▁             ▁ │6.79      0.89 +- 0.27
    P1                  : 1.0   │▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▁▁▁▁▁▁▁▁▁│100.0     3.8 +- 6.3
    t1                  : 0.00  │▇▆▅▄▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▄▆▇│1.00      0.47 +- 0.44

[8]:
result0 = sampler0.run(min_num_live_points=400)
sampler0.print_results()
[ultranest] Sampling 400 live points from prior ...
[ultranest] Explored until L=-3e+01
[ultranest] Likelihood function evaluations: 3406
[ultranest]   logZ = -35.6 +- 0.06298
[ultranest] Effective samples strategy satisfied (ESS = 1258.2, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.07 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.07, need <0.5)
[ultranest]   logZ error budget: single: 0.09 bs:0.06 tail:0.03 total:0.07 required:<0.50
[ultranest] done iterating.

logZ = -35.617 +- 0.167
  single instance: logZ = -35.617 +- 0.092
  bootstrapped   : logZ = -35.602 +- 0.164
  tail           : logZ = +- 0.033
insert order U test : converged: True correlation: inf iterations

    B                   : 0.57  │ ▁  ▁▁▁▁▁▁▂▂▂▂▄▆▆▇▇▇▇▇▆▆▅▅▄▂▂▂▁▁▁▁▁▁▁▁ │1.69      1.15 +- 0.14

Plot the parameter posterior probability distribution

A classic corner plot of the parameter pairs and the marginal distributions:

[9]:
from ultranest.plot import cornerplot
cornerplot(result1)
_images/example-sine-modelcomparison_16_0.svg
[10]:
cornerplot(result0)
_images/example-sine-modelcomparison_17_0.svg

If you want, you can also play with the posterior as a pandas frame:

[11]:
import pandas as pd
df = pd.DataFrame(data=result1['samples'], columns=result1['paramnames'])
df.describe()
[11]:
B A1 P1 t1
count 7073.000000 7073.000000 7073.000000 7073.000000
mean 1.015893 0.885090 3.776838 0.467878
std 0.204398 0.270386 6.260011 0.442664
min -2.188261 0.100294 1.129196 0.000017
25% 0.911866 0.731955 2.889296 0.045377
50% 1.007681 0.882176 3.031760 0.140482
75% 1.114623 1.042846 3.177238 0.953596
max 7.049534 6.467633 99.410089 0.999911

Plot the fit:

To evaluate whether the results make any sense, we want to look whether the fitted function goes through the data points.

[12]:
plt.figure()
plt.title("1-sine fit")
plt.xlabel('x')
plt.ylabel('y')
plt.errorbar(x=t, y=y, yerr=yerr,
             marker='o', ls=' ', color='orange')


t_grid = np.linspace(0, 5, 400)

from ultranest.plot import PredictionBand
band = PredictionBand(t_grid)

# go through the solutions
for B, A1, P1, t1 in sampler1.results['samples']:
    # compute for each time the y value
    band.add(sine_model1(t_grid, B=B, A1=A1, P1=P1, t1=t1))

band.line(color='k')
# add 1 sigma quantile
band.shade(color='k', alpha=0.3)
# add wider quantile (0.01 .. 0.99)
band.shade(q=0.49, color='gray', alpha=0.2)


[12]:
<matplotlib.collections.PolyCollection at 0x7f0fb745d570>
_images/example-sine-modelcomparison_22_1.svg

Model comparison methods

We now want to know:

Is the model with 2 components better than the model with one component?

What do we mean by “better” (“it fits better”, “the component is significant”)?

  1. Which model is better at predicting data it has not seen yet?

  2. Which model is more probably the true one, given this data, these models, and their parameter spaces?

  3. Which model is simplest, but complex enough to capture the information complexity of the data?

Bayesian model comparison

Here we will focus on b, and apply Bayesian model comparison.

For simplicity, we will assume equal a-prior model probabilities.

The Bayes factor is:

[13]:
K = np.exp(result1['logz'] - result0['logz'])
print("K = %.2f" % K)
print("The 1-sine model is %.2f times more probable than the no-signal model" % K)
print("assuming the models are equally probable a priori.")
K = 17.92
The 1-sine model is 17.92 times more probable than the no-signal model
assuming the models are equally probable a priori.

N.B.: Bayes factors are influenced by parameter and model priors. It is a good idea to vary them and see how sensitive the result is.

For making decisions, thresholds are needed. They can be calibrated to desired low false decisions rates with simulations (generate data under the simpler model, look at K distribution).

Calibrating Bayes factor thresholds

Lets generate some data sets under the null hypothesis (noise-only model) and see how often we would get a large Bayes factor. For this, we need to fit with both models.

[14]:
import logging
logging.getLogger('ultranest').setLevel(logging.FATAL)
[15]:
K_simulated = []

import logging
logging.getLogger('ultranest').handlers[-1].setLevel(logging.FATAL)

# go through 100 plausible parameters
for B in sampler0.results['samples'][:10]:
    # generate new data
    y = np.random.normal(sine_model0(t, B=1.0), yerr)

    # analyse with sine model
    sampler1 = ultranest.ReactiveNestedSampler(parameters1, log_likelihood1, prior_transform1)
    Z1 = sampler1.run(viz_callback=False)['logz']
    # analyse with noise-only model
    sampler0 = ultranest.ReactiveNestedSampler(parameters0, log_likelihood0, prior_transform0)
    Z0 = sampler0.run(viz_callback=False)['logz']
    # store Bayes factor
    K_here = Z1 - Z0
    K_simulated.append(K_here)
    print()
    print("Bayes factor: %.2f" % np.exp(K_here))

Z=-31.7(82.75%) | Like=-25.70..-25.22 [-25.7024..-25.7019]*| it/evals=3089/144234 eff=2.1476% N=400
/home/user/Downloads/UltraNest/ultranest/integrator.py:1732: UserWarning: Sampling from region seems inefficient (0/40 accepted in iteration 2500). To improve efficiency, modify the transformation so that the current live points are ellipsoidal, or use a stepsampler, or set frac_remain to a lower number (e.g., 0.5) to terminate earlier.
  warnings.warn(warning_message)
Z=-30.5(95.84%) | Like=-26.48..-26.48 [-26.4768..-26.4768]*| it/evals=2880/3376 eff=96.7742% N=400
Bayes factor: 0.38
Z=-21.2(96.52%) | Like=-17.10..-17.10 [-17.0996..-17.0996]*| it/evals=2970/3481 eff=96.3973% N=400
Bayes factor: 0.26
Z=-24.9(95.65%) | Like=-20.91..-20.91 [-20.9139..-20.9139]*| it/evals=2840/3341 eff=96.5658% N=400
Bayes factor: 0.23
Z=-31.7(87.95%) | Like=-26.13..-24.62 [-26.1311..-26.1309]*| it/evals=3140/172663 eff=1.8228% N=400
/home/user/Downloads/UltraNest/ultranest/integrator.py:1732: UserWarning: Sampling from region seems inefficient (0/40 accepted in iteration 2500). To improve efficiency, modify the transformation so that the current live points are ellipsoidal, or use a stepsampler, or set frac_remain to a lower number (e.g., 0.5) to terminate earlier.
  warnings.warn(warning_message)
Z=-30.4(95.61%) | Like=-26.26..-26.25 [-26.2561..-26.2561]*| it/evals=2880/3376 eff=96.7742% N=400
Bayes factor: 0.30
Z=-30.4(92.05%) | Like=-24.85..-23.40 [-24.8471..-24.8470]*| it/evals=3265/228929 eff=1.4287% N=400
/home/user/Downloads/UltraNest/ultranest/integrator.py:1732: UserWarning: Sampling from region seems inefficient (0/40 accepted in iteration 2500). To improve efficiency, modify the transformation so that the current live points are ellipsoidal, or use a stepsampler, or set frac_remain to a lower number (e.g., 0.5) to terminate earlier.
  warnings.warn(warning_message)
Z=-29.1(96.33%) | Like=-24.97..-24.97 [-24.9683..-24.9683]*| it/evals=2960/3457 eff=96.8270% N=400
Bayes factor: 0.29
Z=-23.8(95.93%) | Like=-19.81..-19.81 [-19.8090..-19.8090]*| it/evals=2880/3378 eff=96.7092% N=400
Bayes factor: 0.27
Z=-32.4(96.32%) | Like=-28.37..-28.37 [-28.3677..-28.3677]*| it/evals=2920/3416 eff=96.8170% N=400
Bayes factor: 0.26
Z=-37.2(96.83%) | Like=-33.21..-33.21 [-33.2074..-33.2074]*| it/evals=2960/3473 eff=96.3228% N=400
Bayes factor: 0.33
Z=-26.7(88.64%) | Like=-20.16..-18.14 [-20.1598..-20.1592]*| it/evals=3725/274699 eff=1.3580% N=400
/home/user/Downloads/UltraNest/ultranest/integrator.py:1732: UserWarning: Sampling from region seems inefficient (0/40 accepted in iteration 2500). To improve efficiency, modify the transformation so that the current live points are ellipsoidal, or use a stepsampler, or set frac_remain to a lower number (e.g., 0.5) to terminate earlier.
  warnings.warn(warning_message)
Z=-25.8(96.30%) | Like=-21.79..-21.79 [-21.7871..-21.7871]*| it/evals=2920/3443 eff=95.9579% N=400
Bayes factor: 0.45
Z=-26.6(96.24%) | Like=-22.38..-22.38 [-22.3770..-22.3770]*| it/evals=2970/3490 eff=96.1165% N=400
Bayes factor: 1.15
[16]:
plt.figure()
plt.hist(np.exp(K_simulated), histtype='step', label='From simulated noise data')
ylo, yhi = plt.ylim()
plt.vlines(K, ylo, yhi, label='From our real data')
plt.xscale('log')
plt.xlabel('Bayes factor')
plt.ylabel('Frequency')
plt.legend(loc='upper center');
_images/example-sine-modelcomparison_30_0.svg

If we run this a bit longer, we will fill in the simulation histogram better. But already now we can see:

We are using simulations to measure how often, by chance, we would see a Bayes factor higher than the one we observe. By building up a histogram, we can get a p-value, telling us our false decision rate for any Bayes factor threshold. Thus, we are putting a frequentist property on our Bayesian inference-based decision.

So I would say: Pure noise does not produce as high a Bayes factor as we see it in the real data.

Calibrating Bayes factor thresholds reduces the dependence on model priors and model parameter priors.

Further reading

[ ]: