
We use the \(\chi ^2\)  function, given by the sum of the individual \(\chi _i^2\)  of each measured observable \(i\) ,
\(\chi ^2(\vec{x})=\sum _i \chi _i^2,\) 
where \(\vec{x}\)  is a vector consisting of the input parameters of model 1 or 2, listed in Eq. (REF ) and (REF ), respectively. For the observables \(\theta _{23}^{\text{PMNS}}\)  and \(\delta ^{\text{PMNS}}\)  we use the exact \(\Delta \chi ^2\)  function provided by NuFIT 5.0 [1]}. For all other observables we compute the individual \(\chi _i^2\)  function as
\(\chi _i^2=\frac{(x_i^{\text{pred}}-x_i^{\text{exp}})^2}{\sigma _{i\pm }^2},\) 
where \(x_i^{\text{pred}}\)  is our prediction for the observable \(i\)  and \(x_i^{\text{exp}}\)  is its central experimental value, whereas \(\sigma _i\)  denotes the corresponding standard deviation which may be asymmetric, i.e. \(\sigma _{i+}\)  for \(x_i^{\text{exp}}>x_i^{\text{pred}}\)  and \(\sigma _{i-}\)  for \(x_i^{\text{exp}}>x_i^{\text{pred}}\) , respectively.

We use a differential evolution algorithm to calculate the best-fit point of a certain model.

Using an adaptive Metropolis-Hastings algorithm [2]} to perform an MCMC analysis we determine the posterior density of the observables of a given model. With flat prior probability distributions starting from the best-fit point we compute 8 independent chains with \(1.25\cdot 10^5\)  data points giving us a total of \(10^6\)  points for each dataset.

