Common Statistical Tools¶

Author: Yannick Copin y.copin@ipnl.in2p3.fr $ \renewcommand{\d}{\mathrm{d}} \renewcommand{\v}[1]{\boldsymbol{#1}} % Vector \newcommand{\m}[1]{\mathsf{#1}} % Matrix \DeclareMathOperator{\Var}{Var} \DeclareMathOperator{\Cov}{Cov} \DeclareMathOperator{\med}{med} $

Abstract: A restricted number of statistical tools are presented, with emphasis on details sometimes overlooked in literature. This may be too simple, or too technical, or whatever: comments and suggestions are welcome.

Sample statistics¶

We define two datasets:

Case without intrinsic dispersion: $\v{x} = (x_1,\ldots, x_N)$ is a sample of $N$ measurements of the same null quantity, with (usually estimated) gaussian measurement uncertainties $\d x_{i}^2$ following a log-normal distribution (here, $\log\d x^2$ follows a $\mathcal{N}(0,1)$ distribution).

Case with intrinsic dispersion: $\v{y} = (y_1,\ldots, y_N)$ is an intrinsically dispersed $\mathcal{N}(\mu=0,\sigma^2)$ sample, measured with the same gaussian measurement uncertainties $\d x_{i}$ as above. Unfortunately, the intrinsic dispersion $\sigma$ is usually not known a priori.

Maximum-likelihood estimators¶

Sample location¶

Arithmetic average: $$\begin{align} \hat{\mu} &= \frac{1}{N}\sum_i x_i \\ \Var(\hat{\mu}) &= \frac{\hat{\sigma}^{2}}{N} \end{align}$$

Weighted ("optimal") average: $$\begin{align} \hat{\mu}_w &= \frac{\sum_i w_i x_i}{\sum_i w_i} \quad\text{with}\quad w_i = \frac{1}{\d x_{i}^{2} + \sigma^2} \\ \Var(\hat{\mu}_w) &= \frac{1}{\sum_i w_i} \end{align}$$

Attention:

• In the optimally weighted mean, $\sigma_i^2 = 1/w_{i}$ is supposed to be the total variance of normally-distributed $x_i$, therefore including both measurement error $\d x_i$ and intrinsic distribution dispersion $\sigma$ if any: $\sigma_i^2 = \d x_i^2 + \sigma^2$.

• If the weighting is not the true inverse variance $1/\sigma_i^2$, the weighted average is optimal (with minimal variance) but biased.

• The reduced $\chi^2$ around estimated mean, defined as $\chi_\nu^2 = 1/(N-1) \times \sum_i (x_i - \mu_w)^2 / \sigma_i^2$, should not be significantly different from one for $\sigma_i$ compatible to the observed dispersion. If this is not the case, $\sigma_i$ has probably been misestimated, either because it ignored the intrinsic dispersion $\sigma$ or because the measurement uncertainties $\d x_i$ are misestimated, or both.

Sample scale¶

Unweighted: The supposedly unbiased sample variance wrt. unknown mean $\mu$ (but see discussion in Oliphant, 2006) is: $$\begin{align} \hat{\sigma}^2 &= \frac{1}{N-1} \sum_i (x_i - \hat{\mu})^2 \\ \Var(\hat{\sigma}^2) &= \frac{1}{N}\left(m^4 - \frac{N-3}{N-1}\sigma^4\right) \\ &\simeq \frac{1}{N} \times 2\hat{\sigma}^4 \quad\text{in the normal approximation} \\ \text{therefore}\quad \Var(\hat{\sigma}) &\simeq \frac{1}{2N}\hat{\sigma}^2 \\ \text{Ahn & Fessler (2003):}\quad \Var(\hat{\sigma}) &= \frac{1}{2(N-1)}\hat{\sigma}^2 \quad\text{pour}\quad n > 10. \end{align}$$

Weighted: $$ \hat{\sigma}_w^2 = \frac{\sum w_i}{\left(\sum w_i\right)^2 - \sum w_i^2} \times \sum_i w_i (x_i - \hat{\mu}_w)^2. $$

Weighted RMS is usually defined as:

This strictly corresponds to the biased weighted variance. If using inverse-variance weighting $w_i = 1/\sigma_i^2$, one has $\mathrm{WRMS}^2 = \chi^2 \times \sigma_{\mu_w}^2$. Under the assumption that the intrinsic scatter is null ($\sigma\equiv 0$), the WRMS can be used as a proxy for the rescaled uncertainty on the weighted mean when the errors $\d x_i$ were misestimated by a multiplicative factor: $$ \sigma_{\mu_w}^{*2} = \chi^2_\nu \times \sigma_{\mu_w}^2 = \mathrm{WRMS}^2/(N-1). $$

Note: the weighted mean $\hat{\mu}_{w}$ is unaffected by a scaling error on the errors $\d x_{i}$.

Bayesian estimates¶

They are studied in Oliphant, 2006, and implemented in scipy.stats.mvsdist (bayesian estimate of the sample mean, variance and standard deviation PDFs) or scipy.stats.bayes_mvs (bayesian confidence intervals).

Robust estimators¶

Median-based estimators are robust (50% of the points can be wrong), but have a low accuracy. More pertinent estimators can be defined.

Location¶

The variance of the median $m$ for large samples and normal distributions is: $$ \Var(m) = \frac{\pi}{2}\frac{s^{2}}{n}, $$ where $s$ is an (presumably robust) estimate of the scale.

Scale¶

Median Absolute Deviation: $$ \mathrm{MAD} = \med(| x - \med(x) |) $$

Inter-Quartile Range: this is the 25%-trimmed range, i.e. the difference between the 75th- and the 25th-percentiles of a sample.

Both estimators have a optimal breakpoint of 0.5: up to half the points of the sample can be arbitrarily large before the estimator gives an arbitrarily large estimate. (By comparison, the sample standard deviation has a breakpoint of 0.)

For a normal distribution, these are related to standard deviation $\sigma$ by: $$\begin{align} \hat{\sigma} &= \frac{1}{\Phi^{-1}(3/4)}\mathrm{MAD} \simeq 1.4826\times\mathrm{MAD} \\ % 1/SS.norm.ppf(0.75) = 1.4826 &= \frac{1}{2\Phi^{-1}(3/4)} \mathrm{IQR} \simeq 0.7413\times\mathrm{IQR} % 2*sqrt(2)*scipy.special.erfinv(0.5) = 2*SS.norm.ppf(0.75) = 1.349 \end{align}$$ where $\Phi$ is the standard normal cumulative distribution function (scipy.stats.norm.cdf), and its inverse $\Phi^{-1}$ is known as the normal quantile function, or probit function (scipy.stats.norm.ppf). Note that MAD is biased for low number of points (see Croux & Rousseeuw, 1992).

Unfortunately, $\mathrm{MAD}$ is only 37% as efficient as the sample standard deviation.

More robust estimators of scale with maximal breakdown point are presented in Croux and Rousseeuw (1992).

Absolute Pairwise Differences: defined in Rousseeuw and Croux (1993) $$\begin{align} S_n &= 1.1926 \med_i \left( \med_j(\left| x_i - x_j \right|) \right), \\ Q_n &= c_n \text{first quartile of} \left( \left| x_i - x_j \right| : i < j \right), \end{align}$$

Covariance matters¶

$\chi^2$ minimisation¶

Consider $\v{y}$ a column-vector of $N$ measurements $y_i$, each of which assumed to be Gaussian distributed with a mean $F_i(\v{\theta})$ (hereafter the model $\v{F}(\v{\theta})$, depending on $M$ parameters $\theta_k$). The maximum-likelyhood estimate $\hat{\v{\theta}}$ of set of parameters $\v{\theta}$ is the one minimizing the $\chi^2$ defined by: $$ \chi^2(\v{\theta}) = (\v{y} - \v{F}(\v{\theta}))^T \cdot \m{V}^{-1} \cdot (\v{y} - \v{F}(\v{\theta})), $$ where $\m{V}$ is the (supposedly known) covariance matrix between the measurements: $V_{ij} = \Cov(y_i, y_j)$.

In that case, $\chi^2(\v{\theta})$ follows a $\chi^2$-distribution with $k=N-M$ degrees of freedom. In particular:

• The goodness-of-fit (level of agreement between the measurements and the model) can be assessed from the minimum value $\chi^2_\min = \chi^2(\hat{\v{\theta}})$.
• $\chi^2 = \chi^2_\min + 1$ defines the $1\,\sigma$-isocontours around the maximum-likelyhood estimate $\hat{\v{\theta}}$. Using a Taylor-series quadratic development of $\chi^2(\v{\theta})$ around its minimum, one gets an approximation of the best-fit paramater covariance matrix:

$$ \Cov(\hat{\v{\theta}}) \simeq 2\,\m{H}^{-1} $$

where $\m{H}$ the Hessian matrix of the $\chi^2$ estimated at $\hat{\v{\theta}}$:

\begin{align} H_{ij} &= \left.\frac{\partial^2\chi^2}{\partial\theta_i\partial\theta_j}\right|_{\v{\theta}=\hat{\v{\theta}}} \\ &= 2 \left.\frac{\partial\v{F}}{\partial\theta_i}\right|_{\hat{\v{\theta}}}^T \cdot \m{V}^{-1} \cdot \left.\frac{\partial\v{F}}{\partial\theta_j}\right|_{\hat{\v{\theta}}}. \end{align}

If the measurements are uncorrelated, the covariance matrix $\m{V}$ is purely diagonal, with $V_{ii} = \sigma_i^2$ the (supposedly known) variance of measurement $y_i$, and the previous equations can be rewritten:

\begin{align} \chi^2(\v{\theta}) &= \sum_{i=1}^N \frac{(y_i - F_i(\v{\theta}))^2}{\sigma_i^2}, \\ H_{ij} &= 2 \sum_{k=1}^N \frac{1}{\sigma_k^2} \frac{\partial F_k}{\partial\theta_i}\frac{\partial F_k}{\partial\theta_j}. \end{align}

Example: weighted mean¶

We define two datasets:

• an uncorrelated repeated measurement sample $\v{x}$ ($N=100$) of a null-quantity with log-normal measurement uncertainties,
• a correlated sample $\v{y}$, derived from $\v{x}$ with an exponential correlation length of $\tau = 5$ px and same measurement uncertainties.

The goodness-of-fit with respect to constant value $\mu = 0$ is given by:

In this particular case, the fact of disregarding the covariance for the correlated dataset leads to a spurious 2.4-sigma tension with the true average $\mu = 0$.

The objective is to derive the maximum-likelihood estimate of the average values, the model is simply the 1-parameter function $\v{F} = (\mu, \ldots, \mu)$. In presence of a covariance matrix $\m{V}$, one has then $\sigma_{\hat{\mu}}^2 = \left(\sum_{ij} (\m{V}^{-1})_{ij}\right)^{-1}$. For uncorrelated measurements, this gives $\sigma_{\hat{\mu}}^2 = \left(\sum_{i} 1/\sigma_i^2\right)^{-1}$ (uncertainty on variance weighted mean, see above).

The effect of neglecting the covariance in correlated data is even stronger on ML-estimate of mean: the estimated value is more than 7-sigma away from the true value, with an unrealistically good $\chi^2_\min$.

Propagation of uncertainties¶

Consider a set of $N$ quantities (measurements, parameter estimates, etc.) $\v{y} = (y_1, \ldots y_N)$ associated to a covariance matrix $V_{ij} = \Cov(y_i, y_j)$ (estimated from measurement or adjustment procedures), and a set of $M$ functions $\v{f} = (f_1, \ldots, f_M)$. The covariance matrix among the propagated values $\v{f}(\v{y})$ can be estimated from a Taylor expansion around $\v{y}$. To first order, one has: $$ \Cov(\v{f}(\v{y})) \simeq \m{J}\cdot\m{V}\cdot{\m{J}}^T $$ where the matrix of derivative ("jacobian") is given by: $J_{ij} = \left.\frac{\partial f_i}{\partial y_j}\right|_{\v{y}}$.

This 1st-order expression is exact if $\v{f}$ is linear: $\v{f}(\v{y}) = \m{J}\cdot\v{y}$. If this is not the case, one has to make sure the functions are reasonably linear over the domain $\v{y} \pm \v{\sigma}_{\v{y}}$, or expand the Taylor series to higher orders.

The 1st-order expression be rewritten in a non-matricial form: $$ \Cov(f_i, f_j) \simeq \sum_{k,l = 1}^{N} \frac{\partial f_i}{\partial y_k}\frac{\partial f_j}{\partial y_l} V_{kl}. $$ In the case of uncorrelated input quantities, $V_{ij} = \sigma_i^2 \delta_{ij}$ and one recovers the traditional uncertainty propagation expression: $$ \sigma_{f_i}^2 = \sum_{k=1}^{N} \left(\frac{\partial f_i}{\partial y_k}\right)^2 \sigma_k^2. $$

Example: synthetic photometry¶

Consider a spectrum made of $N$ flux measurements $\v{y} = (y_1, \ldots y_N)$ associated to a covariance matrix $V_{ij} = \Cov(y_i, y_j)$. One wants to compute the synthetic photometry (i.e. integrated flux) in bands $A$ and $B$, \begin{align} A &= \sum_i \alpha_i y_i \\ B &= \sum_i \beta_i y_i \end{align} where the fixed values of $\alpha$ and $\beta$ coefficients define the filter bandpass, as well as resulting variances and covariances.

In this simple linear case, fluctuations $\delta A$ are directly related to flux fluctuations $\delta y_i$: $$ \delta A = \sum_i \alpha_i \delta y_i $$ and therefore \begin{align} \sigma_A^2 &= \langle \delta^2 A \rangle \\ &= \sum_{i,j} \alpha_i \alpha_j \langle \delta y_i \delta y_j \rangle \\ &= \v{\alpha} \cdot \m{V} \cdot \v{\alpha}^T \\ \Cov(A, B) &= \langle \delta A \delta B \rangle \\ &= \sum_{i,j} \alpha_i \beta_j \langle \delta y_i \delta y_j \rangle \\ &= \v{\alpha} \cdot \m{V} \cdot \v{\beta}^T \end{align}

Regression coefficient¶

Correlation coefficients measure the relationship between two datasets $\{x_{i}\}$ and $\{y_{j}\}$.

• Pearson correlation coefficient
• Spearman correlation coefficient

Plot tests¶

Histograms¶

Model adjustment¶

Linear regression¶

Least-squares¶

Chi2¶

Cash statistics¶

Likelihood¶

Resampling techniques¶

Pull¶

• Everything you always wanted to know about pulls (Demortier & Lyons, 2002)

References¶

Besides standard statistics handbooks,

• The Review of Particle Physics - Statistics
• NIST/SEMATECH e-Handbook of Statistical Methods
• Hogg et al. (2010)
• Andrae (2010)