Giù Marcer 3f6a54c456 notes: correct integrals spacing

2020-07-05 11:36:29 +02:00

14 KiB

Raw Blame History

Exercise 1

Random numbers following the Landau distribution

The Landau distribution is a probability density function which can be defined as follows:


  f(x) = \int \limits_{0}^{+ \infty} dt \, e^{-t \log(t) -xt} \sin (\pi t)

{width=50%}

The GNU Scientific Library (GSL) provides a number of functions for generating random variates following tens of probability distributions. Thus, the function for generating numbers from the Landau distribution, namely gsl_ran_landau(), was used.
For the purpose of visualizing the resulting sample, the data was put into an histogram and plotted with matplotlib. The result is shown in @fig:landau.

{#fig:landau}

Randomness testing of the generated sample

Kolmogorov-Smirnov test

In order to compare the sample with the Landau distribution, the Kolmogorov-Smirnov (KS) test was applied. This test can be used to statistically quantifies the distance between the cumulative distribution function of the Landau distribution and the one of the sample. The null hypothesis is that the sample was drawn from the reference distribution.
The KS statistic for a given cumulative distribution function F(x) is:


  D_N = \text{sup}_x |F_N(x) - F(x)|

where:

x runs over the sample,
F(x) is the Landau cumulative distribution function,
F_N(x) is the empirical cumulative distribution function of the sample.

If N numbers were generated, for every point x, F_N(x) is simply given by the number of points preceding the point (itself included) normalized by N, once the sample is sorted in ascending order.
F(x) was computed numerically from the Landau distribution with a maximum relative error of 10^{-6}, using the function gsl_integration_qagiu(), found in GSL.

Under the null hypothesis, the distribution of D_N is expected to asymptotically approach a Kolmogorov distribution:


  \sqrt{N}D_N \xrightarrow{N \rightarrow + \infty} K

where K is the Kolmogorov variable, with cumulative distribution function given by [@marsaglia03]:


  P(K \leqslant K_0) = 1 - p = \frac{\sqrt{2 \pi}}{K_0}
  \sum_{j = 1}^{+ \infty} e^{-(2j - 1)^2 \pi^2 / 8 K_0^2}

Plugging the observed value \sqrt{N}D_N in K_0, the $p$-value can be computed. At 95% confidence level (which is the probability of confirming the null hypothesis when correct) the compatibility with the Landau distribution cannot be disproved if p > α = 0.05.
To approximate the series, the convergence was accelerated using the Levin $u$-transform with the gsl_sum_levin_utrunc_accel() function. The algorithm terminates when the difference between two successive extrapolations reaches a minimum.

For N = 50000, the following results were obtained:

D = 0.004
p = 0.38

Hence, the data were reasonably sampled from a Landau distribution.

Note: Contrary to what one would expect, the \chi^2 test on a histogram is not very useful in this case. For the test to be significant, the data have to be binned such that at least several points fall in each bin. However, it can be seen in @fig:landau that many bins are empty both in the right and left side of the distribution, so it would be necessary to fit only the region where the points cluster or use very large bins in the others, making the \chi^2 test unpractical.

Parameters comparison

When a sample of points is generated in a given range, different tests can be applied in order to check whether they follow a given distribution or not. The idea which lies beneath most of them is to measure how far the parameters of the distribution are from the ones measured in the sample.
The same principle can be used to verify if the generated sample effectively follows the Landau distribution. Since it turns out to be a very pathological PDF, very few parameters can be easily checked: mode, median and full width at half maximum (FWHM).

Mode

The mode of a set of data values is defined as the value that appears most often, namely: it is the maximum of the PDF. Since there is no closed form for the mode of the Landau PDF, it was computed numerically by the Brent algorithm (gsl_min_fminimizer_brent in GSL), applied to -f with a relative tolerance of 10^{-7}, giving:


  \text{expected mode: } m_e = \num{-0.22278298 \pm 0.00000006}

This is a minimization algorithm that begins with a bounded region known to contain a minimum. The region is described by a lower bound x_\text{min} and an upper bound x_\text{max}, with an estimate of the location of the minimum x_e. The value of the function at x_e must be less than the value of the function at the ends of the interval, in order to guarantee that a minimum is contained somewhere within the interval:


  f(x_\text{min}) > f(x_e) < f(x_\text{max})

On each iteration the function is interpolated by a parabola passing though the points x_\text{min}, x_e, x_\text{max} and the minimum is computed as the vertex of the parabola. If this point is found to be inside the interval, it is taken as a guess for the true minimum; otherwise the method falls back to a g olden section (using the ratio (3 - \sqrt{5})/2 \approx 0.3819660 proven to be optimal) of the interval. The value of the function at this new point x' is calculated. In any case, if the new point is a better estimate of the minimum, namely if f(x') < f(x_e), then the current estimate of the minimum is updated.
The new point allows the size of the bounded interval to be reduced, by choosing the most compact set of points which satisfies the constraint $f(a) > f(x') < f(b)$ between f(x_\text{min}), f(x_\text{min}) and f(x_e). The interval is reduced until it encloses the true minimum to a desired tolerance. The error of the result is estimated by the length of the final interval.

On the other hand, to compute the mode of the sample, the half-sample mode (HSM) or Robertson-Cryer estimator was used. This estimator was chosen because makes no assumptions on the underlying distribution and is not computationally expensive. The HSM is obtained by iteratively identifying the half modal interval, which is the smallest interval containing half of the observation. Once the sample is reduced to less than three points the mode is computed as the average. The special case n=3 is dealt with by averaging the two closer points [@robertson74].

To obtain a better estimate of the mode and its error, the above procedure was bootstrapped. The original sample was treated as a population and used to build 100 other samples of the same size, by sampling with replacements. For each one of the new samples, the above statistic was computed. By simply taking the mean of these statistics the following estimate was obtained:


  \text{observed mode: } m_o = \num{-0.29 \pm 0.19}

In order to compare the values m_e and m_0, the following compatibility $t$-test was applied:


  p = 1 - \text{erf}\left(\frac{t}{\sqrt{2}}\right)\ \with
  t = \frac{|m_e - m_o|}{\sqrt{\sigma_e^2 + \sigma_o^2}}

where \sigma_e and \sigma_o are the absolute errors of m_e and $m_o$ respectively. At 95% confidence level, the values are compatible if p > 0.05. In this case:

t = 1.012
p = 0.311

Thus, the observed mode is compatible with the mode of the Landau distribution, although the result is quite imprecise.

Median

The median is a central tendency statistics that, unlike the mean, is not very sensitive to extreme values, albeit less indicative. For this reason is well suited as test statistic in a pathological case such as the Landau distribution.
The median of a probability distribution is defined as the value such that its cumulative probability is 1/2. In other words, the median partitions the probability in two (connected) halves. The median of a sample, once sorted, is given by its middle element if the sample size is odd, or the average of the two middle elements otherwise.

The expected median was derived from the quantile function (QDF) of the Landau distribution¹. Once this is know, the median is simply given by \text{QDF}(1/2). Since both the CDF and QDF have no known closed form, they must be computed numerically. The cumulative probability was computed by quadrature-based numerical integration of the PDF (gsl_integration_qagiu() function in GSL). The function calculate an approximation of the integral:


  I(x) = \int\limits_x^{+\infty} f(t)dt

The CDF is then given by p(x) = 1 - I(x). This was done to avoid the left tail of the distribution, where the integration can sometimes fail. The integral I was actually mapped beforehand onto (0, 1] by the change of variable t = x + (1-u)/u, because the integration routine works on definite integrals. The result should satisfy the following accuracy requirement:


  |\text{result} - I| \le \max(\varepsilon_\text{abs}, \varepsilon_\text{rel}I)

where the absolute and relative tolerances \varepsilon_\text{abs} and \varepsilon_\text{rel} were set to \num{1e-10} and \num{1e-6}, respectively.
As for the QDF, this was implemented by numerically inverting the CDF. This was done by solving the equation;


  p(x) = p_0

for x, given a probability value p_0, where p(x) is the CDF. The (unique) root of this equation was found by a root-finding routine (gsl_root_fsolver_brent in GSL) based on the Brent-Dekker method. The following condition was checked for convergence:


  |a - b| < \varepsilon_\text{abs} + \varepsilon_\text{rel} \min(|a|, |b|)

where a,b are the current interval bounds. The condition immediately gives an upper bound on the error of the root as \varepsilon = |a-b|. The tolerances here were set to 0 and \num{1e-3}.
The result of the numerical computation is:


  \text{expected median: } m_e = \num{1.3557804 \pm 0.0000091}

while the sample median, obtained again by bootstrapping, was found to be:


  \text{observed median: } m_e = \num{1.3605 \pm 0.0062}

As stated above, the median is less sensitive to extreme values with respect to the mode: this lead the result to be much more precise. Applying again the aforementioned $t$-test to this statistic:

t=0.761
p=0.446

Hence, the two values show a good agreement.

FWHM

For a unimodal distribution (having a single peak) this statistic is defined as the distance between the two points at which the PDF attains half the maximum value. For the Landau distribution, again, there is no analytic expression known, thus the FWHM was computed numerically as follows. First of all, some definitions must be given:


  f_{\text{max}} = f(m_e) \et \text{FWHM} = x_+ - x_- \with
  f(x_\pm) = \frac{f_\text{max}}{2}

The function derivative f'(x) was minimized using the same minimization method used for finding m_e. Once f_\text{max} was known, the equation:


  f(x) = \frac{f_\text{max}}{2}

was solved by performing the Brent-Dekker method in the ranges $[x_\text{min}, m_e]$ and [m_e, x_\text{max}], where x_\text{min} and x_\text{max} are the first and last sampled point respectively, once all the points are sorted in ascending order. This lead to the two solutions x_\pm.
With a relative tolerance of \num{1e-7}, the following result was obtained:


\text{expected FWHM: } w_e = \num{4.0186457 \pm 0.0000001}

\vspace{-1em}

On the other hand, obtaining a good estimate of the FWHM from a sample is much more difficult. In principle, it could be measured by binning the data and applying the definition to the discretized values, however this yields very poor results and depends on an completely arbitrary parameter: the bin width.
A more refined method to construct a nonparametric empirical PDF function from the sample is a kernel density estimation (KDE). This method consist in convolving the (ordered) data with a smooth symmetrical kernel: in this case a standard Gaussian function. Given a sample of values \{x_i\}_{i=1}^N, the empirical PDF is defined as:


  f_\varepsilon(x) = \frac{1}{N\varepsilon} \sum_{i = 1}^N
    \mathcal{N}\left(\frac{x-x_i}{\varepsilon}\right)

where \mathcal{N} is the kernel and the parameter \varepsilon, called bandwidth, controls the strength of the smoothing. This parameter can be determined in several ways. For simplicity, it was chosen to use Silverman's rule of thumb [@silverman86], which gives:


  \varepsilon = 0.63 \, S_N
    \left(\frac{d + 2}{4}N\right)^{-1/(d + 4)}

where the 0.63 factor was chosen to compensate for the distortion that systematically reduces the peaks height, which affects the estimation of the mode, and:

S_N is the sample standard deviation;
d is the number of dimensions, in this case d=1.

With the empirical density estimation at hand, the FWHM can be computed by the same numerical method described for the true PDF. Again this was bootstrapped to estimate the standard error giving:


  \text{observed FWHM: } w_o = \num{4.06 \pm 0.08}

Applying the $t$-test to these two values gives

t=0.495
p=0.620

which shows a very good agreement and proves the estimator is robust. For reference, the initial estimation based on an histogram gave a rather inadequate \si{4 \pm 2}.

This is neither necessary nor the easiest way: it was chosen simply because the quantile had been already implemented and was initially used for reverse sampling. ↩︎

14 KiB Raw Blame History Unescape Escape