Simple example

This example is supposed to serve as a quick start to the ModelFitting module of FlavorPy.

Import FlavorPy

After installing FlavorPy with pip install flavorpy, import the modelfitting module.

# Import the modelfitting module of FlavorPy
import flavorpy.modelfitting as mf

# We will also need numpy and pandas
import numpy as np
import pandas as pd

Defining mass matrices

To define a model of leptons, we start by defining its mass matrices

\(M_e = \begin{pmatrix} 0.0048\times0.0565 & 0 & 0 \\ 0 & 0.0565 & 0 \\ 0 & 0 & 1\end{pmatrix} \quad\) and \(\quad M_n = \begin{pmatrix} 1 & v_2 & v_3 \\ v_1 & 1 & v_3 \\ v_1 & v_2 & 1\end{pmatrix}\)

For this example we have:

# Charged lepton mass matrix
def Me(params):
    v1, v2, v3 = params['v1'], params['v2'], params['v3']
    return np.array([[0.0048*0.0565, 0 ,0], [0, 0.0565, 0], [0, 0, 1]])

# Neutrino mass matrix
def Mn(params):
    v1, v2, v3 = params['v1'], params['v2'], params['v3']
    return np.array([[1, v2, v3], [v1, 1, v3], [v1, v2, 1]])

Defining parameterspace

Next, we define the parameterspace of our model. We therefore construct an empty parameter space and add the parameters to it. When fitting, we will draw random points within this parameter space.

ParamSpace = mf.ParameterSpace()
ParamSpace.add_dim(name='v1')
ParamSpace.add_dim(name='v2')
ParamSpace.add_dim(name='v3')

Constructing Model

Then we can construct the lepton model as follows:

MyModel = mf.FlavorModel(mass_matrix_e=Me, mass_matrix_n=Mn, parameterspace=ParamSpace, ordering='NO')

Now we can determine the masses and mixing observables of a given point in parameter space by:

MyModel.get_obs({'v1': 1.5, 'v2': 1.1, 'v3': 1.3})

{'me/mu': 0.0048,
 'mu/mt': 0.0565,
 's12^2': 0.9053720503789906,
 's13^2': 0.2893377761696659,
 's23^2': 0.5280465699834944,
 'd/pi': 0.0,
 'r': 0.010306101829667999,
 'm21^2': 4.9968698643488846e-05,
 'm3l^2': 0.0048484576874298696,
 'm1': 0.0033969202435270816,
 'm2': 0.007842688683377208,
 'm3': 0.06971367695488995,
 'eta1': 1.0,
 'eta2': 1.0,
 'J': 0.0,
 'Jmax': 0.05585663171203268,
 'Sum(m_i)': 0.08095328588179423,
 'm_b': 0.03822289118358137,
 'm_bb': 0.025548760738177294,
 'nscale': 0.019258907884260646}

Here, ‘me/mu’ is the mass ratio of electron mass divided by muon mass, ‘sij^2’ refers to the mixing angles \(\sin^2(\theta_{ij})\), ‘d/pi’ is the cp violating phase in the PMNS matrix divided by \(\pi\), ‘m21^2’ and ‘m3l^2’ and the squared neutrino mass differences, i.e. mij^2 = m_i^2 - m_j^2, ‘r’ is their quotient r = m21^2 / m3l^2, ‘m1’ and ‘m2’ and ‘m3’ are the neutrino masses, ‘eta1’ and ‘eta2’ are the majorana phases, ‘J’ is the Jarskog determinant, ‘m_b’ and ‘m_bb’ are the effective neutrino masses for beta decay and neutrinoless double beta decay, respectively.

Fitting model to experimental data

Let us now fit this model to a specific experimental data set. As a default the NuFit v5.3 for NO with SK data is used. To fit this model we choose for example 3 randomly drawn points in the parameter space and apply minimization algorithms to these points, in order to find a point that matches the experimental data well. Note that by default 4 minimization algorithms are applied consecutively to all 3 random points such that we get 12 points in the end.

pd.set_option('display.max_columns', None)  # This pandas setting allows us to see all columns

df = MyModel.make_fit(points=5)
df

	chisq	chisq_dimless	v1	v2	v3	n_scale	me/mu	mu/mt	s12^2	s13^2	s23^2	d/pi	r	m21^2	m3l^2	m1	m2	m3	eta1	eta2	J	Jmax	Sum(m_i)	m_b	m_bb	nscale
