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 0.3 to 30
    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]:
import scipy.stats

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: 141014
[ultranest]   logZ = -32.89 +- 0.09381
[ultranest] Effective samples strategy satisfied (ESS = 2581.2, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.06 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.25, need <0.5)
[ultranest]   logZ error budget: single: 0.16 bs:0.09 tail:0.01 total:0.09 required:<0.50
[ultranest] done iterating.

logZ = -32.898 +- 0.147
  single instance: logZ = -32.898 +- 0.155
  bootstrapped   : logZ = -32.885 +- 0.147
  tail           : logZ = +- 0.010

    B                   1.02 +- 0.28
    A1                  0.87 +- 0.32
    P1                  4.0 +- 7.1
    t1                  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: 3563
[ultranest]   logZ = -35.88 +- 0.08832
[ultranest] Effective samples strategy satisfied (ESS = 1302.7, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.09 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.22, need <0.5)
[ultranest]   logZ error budget: single: 0.10 bs:0.09 tail:0.03 total:0.09 required:<0.50
[ultranest] done iterating.

logZ = -35.878 +- 0.194
  single instance: logZ = -35.878 +- 0.095
  bootstrapped   : logZ = -35.880 +- 0.191
  tail           : logZ = +- 0.031

    B                   1.15 +- 0.15

Plot the parameter posterior probability distribution

A classic corner plot:

[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 7131.000000 7131.000000 7131.000000 7131.000000
mean 1.019333 0.869917 3.992960 0.469300
std 0.276567 0.322632 7.055792 0.438929
min -4.892674 0.102989 1.085266 0.000007
25% 0.903343 0.689195 2.882924 0.047505
50% 1.013807 0.862878 3.040515 0.157557
75% 1.117446 1.030139 3.205770 0.949884
max 7.056695 6.645766 98.880853 0.999913

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 0x7f9ccb9d6ac8>
_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 = 19.68
The 1-sine model is 19.68 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=-27.9(95.67%) | Like=-23.92..-23.92 [-23.9174..-23.9173]*| it/evals=2840/3334 eff=96.7962% N=400
Bayes factor: 0.31
Z=-35.0(93.32%) | Like=-28.70..-26.34 [-28.6969..-28.6960]*| it/evals=3768/290813 eff=1.2975% N=400
/home/user/.local/lib/python3.6/site-packages/ultranest/integrator.py:1243: UserWarning: Sampling from region seems inefficient. You can try increasing nlive, frac_remain, dlogz, dKL, decrease min_ess). [0/40 accepted, it=2500]
  " [%d/%d accepted, it=%d]" % (accepted.sum(), ndraw, nit))
Z=-33.5(95.91%) | Like=-29.56..-29.56 [-29.5582..-29.5582]*| it/evals=2840/3342 eff=96.5330% N=400
Bayes factor: 0.23
Z=-35.7(91.13%) | Like=-29.52..-28.42 [-29.5236..-29.5208]*| it/evals=3519/257142 eff=1.3706% N=400
/home/user/.local/lib/python3.6/site-packages/ultranest/integrator.py:1243: UserWarning: Sampling from region seems inefficient. You can try increasing nlive, frac_remain, dlogz, dKL, decrease min_ess). [0/40 accepted, it=2500]
  " [%d/%d accepted, it=%d]" % (accepted.sum(), ndraw, nit))
Z=-34.0(96.63%) | Like=-30.16..-30.16 [-30.1597..-30.1597]*| it/evals=2880/3398 eff=96.0640% N=400
Bayes factor: 0.20
Z=-27.6(95.84%) | Like=-23.71..-23.71 [-23.7145..-23.7145]*| it/evals=2800/3288 eff=96.9529% N=400
Bayes factor: 0.23
Z=-22.6(96.01%) | Like=-18.68..-18.68 [-18.6827..-18.6827]*| it/evals=2840/3339 eff=96.6315% N=400
Bayes factor: 0.27
Z=-29.6(95.79%) | Like=-25.53..-25.52 [-25.5260..-25.5260]*| it/evals=2880/3389 eff=96.3533% N=400
Bayes factor: 0.35
Z=-31.4(81.90%) | Like=-25.83..-24.37 [-25.8312..-25.8296]*| it/evals=2960/246811 eff=1.2012% N=400
/home/user/.local/lib/python3.6/site-packages/ultranest/integrator.py:1243: UserWarning: Sampling from region seems inefficient. You can try increasing nlive, frac_remain, dlogz, dKL, decrease min_ess). [0/40 accepted, it=2500]
  " [%d/%d accepted, it=%d]" % (accepted.sum(), ndraw, nit))
Z=-29.8(96.99%) | Like=-25.81..-25.81 [-25.8148..-25.8148]*| it/evals=3000/3505 eff=96.6184% N=400
Bayes factor: 0.24
Z=-26.2(85.93%) | Like=-19.56..-18.33 [-19.5597..-19.5587]*| it/evals=3440/159903 eff=2.1567% N=400
/home/user/.local/lib/python3.6/site-packages/ultranest/integrator.py:1243: UserWarning: Sampling from region seems inefficient. You can try increasing nlive, frac_remain, dlogz, dKL, decrease min_ess). [0/40 accepted, it=2500]
  " [%d/%d accepted, it=%d]" % (accepted.sum(), ndraw, nit))
Z=-25.8(96.37%) | Like=-21.71..-21.71 [-21.7134..-21.7134]*| it/evals=2960/3457 eff=96.8270% N=400
Bayes factor: 0.80
Z=-27.7(94.59%) | Like=-21.83..-20.66 [-21.8326..-21.8325]*| it/evals=3615/237162 eff=1.5268% N=400
/home/user/.local/lib/python3.6/site-packages/ultranest/integrator.py:1243: UserWarning: Sampling from region seems inefficient. You can try increasing nlive, frac_remain, dlogz, dKL, decrease min_ess). [0/40 accepted, it=2500]
  " [%d/%d accepted, it=%d]" % (accepted.sum(), ndraw, nit))
Z=-26.4(96.29%) | Like=-22.29..-22.29 [-22.2894..-22.2894]*| it/evals=2960/3464 eff=96.6057% N=400
Bayes factor: 0.27
Z=-25.5(90.09%) | Like=-19.93..-18.02 [-19.9288..-19.9284]*| it/evals=3267/193277 eff=1.6938% N=400
/home/user/.local/lib/python3.6/site-packages/ultranest/integrator.py:1243: UserWarning: Sampling from region seems inefficient. You can try increasing nlive, frac_remain, dlogz, dKL, decrease min_ess). [0/40 accepted, it=2500]
  " [%d/%d accepted, it=%d]" % (accepted.sum(), ndraw, nit))
Z=-24.2(96.62%) | Like=-20.14..-20.14 [-20.1402..-20.1401]*| it/evals=2960/3469 eff=96.4484% N=400
Bayes factor: 0.30
[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

[ ]: