Simple example ============== This example is supposed to serve as a quick start to the ModelFitting module of FlavorPy. Import FlavorPy --------------- After installing FlavorPy with :code:`pip install flavorpy`, import the modelfitting module. .. code:: ipython3 # 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 :math:`M_e = \begin{pmatrix} 0.0048\times0.0565 & 0 & 0 \\ 0 & 0.0565 & 0 \\ 0 & 0 & 1\end{pmatrix} \quad` and :math:`\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: .. code:: ipython3 # 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. .. code:: ipython3 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: .. code:: ipython3 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: .. code:: ipython3 MyModel.get_obs({'v1': 1.5, 'v2': 1.1, 'v3': 1.3}) .. parsed-literal:: {'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 :math:`\sin^2(\theta_{ij})`, ‘d/pi’ is the cp violating phase in the PMNS matrix divided by :math:`\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. .. code:: ipython3 pd.set_option('display.max_columns', None) # This pandas setting allows us to see all columns df = MyModel.make_fit(points=5) df .. raw:: html
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 :math:`\chi^2` arround 5. Since :math:`\chi^2=x` can be interpreted as the specific point lying in the :math:`\sqrt{x}\,\sigma` confidence level region, this means that our point is outside the 2\ :math:`\sigma` but inside the 3\ :math:`\sigma` region of the experimental data. Let us take a look at the individual contributions to :math:`\chi^2` for the best-fit point by .. code:: ipython3 MyModel.print_chisq(df.loc[0]) .. parsed-literal:: '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 :math:`\sin^2\theta_{12}`, :math:`\sin^2\theta_{13}`, :math:`\Delta m_{21}^2`, and :math:`\Delta m_{3\ell}^2` are within their experimental 1\ :math:`\sigma` intervall. However, :math:`\sin^2\theta_{23}` and the CP phase :math:`\delta_{\mathrm{CP}}/\pi` are only within their experimental 2\ :math:`\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\ :math:`\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. .. code:: ipython3 df_mcmc = MyModel.mcmc_fit(df.loc[[0]], mcmc_steps=2000) df_mcmc = MyModel.complete_fit(df_mcmc) .. parsed-literal:: 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\ :math:`\sigma` CL region. The colormap is scaled in such a way that the 1\ :math:`\sigma` region, i.e. \ :math:`\chi^2<1`, is green, the 2\ :math:`\sigma` region with :math:`1<\chi^2<4` is yellow, the 3\ :math:`\sigma` region with :math:`4<\chi^2<9` is orange, and anything white lies outside of the :math:`5\sigma` region. .. code:: ipython3 mf.plot(df_mcmc, x='v1', y='v2', vmin=0, vmax=25); .. image:: simpleexample_modelfitting_v1_v2.png We can also plot the observables and their corresponding 1\ :math:`\sigma` and 3\ :math:`\sigma` CL bounds from NuFit v5.3, i.e. http://www.nu-fit.org/?q=node/278, .. code:: ipython3 mf.plot(df_mcmc, x='s12^2', y='m21^2', show_exp='1dim', vmin=0, vmax=25); .. image:: simpleexample_modelfitting_s12_m21.png but also the two dimensional :math:`\chi^2`-profiles to get a better estimate of the CL regions for correlated, non-gaussian errors, e.g. .. code:: ipython3 mf.plot(df_mcmc, x='s23^2', y='d/pi', show_exp='2dim', vmin=0, vmax=25, gridsize=8); .. image:: simpleexample_modelfitting_s23_dpi.png Note that our model spreads out the full 1\ :math:`\sigma` region of the well measured :math:`\sin^2\theta_{12}` and :math:`\sin^2\theta_{13}` while having a clear prediction for :math:`\sin^2\theta_{23}` and the CP-angle :math:`\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: .. code:: ipython3 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'); .. image:: simpleexample_modelfitting_lobster.png