425 lines
16 KiB
Markdown
425 lines
16 KiB
Markdown
# Exercise 3 {#sec:3}
|
|
|
|
## Meson decay events generation
|
|
|
|
A number of $N = 50'000$ points on the unit sphere, each representing a
|
|
particle detection event, must be generated according to then angular
|
|
probability distribution function $F_0$:
|
|
|
|
\begin{align*}
|
|
F_0 (\theta, \phi) = \frac{3}{4 \pi} \Bigg[
|
|
\frac{1}{2} (1 - \alpha_0) + \frac{1}{2} (3\alpha_0 - 1) \cos^2(\theta)& \\
|
|
- \beta_0 \sin^2(\theta) \cos(2 \phi)
|
|
- \gamma_0 \sqrt{2} \sin(2 \theta) \cos(\phi)& \Bigg]
|
|
\end{align*}
|
|
|
|
where $\theta$ and $\phi$ are, respectively, the polar and azimuthal angles, and
|
|
|
|
$$
|
|
\alpha_0 = 0.65 \et \beta_0 = 0.06 \et \gamma_0 = -0.18
|
|
$$
|
|
|
|
To generate the points a *try and catch* method was employed:
|
|
|
|
1. generate a point $(\theta , \phi)$ uniformly on a unit sphere
|
|
2. generate a third value $Y$ uniformly in $[0, 1]$
|
|
3. if @eq:requirement is satisfied save the point
|
|
4. repeat from 1.
|
|
|
|
$$
|
|
Y < F_0(\theta, \phi)
|
|
$$ {#eq:requirement}
|
|
|
|
To increase the efficiency of the procedure, the $Y$ values were actually
|
|
generated in $[0, 0.2]$, since $\max F_0 = 0.179$ for the given parameters
|
|
(an other option would have been generating the numbers in the range $[0 ; 1]$,
|
|
but all the points with $Y_{\theta, \phi} > 0.179$ would have been rejected with
|
|
the same probability, namely 1).
|
|
|
|
While $\phi$ can simply be generated uniformly between 0 and $2 \pi$, for
|
|
$\theta$ one has to be more careful: by generating evenly distributed numbers
|
|
between 0 and $\pi$, the points turn out to cluster on the poles. The point
|
|
$(\theta,\phi)$ is *not* uniform on the sphere but is uniform on a cylindrical
|
|
projection of it. (The difference can be appreciated in @fig:compare).
|
|
|
|
The proper distribution of $\theta$ is obtained by applying a transform
|
|
to a variable $x$ with uniform distribution in $[0, 1)$,
|
|
easily generated by the GSL function `gsl_rng_uniform()`.
|
|
|
|
<div id="fig:compare">
|
|
{width=45%}
|
|
{width=45%}
|
|
|
|
{width=45%}
|
|
{width=45%}
|
|
|
|
Examples of samples. On the left, points with $\theta$ evenly distributed
|
|
between 0 and $\pi$; on the right, points with $\theta$ properly distributed.
|
|
</div>
|
|
|
|
The transformation can be found by imposing the angular PDF to be a constant:
|
|
|
|
\begin{align*}
|
|
\frac{d^2 P}{d \Omega^2} = \text{const} = \frac{1}{4 \pi}
|
|
\hspace{20pt} &\Longrightarrow \hspace{20pt}
|
|
d^2P = \frac{1}{4 \pi} \, d \Omega^2 = \frac{1}{4 \pi} \,
|
|
d \phi \sin(\theta) d \theta
|
|
\\
|
|
&\Longrightarrow \hspace{20pt}
|
|
\frac{dP}{d \theta} = \int_0^{2 \pi} d \phi \, \frac{1}{4 \pi} \sin(\theta)=
|
|
\frac{1}{4 \pi} \sin(\theta) 2 \pi = \frac{1}{2} \sin(\theta)
|
|
\end{align*}
|
|
\begin{align*}
|
|
\theta = \theta (x)
|
|
\hspace{20pt} &\Longrightarrow \hspace{20pt}
|
|
\frac{dP}{d \theta} = \frac{dP}{dx} \cdot \left| \frac{dx}{d \theta} \right|
|
|
= \left. \frac{dP}{dx} \middle/ \, \left| \frac{d \theta}{dx} \right| \right.
|
|
\\
|
|
&\Longrightarrow \hspace{20pt}
|
|
\frac{1}{2} \sin(\theta) = \left. 1 \middle/ \, \left|
|
|
\frac{d \theta}{dx} \right| \right.
|
|
\end{align*}
|
|
|
|
If $\theta$ is chosen to grew together with $x$, then the absolute value can be
|
|
omitted:
|
|
|
|
\begin{align*}
|
|
\frac{d \theta}{dx} = \frac{2}{\sin(\theta)}
|
|
\hspace{20pt} &\Longrightarrow \hspace{20pt}
|
|
d \theta \sin(\theta) = 2 dx
|
|
\\
|
|
&\Longrightarrow \hspace{20pt}
|
|
- \cos (\theta') |_{0}^{\theta} = 2[x(\theta) - x(0)] = 2(x - 0) = 2x
|
|
\\
|
|
&\Longrightarrow \hspace{20pt}
|
|
- \cos(\theta) + 1 = 2x
|
|
\\
|
|
&\Longrightarrow \hspace{20pt}
|
|
\theta = \text{acos} (1 -2x)
|
|
\end{align*}
|
|
|
|
Since 1 is not in $[0 , 1)$, $\theta = \pi$ can't be generated. Being it just
|
|
a single point, the effect of this omission is negligible.
|
|
|
|
\vspace{30pt}
|
|
|
|
|
|
## Parameter estimation
|
|
|
|
The sample must now be used to estimate the parameters $\alpha_0$,
|
|
$\beta_0$ and $\gamma_0$ of the angular distribution $F_0$.
|
|
|
|
|
|
### Maximum Likelihood method
|
|
|
|
The Maximum Likelihood method (MLM) is based on the observation that the best
|
|
estimate {$\bar{\alpha}$, $\bar{\beta}$, $\bar{\gamma}$} of the parameters
|
|
{$\alpha$, $\beta$, $\gamma$} should maximize the Likelihood function $L$
|
|
defined in @eq:Like0, where the index $i$ runs over the sample and $F$ is the
|
|
function $F_0$ with free parameters $\alpha$, $\beta$ and $\gamma$:
|
|
|
|
$$
|
|
L = \prod_i F(\theta_i, \phi_i) = \prod_i F_i
|
|
$$ {#eq:Like0}
|
|
|
|
The reason is that the points were generated according to the probability
|
|
distribution function $F$ with parameters $\alpha_0$, $\beta_0$ and $\gamma_0$,
|
|
thus the probability $P$ of sampling this special sample is, by definition, the
|
|
product of their probabilities $F_i$ and it must be maximum for {$\alpha$,
|
|
$\beta$, $\gamma$} = {$\alpha_0$, $\beta_0$, $\gamma_0$}. Thus, by viewing $F_i$
|
|
as a function of the parameters, by maximizing $P = L$ with respect
|
|
to $\alpha$, $\beta$ and $\gamma$, a good estimate should be found.
|
|
Instead of actually maximising $L$, the function $-\ln(L)$
|
|
was minimised, as minimisation methods as usually described in literature
|
|
and $\ln$ since it simplifies the math:
|
|
|
|
$$
|
|
L = \prod_i F_i \thus - \ln(L) = - \sum_{i} \ln(F_i)
|
|
$$ {#eq:Like}
|
|
|
|
|
|
The problem of multidimensional minimisation was tackled using the GSL library
|
|
`gsl_multimin.h`, which implements a variety of methods, including the
|
|
conjugate gradient Fletcher-Reeves algorithm, used in the solution, which requires
|
|
the implementation of $-\ln(L)$ and its gradient.
|
|
|
|
To minimise a function $f$, given an initial guess point $x_0$, the gradient
|
|
$\nabla f$ is used to find the initial direction $v_0$ (the steepest descent) along
|
|
which a line minimisation is performed: the minimum occurs where the
|
|
directional derivative along $v_0$ is zero:
|
|
$$
|
|
\frac{\partial f}{\partial v_0} = \nabla f \cdot v_0 = 0
|
|
$$
|
|
or, equivalently, at the point where the gradient is orthogonal to $v_0$.
|
|
After this first step the following procedure is iterated:
|
|
|
|
1. Find the steepest descent $v_n$.
|
|
2. Compute the *conjugate* direction: $w_n = v_n + \beta_n w_{n-1}$
|
|
3. Perform a line minimisation along $w_n$
|
|
4. Update the position $x_{n+1} = x_n + \alpha_n w_n$
|
|
|
|
where $\alpha_n$ is line minimum and $\beta_n$ is given by the Fletcher-Reeves
|
|
formula:
|
|
$$
|
|
\beta = \frac{\|\nabla f(x_n)\|^2}{\|\nabla f(x_{n-1})\|^2}
|
|
$$
|
|
|
|
The accuracy of each line minimisation is controlled the parameter `tol`,
|
|
meaning:
|
|
|
|
$$
|
|
w \cdot \nabla f < \text{tol} \, \|w\| \, \|\nabla f\|
|
|
$$
|
|
|
|
The minimisation is quite unstable and this forced the starting point to be
|
|
taken very close to the solution, namely:
|
|
|
|
$$
|
|
\alpha_{sp} = 0.79 \et
|
|
\beta_{sp} = 0.02 \et
|
|
\gamma_{sp} = -0.17
|
|
$$
|
|
|
|
The overall minimisation ends when the gradient module is smaller than $10^{-4}$.
|
|
The program took 25 of the above iterations to reach the result shown in @eq:Like_res.
|
|
|
|
The Cramér-Rao bound states that the covariance matrix of parameters estimated
|
|
by MLM is greater than the of the inverse of the Hessian matrix of $-\log(L)$
|
|
at the minimum. Thus, the Hessian matrix was computed
|
|
analytically and inverted by a Cholesky decomposition, which states that
|
|
every positive definite symmetric square matrix $H$ can be written as the
|
|
product of a lower triangular matrix $L$ and its transpose $L^T$:
|
|
|
|
$$
|
|
H = L \cdot L^T
|
|
$$
|
|
|
|
The Hessian matrix is, indeed, both symmetric and positive definite, when
|
|
computed in a minimum. Once decomposed the inverse is given by
|
|
|
|
$$
|
|
H^{-1} = (L \cdot L^T)^{-1} = (L^{-1})^T \cdot L^{-1}
|
|
$$
|
|
The inverse of $L$ is easily computed since it's a triangular matrix.
|
|
|
|
The GSL library provides two functions `gsl_linalg_cholesky_decomp()` and
|
|
`gsl_linalg_cholesky_invert()` to decompose and invert in-place a matrix.
|
|
The result is shown below:
|
|
|
|
$$
|
|
\bar{\alpha_L} = 0.6541 \et
|
|
\bar{\beta_L} = 0.06125 \et
|
|
\bar{\gamma_L} = -0.1763
|
|
$$ {#eq:Like_res}
|
|
|
|
$$
|
|
\text{cov}_L = \begin{pmatrix}
|
|
9.225 \cdot 10^{-6} & -1.785 \cdot 10^{-6} & 2.41 \cdot 10^{-7} \\
|
|
-1.785 \cdot 10^{-6} & 4.413 \cdot 10^{-6} & 1.58 \cdot 10^{-6} \\
|
|
2.41 \cdot 10^{-7} & 1.58 \cdot 10^{-6} & 2.602 \cdot 10^{-6}
|
|
\end{pmatrix}
|
|
$$ {#eq:Like_cov}
|
|
|
|
Hence:
|
|
\begin{align*}
|
|
\alpha_L &= 0.654 \pm 0.003 \\
|
|
\beta_L &= 0.0613 \pm 0.0021 \\
|
|
\gamma_L &= -0.1763 \pm 0.0016
|
|
\end{align*}
|
|
|
|
|
|
### Least squares estimation
|
|
|
|
Another method that can be used to estimate {$\alpha$, $\beta$, $\gamma$} are
|
|
the least squares (LSQ), which is a multidimensional minimisation
|
|
specialised to the case of sum of squared functions called residuals.
|
|
The residuals can be anything in general but are usually interpreted as the
|
|
difference between an expected (fit) and an observed value. The least squares
|
|
then correspond to the expected values that best fit the observations.
|
|
|
|
To apply the LSQ method, the data must be properly binned, meaning that each
|
|
bin should contain a significant number of events (say greater than four or
|
|
five): in this case, 30 bins for $\theta$ and 60 for $\phi$ turned out
|
|
to be satisfactory. The expected values were given as the product of $N$, $F$
|
|
computed in the geometric center of the bin and the solid angle enclosed by the
|
|
bin itself, namely:
|
|
|
|
$$
|
|
E_{\theta, \phi} = N F(\theta, \phi) \, \Delta \theta \, \Delta \phi \,
|
|
\sin{\theta}
|
|
$$
|
|
|
|
Once the data are binned, the number of events in each bin plays the role of the
|
|
observed value, while the expected one is given by the probability of finding
|
|
an event in that bin multiplied by the total number of points, $N$. Then, the
|
|
$\chi^2$, defined as follow, must be minimized with respect to the three
|
|
parameters to be estimated:
|
|
|
|
$$
|
|
\chi^2 = \sum_i f_i^2 = \|\vec f\|^2 \with f_i = \frac{O_i - E_i}{\sqrt{E_i}}
|
|
$$
|
|
|
|
where the sum runs over the bins, $O_i$ are the observed values and $E_i$ the
|
|
expected ones. Besides the more generic `gsl_multimin.h` library, GSL provides a
|
|
special one for least square minimisation, called `gsl_multifit_nlinear.h`. A
|
|
trust region method was used to minimize the $\chi^2$.
|
|
In such methods, the $\chi^2$, its gradient and its Hessian matrix are computed
|
|
in a series of points in the parameters space till the gradient and the matrix
|
|
reach given proper values, meaning that the minimum has been well approached,
|
|
where "well" is related to that values.
|
|
If {$x_1 \dots x_p$} are the parameters which must be estimated, first $\chi^2$
|
|
is computed in the starting point $\vec{x_k}$, then its value in a point
|
|
$\vec{x}_k + \vec{\delta}$ is modelled by a function $m_k (\vec{\delta})$ which
|
|
is the second order Taylor expansion around the point $\vec{x}_k$, namely:
|
|
|
|
$$
|
|
\chi^2 (\vec{x}_k + \vec{\delta}) \sim m_k (\vec{\delta}) =
|
|
\chi^2 (\vec{x}_k) + \nabla_k^T \vec{\delta} + \frac{1}{2}
|
|
\vec{\delta}^T H_k \vec{\delta}
|
|
$$
|
|
|
|
where $\nabla_k$ and $H_k$ are the gradient and the Hessian matrix of $\chi^2$
|
|
in the point $\vec{x}_k$ and the superscript $T$ stands for the transpose. The
|
|
initial problem is reduced the so called trust-region subproblem (TRS), which
|
|
is the minimisation of $m_k(\vec{\delta})$ in a region where
|
|
$m_k(\vec{\delta})$ should be a good approximation of $\chi^2 (\vec{x}_k
|
|
+ \vec{\delta})$, given by the condition:
|
|
$$
|
|
\|D_k \vec\delta\| < \Delta_k
|
|
$$
|
|
If $D_k = I$ the region is a hypersphere of radius $\Delta_k$
|
|
centered at $\vec{x}_k$ but can be deformed according to $D_k$. This is
|
|
necessary in the case one or more parameters have scale very different scales.
|
|
|
|
Given an initial point $x_0$, radius $\Delta_0$, scale matrix $D_0$
|
|
and step $\vec\delta_0$ the LSQ algorithm consist in the following iteration:
|
|
|
|
1. Construct the function $m_k(\vec\delta)$.
|
|
2. Solve the TRS for step $\vec\delta_k$.
|
|
3. Check whether the solution actually decreases $\chi^2$
|
|
1. If positive increase the trust region radius $\Delta_{k+1} =
|
|
\alpha\Delta_k$ and shift the position $\vec x_{k+1} = \vec x_k +
|
|
\vec \delta_k$.
|
|
2. If negative decrease the radius $\Delta_{k+1} = \Delta_k/\beta$.
|
|
4. Repeat
|
|
|
|
This method is advantageous compared to a general minimisation because
|
|
the TRS usually amounts to solving a linear system. There are many algorithms
|
|
to solve the problem, in this case the Levenberg-Marquardt was used. It is based
|
|
on a theorem that proves the existence of a number $\mu_k$ such that
|
|
\begin{align*}
|
|
\Delta_k \mu_k = \|D_k \vec\delta_k\| &&
|
|
(H_k + \mu_k D_k^TD_k) \vec\delta_k = -\nabla_k
|
|
\end{align*}
|
|
Using the approximation[^1] $H\approx J^TJ$, obtained by computing the Hessian
|
|
of the first-order Taylor expansion of $\chi^2$, $\vec\delta_k$ can
|
|
be found by solving the system
|
|
|
|
$$
|
|
\begin{cases}
|
|
J_k \vec{\delta_k} + \vec{f_k} = 0 \\
|
|
\sqrt{\mu_k} D_k \vec{\delta_k} = 0
|
|
\end{cases}
|
|
$$
|
|
|
|
[^1]: Here $J_{ij} = \partial f_i/\partial x_j$ is the Jacobian matrix of the
|
|
vector-valued function $\vec f(\vec x)$.
|
|
|
|
The algorithm terminates if on of the following condition are satisfied:
|
|
|
|
1. $|\delta_i| \leq \text{xtol} (|x_i| + \text{xtol})$ for every component
|
|
of $\vec \delta$.
|
|
2. $\max_i |\nabla_i \cdot \max(x_i, 1)| \leq \text{gtol}\cdot \max(\chi^2, 1)$
|
|
3. $\|\vec f(\vec x+ \vec\delta) - \vec f(\vec x)|| \leq \text{ftol}
|
|
\cdot \max(\|\vec f(\vec x)\|, 1)$
|
|
|
|
These tolerance have all been set to \SI{1e-8}{}. The program converged in 7
|
|
iterations giving the results below. The covariance of the parameters can again
|
|
been estimated through the Hessian matrix at the minimum. The following results
|
|
were obtained:
|
|
|
|
$$
|
|
\bar{\alpha_{\chi}} = 0.6537 \et
|
|
\bar{\beta_{\chi}} = 0.05972 \et
|
|
\bar{\gamma_{\chi}} = -0.1736
|
|
$$ {#eq:lsq-res}
|
|
|
|
$$
|
|
\text{cov}_{\chi} = \begin{pmatrix}
|
|
9.244 \cdot 10^{-6} & -1.66 \cdot 10^{-6} & 2.383 \cdot 10^{-7} \\
|
|
-1.66 \cdot 10^{-6} & 4.495 \cdot 10^{-6} & 1.505 \cdot 10^{-6} \\
|
|
2.383 \cdot 10^{-7} & 1.505 \cdot 10^{-6} & 2.681 \cdot 10^{-6}
|
|
\end{pmatrix}
|
|
$$
|
|
|
|
Hence:
|
|
\begin{align*}
|
|
&\alpha_{\chi} = 0.654 \pm 0.003 \\
|
|
&\beta_{\chi} = 0.0597 \pm 0.0021 \\
|
|
&\gamma_{\chi} = -0.1736 \pm 0.0016
|
|
\end{align*}
|
|
|
|
|
|
### Results compatibility
|
|
|
|
In order to compare the values $x_L$ and $x_{\chi}$ obtained from both methods
|
|
with the correct ones (namely {$\alpha_0$, $\beta_0$, $\gamma_0$}), the
|
|
following compatibility t-test was applied:
|
|
|
|
$$
|
|
p = 1 - \text{erf}\left(\frac{t}{\sqrt{2}}\right)\ \with
|
|
t = \frac{x_i - x_0}{\sqrt{\sigma_{x_i}^2 + \sigma_{x_0}^2}}
|
|
= \frac{x_i - x_0}{\sigma_{x_i}}
|
|
$$
|
|
|
|
where $i$ stands either for the MLM or LSQ parameters and $\sigma_{x_i}$ is the
|
|
uncertainty of the value $x_i$. At 95% confidence level, the values are
|
|
compatible if $p > 0.05$.
|
|
|
|
Likelihood results:
|
|
|
|
----------------------------
|
|
par $p$-value
|
|
------------ ---------------
|
|
$\alpha_L$ 0.18
|
|
|
|
$\beta_L$ 0.55
|
|
|
|
$\gamma_L$ 0.023
|
|
----------------------------
|
|
|
|
Table: Likelihood results compatibility.
|
|
|
|
$\chi^2$ results:
|
|
|
|
---------------------------------
|
|
par $p$-value
|
|
----------------- ---------------
|
|
$\alpha_{\chi}$ 0.22
|
|
|
|
$\beta_{\chi}$ 0.89
|
|
|
|
$\gamma_{\chi}$ 0.0001
|
|
---------------------------------
|
|
|
|
Table: $\chi^2$ results compatibility.
|
|
|
|
It can be concluded that only the third parameter, $\gamma$ is not compatible
|
|
with the expected one in both cases. An in-depth analysis of the algebraic
|
|
arrangement of $F$ would be required in order to justify this outcome.
|
|
|
|
\vspace{30pt}
|
|
|
|
## Isotropic hypothesis testing
|
|
|
|
What if the probability distribution function was isotropic? Could it be
|
|
compatible with the found results? If $F$ was isotropic, then $\alpha_I$,
|
|
$\beta_I$ and $\gamma_I$ would be $1/3$, 0, and 0 respectively, since this
|
|
gives $F_I = 1/{4 \pi}$. The t-test gives a $p$-value approximately zero for all
|
|
the three parameters, meaning that there is no compatibility at all with this
|
|
hypothesis.
|