Higher-dimensional fitting

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 there is another sine component present:

\[y = A_1 \sin\left(2 \pi \left(\frac{t}{P_1} + t_1\right)\right) + A_2 \sin\left(2 \pi \left(\frac{t}{P_2} + t_2\right)\right) + B + \epsilon\]
[1]:
import numpy as np
from numpy import sin, pi

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

def sine_model2(t, B, A1, P1, t1, A2, P2, t2):
    return A1 * sin((t / P1 + t1) * 2 * pi) + A2 * sin((t / P2 + t2) * 2 * pi) + B

The model has four unknown parameters per component:

  • the signal offset \(B\)

  • the amplitude \(A\)

  • the period \(P\)

  • the time offset \(t_0\)

As we will see, the second component makes the 7-dimensional parameter space already quite challenging to explore.

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_model2(t, B=1.0, A1=4.2, P1=3, t1=0, A2=1.2, P2=1.2, t2=1.2), 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-highd_6_0.svg

A beautiful noisy data set, but we can see the 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

parameters2 = ['B', 'A1', 'P1', 't1', 'A2', 'P2', 't2']

def prior_transform2(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]

    # let amplitude go from 0.01 to 100
    params[4] = 10**(cube[4] * 3 - 1)
    # let period go from 0.3 to 30
    params[5] = 10**(cube[5] * 2)
    # let time go from 0 to 1
    params[6] = cube[6]
    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_likelihood2(params):
    # unpack the current parameters:
    B, A1, P1, t1, A2, P2, t2 = params

    # avoid unnecessary multiple solutions:
    #    force ordering by period from large to small
    if P1 < P2:
        return -1e300

    # compute for each x point, where it should lie in y
    y_model = sine_model2(t, B=B, A1=A1, P1=P1, t1=t1, A2=A2, P2=P2, t2=t2)
    # 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,
    wrapped_params=[False, False, False, True],
)

sampler2 = ultranest.ReactiveNestedSampler(
    parameters2,
    log_likelihood2,
    prior_transform2,
    wrapped_params=[False, False, False, True, False, False, True],
)

First, we solve the simpler model:

[7]:
result1 = sampler1.run(min_num_live_points=400)
sampler1.print_results()
[ultranest] Sampling 400 live points from prior ...
[ultranest] Explored until L=-4e+01
[ultranest] Likelihood function evaluations: 54002
[ultranest]   logZ = -52.15 +- 0.1601
[ultranest] Effective samples strategy satisfied (ESS = 2108.9, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.45+-0.06 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.37, need <0.5)
[ultranest]   logZ error budget: single: 0.19 bs:0.16 tail:0.01 total:0.16 required:<0.50
[ultranest] done iterating.

logZ = -52.141 +- 0.369
  single instance: logZ = -52.141 +- 0.191
  bootstrapped   : logZ = -52.150 +- 0.369
  tail           : logZ = +- 0.010

    B                   1.04 +- 0.15
    A1                  3.93 +- 0.21
    P1                  3.074 +- 0.052
    t1                  0.18 +- 0.36

Now lets have a go at the harder problem. We limit the number of evaluations:

[8]:
result2 = sampler2.run(min_num_live_points=400, max_ncalls=400000)
[ultranest] Sampling 400 live points from prior ...
Z=-180.7(0.00%) | Like=-170.47..-37.25 [-197.7326..-67.1018] | it/evals=2894/244107 eff=1.1875% 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))
[ultranest] Explored until L=-4e+01
[ultranest] Likelihood function evaluations: 400343
/home/user/.local/lib/python3.6/site-packages/numpy/core/_methods.py:205: RuntimeWarning: overflow encountered in multiply
  x = um.multiply(x, x, out=x)
[ultranest] Reached maximum number of likelihood calls (400343 > 400000)...
[ultranest] done iterating.

The efficiency is very low. This is not just because of the dimensionality of the problem, but also because of the degeneracies. To make progress, lets use a slice sampler:

[9]:
import ultranest.stepsampler

# have to choose the number of steps the slice sampler should take
# after first results, this should be increased and checked for consistency.

nsteps = 2 * len(parameters2)
# create step sampler:
sampler2.stepsampler = ultranest.stepsampler.RegionSliceSampler(nsteps=nsteps)

# alternatively, we can let the sample identify the number of steps needed on the fly:
# This is done by requiring the point to move at least the typical distance
# between live points, on average.
#sampler2.stepsampler = ultranest.stepsampler.RegionSliceSampler(nsteps=400, adaptive_nsteps='move-distance')


# run again:
result2 = sampler2.run(min_num_live_points=400)
sampler2.print_results()
[ultranest] Widening roots to 400 live points (have 400 already) ...
[ultranest] Explored until L=-2e+01
[ultranest] Likelihood function evaluations: 976364
/home/user/.local/lib/python3.6/site-packages/numpy/core/_methods.py:205: RuntimeWarning: overflow encountered in multiply
  x = um.multiply(x, x, out=x)
[ultranest]   logZ = -44.98 +- 0.1615
[ultranest] Effective samples strategy satisfied (ESS = 2924.5, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.05 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.36, need <0.5)
[ultranest]   logZ error budget: single: 0.23 bs:0.16 tail:0.01 total:0.16 required:<0.50
[ultranest] done iterating.

logZ = -45.019 +- 0.270
  single instance: logZ = -45.019 +- 0.234
  bootstrapped   : logZ = -44.982 +- 0.270
  tail           : logZ = +- 0.010

    B                   1.00 +- 0.15
    A1                  4.21 +- 0.23
    P1                  3.07 +- 0.45
    t1                  0.06 +- 0.19
    A2                  1.16 +- 0.25
    P2                  1.243 +- 0.039
    t2                  0.246 +- 0.049

The efficiency is now constant (at 1/nsteps).

Plot the parameter posterior probability distribution

A classic corner plot:

[10]:
from ultranest.plot import cornerplot
cornerplot(result1)
_images/example-sine-highd_21_0.svg
[11]:
cornerplot(result2)
_images/example-sine-highd_22_0.svg
[12]:
sampler1.ncall
[12]:
54002

Plot the fit

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

[13]:
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)


[13]:
<matplotlib.collections.PolyCollection at 0x7f7d4c364dd8>
_images/example-sine-highd_26_1.svg
[14]:
plt.figure()
plt.title("2-sine fit")
plt.xlabel('x')
plt.ylabel('y')
plt.errorbar(x=t, y=y, yerr=yerr,
             marker='o', ls=' ', color='orange')

band = PredictionBand(t_grid)

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

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)
[14]:
<matplotlib.collections.PolyCollection at 0x7f7d4c27dd68>
_images/example-sine-highd_27_1.svg

Model comparison

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 second 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, and 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:

[15]:
K = np.exp(result2['logz'] - result1['logz'])
K
[15]:
1238.8914927925082

This tells us, assuming both models are equally probable a-priori, that the 2-sine model is 150 times more probable to be the true model than the 1-sine model.

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. Here, the factor is extremely large, so we can be fairly confident that the 2-sine model is correct.

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).