0	5.651948	5.651702	5.036863	-0.528403	0.701039	1	0.0048	0.0565	0.305000	0.022200	0.499957	1.0	0.029596	0.000074	0.002505	0.002326	0.008920	0.050108	0.0	1.0	4.107236e-18	0.033538	0.061353	0.009208	0.005369	0.006852
1	5.651948	5.651702	5.036863	-0.528403	0.701039	1	0.0048	0.0565	0.305000	0.022200	0.499957	1.0	0.029596	0.000074	0.002505	0.002326	0.008920	0.050108	0.0	1.0	4.107236e-18	0.033538	0.061353	0.009208	0.005369	0.006852
2	5.651948	5.651702	5.036863	-0.528403	0.701039	1	0.0048	0.0565	0.305000	0.022200	0.499957	1.0	0.029596	0.000074	0.002505	0.002326	0.008920	0.050108	0.0	1.0	4.107236e-18	0.033538	0.061353	0.009208	0.005369	0.006852
3	5.651948	5.651702	5.036863	-0.528403	0.701039	1	0.0048	0.0565	0.305000	0.022200	0.499957	1.0	0.029596	0.000074	0.002505	0.002326	0.008920	0.050108	0.0	1.0	4.107236e-18	0.033538	0.061353	0.009208	0.005369	0.006852
4	5.655930	5.655555	5.049672	-0.534924	0.713045	1	0.0048	0.0565	0.309623	0.022200	0.500035	1.0	0.029592	0.000074	0.002506	0.002402	0.008939	0.050113	0.0	1.0	4.124460e-18	0.033679	0.061455	0.009246	0.005457	0.006835
5	5.655930	5.655555	5.049672	-0.534924	0.713045	1	0.0048	0.0565	0.309623	0.022200	0.500035	1.0	0.029592	0.000074	0.002506	0.002402	0.008939	0.050113	0.0	1.0	4.124460e-18	0.033679	0.061455	0.009246	0.005457	0.006835
6	5.655930	5.655555	5.049672	-0.534924	0.713045	1	0.0048	0.0565	0.309623	0.022200	0.500035	1.0	0.029592	0.000074	0.002506	0.002402	0.008939	0.050113	0.0	1.0	4.124460e-18	0.033679	0.061455	0.009246	0.005457	0.006835
7	5.657581	5.657194	5.041598	-0.533464	0.714117	1	0.0048	0.0565	0.309623	0.022280	0.500068	1.0	0.029592	0.000074	0.002506	0.002413	0.008942	0.050114	0.0	1.0	4.131556e-18	0.033737	0.061469	0.009260	0.005469	0.006845
8	5.657581	5.657194	5.041598	-0.533464	0.714117	1	0.0048	0.0565	0.309623	0.022280	0.500068	1.0	0.029592	0.000074	0.002506	0.002413	0.008942	0.050114	0.0	1.0	4.131556e-18	0.033737	0.061469	0.009260	0.005469	0.006845
9	5.657581	5.657194	5.041598	-0.533464	0.714117	1	0.0048	0.0565	0.309623	0.022280	0.500068	1.0	0.029592	0.000074	0.002506	0.002413	0.008942	0.050114	0.0	1.0	4.131556e-18	0.033737	0.061469	0.009260	0.005469	0.006845
10	5.657581	5.657194	5.041598	-0.533464	0.714117	1	0.0048	0.0565	0.309623	0.022280	0.500068	1.0	0.029592	0.000074	0.002506	0.002413	0.008942	0.050114	0.0	1.0	4.131556e-18	0.033737	0.061469	0.009260	0.005469	0.006845
11	5.658516	5.658316	5.050711	-0.535313	0.712793	1	0.0048	0.0565	0.309623	0.022189	0.500027	1.0	0.029604	0.000074	0.002505	0.002400	0.008940	0.050108	0.0	1.0	4.123487e-18	0.033671	0.061447	0.009243	0.005455	0.006833
12	14.019848	14.013724	5.028969	0.702363	-0.526542	1	0.0048	0.0565	0.305000	0.022280	0.500002	0.0	0.029570	0.000074	0.002506	0.002339	0.008921	0.050119	0.0	0.0	0.000000e+00	0.033596	0.061379	0.009223	0.005383	0.006864
13	14.019853	14.013724	5.028969	0.702363	-0.526542	1	0.0048	0.0565	0.305000	0.022280	0.500002	0.0	0.029570	0.000074	0.002506	0.002339	0.008921	0.050119	0.0	0.0	0.000000e+00	0.033596	0.061379	0.009223	0.005383	0.006864
14	14.019853	14.013724	5.028969	0.702363	-0.526542	1	0.0048	0.0565	0.305000	0.022280	0.500002	0.0	0.029570	0.000074	0.002506	0.002339	0.008921	0.050119	0.0	0.0	0.000000e+00	0.033596	0.061379	0.009223	0.005383	0.006864
15	14.019853	14.013724	5.028969	0.702363	-0.526542	1	0.0048	0.0565	0.305000	0.022280	0.500002	0.0	0.029570	0.000074	0.002506	0.002339	0.008921	0.050119	0.0	0.0	0.000000e+00	0.033596	0.061379	0.009223	0.005383	0.006864
16	14.020580	14.014997	5.028979	0.702446	-0.526597	1	0.0048	0.0565	0.305033	0.022281	0.500001	0.0	0.029571	0.000074	0.002506	0.002340	0.008921	0.050119	0.0	0.0	0.000000e+00	0.033597	0.061380	0.009224	0.005383	0.006864
17	14.021824	14.013749	5.028994	0.702405	-0.526469	1	0.0048	0.0565	0.305000	0.022280	0.500000	0.0	0.029565	0.000074	0.002507	0.002339	0.008921	0.050121	0.0	0.0	0.000000e+00	0.033596	0.061382	0.009223	0.005383	0.006864
18	14.021824	14.013749	5.028994	0.702405	-0.526469	1	0.0048	0.0565	0.305000	0.022280	0.500000	0.0	0.029565	0.000074	0.002507	0.002339	0.008921	0.050121	0.0	0.0	0.000000e+00	0.033596	0.061382	0.009223	0.005383	0.006864
19	14.022732	14.014400	5.029001	0.702364	-0.526447	1	0.0048	0.0565	0.304985	0.022280	0.500000	0.0	0.029565	0.000074	0.002507	0.002339	0.008921	0.050121	0.0	0.0	0.000000e+00	0.033595	0.061381	0.009223	0.005383	0.006864

