Quickstart

In this simple example we will generate some simulated data, and fit them with 3ML.

Let’s start by generating our dataset:

[1]:
from threeML import *

# Let's generate some data with y = Powerlaw(x)

gen_function = Powerlaw()


# Generate a dataset using the power law, and a
# constant 30% error

x = np.logspace(0, 2, 50)

xyl_generator = XYLike.from_function("sim_data", function = gen_function,
                                     x = x,
                                     yerr = 0.3 * gen_function(x))

y = xyl_generator.y
y_err = xyl_generator.yerr
Configuration read from /home/ndilalla/.threeML/threeML_config.yml
INFO:keyring.backend:Loading KWallet
INFO:keyring.backend:Loading SecretService
INFO:keyring.backend:Loading Windows
INFO:keyring.backend:Loading chainer
INFO:keyring.backend:Loading macOS
Using Gaussian statistic (equivalent to chi^2) with the provided errors.

We can now fit it easily with 3ML:

[2]:
fit_function = Powerlaw()

xyl = XYLike("data", x, y, y_err)

parameters, like_values = xyl.fit(fit_function)
Using Gaussian statistic (equivalent to chi^2) with the provided errors.
Best fit values:

result unit
parameter
source.spectrum.main.Powerlaw.K (10.0 -0.8 +0.9) x 10^-1 1 / (cm2 keV s)
source.spectrum.main.Powerlaw.index -2.023 +/- 0.033

Correlation matrix:

1.00-0.85
-0.851.00

Values of -log(likelihood) at the minimum:

-log(likelihood)
data 24.386648
total 24.386648

Values of statistical measures:

statistical measures
AIC 53.028615
BIC 56.597342

Plot data and model:

[3]:
fig = xyl.plot(x_scale='log', y_scale='log')
../_images/notebooks_Quickstart_6_0.png

Compute the goodness of fit using Monte Carlo simulations (NOTE: if you repeat this exercise from the beginning many time, you should find that the quantity “gof” is a random number distributed uniformly between 0 and 1. That is the expected result if the model is a good representation of the data)

[4]:
gof, all_results, all_like_values = xyl.goodness_of_fit()

print("The null-hypothesis probability from simulations is %.2f" % gof['data'])
The null-hypothesis probability from simulations is 0.45

The procedure outlined above works for any distribution for the data (Gaussian or Poisson). In this case we are using Gaussian data, thus the log(likelihood) is just half of a \(\chi^2\). We can then also use the \(\chi^2\) test, which give a close result without performing simulations:

[5]:
import scipy.stats

# Compute the number of degrees of freedom
n_dof = len(xyl.x) - len(fit_function.free_parameters)

# Get the observed value for chi2
# (the factor of 2 comes from the fact that the Gaussian log-likelihood is half of a chi2)
obs_chi2 = 2 * like_values['-log(likelihood)']['data']

theoretical_gof = scipy.stats.chi2(n_dof).sf(obs_chi2)

print("The null-hypothesis probability from theory is %.2f" % theoretical_gof)
The null-hypothesis probability from theory is 0.44

There are however many settings where a theoretical answer, such as the one provided by the \(\chi^2\) test, does not exist. A simple example is a fit where data follow the Poisson statistic. In that case, the MC computation can provide the answer.