The fit yields a point with \(\chi^2\) arround 5. Since \(\chi^2=x\) can be interpreted as the specific point lying in the \(\sqrt{x}\,\sigma\) confidence level region, this means that our point is outside the 2\(\sigma\) but inside the 3\(\sigma\) region of the experimental data. Let us take a look at the individual contributions to \(\chi^2\) for the best-fit point by

MyModel.print_chisq(df.loc[0])

'me/mu': 0.0048,   chisq: 0.0
'mu/mt': 0.0565,   chisq: 0.0
's12^2': 0.30500000121100895,   chisq: 4.2032172992145e-08
's13^2': 0.02219999995440627,   chisq: 1.0523577958707064e-08
's23^2': 0.4999568484745416,   chisq: 2.9383389552791437
'd/pi': 1.0,   chisq: 2.71318
'm21^2': 7.414848837981547e-05,   chisq: 0.0003285239267222246
'm3l^2': 0.0025053616533515753,   chisq: 0.00010020584059472139
Total chi-square: 5.651947737602212

It looks like the \(\sin^2\theta_{12}\), \(\sin^2\theta_{13}\), \(\Delta m_{21}^2\), and \(\Delta m_{3\ell}^2\) are within their experimental 1\(\sigma\) intervall. However, \(\sin^2\theta_{23}\) and the CP phase \(\delta_{\mathrm{CP}}/\pi\) are only within their experimental 2\(\sigma\) intervall. On a side note, it is no surprise, that the CP phase is off, since it is always CP conserving due to the mass matrices being real. All errors then add up to our best-fit point of the model lying within the 3\(\sigma\) confidence level region.

MCMC fit to explore the minimum of the model

To explore the neighborhood of the minimum we use the emcee marcov chain monte carlo sampler.

df_mcmc = MyModel.mcmc_fit(df.loc[[0]], mcmc_steps=2000)
df_mcmc = MyModel.complete_fit(df_mcmc)

0 : 100%|██████████| 2000/2000 [00:30<00:00, 65.16it/s]

Plotting the results

We can plot the parameterspace and see the contour of the 3\(\sigma\) CL region. The colormap is scaled in such a way that the 1\(\sigma\) region, i.e. \(\chi^2<1\), is green, the 2\(\sigma\) region with \(1<\chi^2<4\) is yellow, the 3\(\sigma\) region with \(4<\chi^2<9\) is orange, and anything white lies outside of the \(5\sigma\) region.

mf.plot(df_mcmc, x='v1', y='v2', vmin=0, vmax=25);

We can also plot the observables and their corresponding 1\(\sigma\) and 3\(\sigma\) CL bounds from NuFit v5.3, i.e. http://www.nu-fit.org/?q=node/278,

mf.plot(df_mcmc, x='s12^2', y='m21^2', show_exp='1dim', vmin=0, vmax=25);

but also the two dimensional \(\chi^2\)-profiles to get a better estimate of the CL regions for correlated, non-gaussian errors, e.g.

mf.plot(df_mcmc, x='s23^2', y='d/pi', show_exp='2dim', vmin=0, vmax=25, gridsize=8);

Note that our model spreads out the full 1\(\sigma\) region of the well measured \(\sin^2\theta_{12}\) and \(\sin^2\theta_{13}\) while having a clear prediction for \(\sin^2\theta_{23}\) and the CP-angle \(\delta_{\mathrm{CP}}^\ell/\pi\), which have not yet been measured that precisely.

Also from the lobster plot we see a prediction for the lightest neutrino mass and the effective neutrion mass for neutrinoless double beta decay:

ax = mf.plot(df_mcmc, x='m1', y='m_bb',
             ordering='NO', show_exp='2dim', xscale='log', yscale='log', vmin=0, vmax=25)
ax.axvspan(0.037,10, facecolor='gray',alpha=0.21); # arXiv: 2009.03287
ax.axhspan(0.156,10, facecolor='gray',alpha=0.21);  # arXiv: 2203.02139
ax.axhspan(0.036,0.156, facecolor='gray',alpha=0.13);  # arXiv: 2203.02139
ax.text(0.1,0.0012,'Cosmology', color='gray');
ax.text(0.0006,0.3,'KamLAND-Zen', color='gray');