Estimation of a multivariate von Mises distribution for contaminated torus data

• Free University of Bozen-Bolzano, Italy [giulia.bertagnolli@unibz.it]

• University G. Fortunato, Benevento, Italy [l.greco@unifortunato.eu]

• University of Trento, Trento, Italy [claudio.agostinelli@unitn.i]

Keywords: Weighted likelihood – Circular data – ICORS – Robustness

1 Contamination in circular data

Observations of wind directions, arm movements and other angles measurements, but also time data over the 24 hours (0:00 corresponding to $0^{\circ }$ and 24:00 to ${360}^{\circ }$), due to their periodicity, can be thought and modelled as circular data, that is, points on the surface of a torus. Its particular topology requires careful adaptation of classical linear statistical methods for the analysis of observations on its surface, as well as specific families of distributions to model them. The von Mises distribution is a well-known distribution for circular data. The $p-$variate density of the von Mises sine distribution is

$$m(𝜽;𝝁,𝜿,\mathrm{\Lambda })=C^{-1}(𝜿,\mathrm{\Lambda })\exp [{𝜿}^{\top }\cos (𝜽-𝝁)-\frac{1}{2}\sin {(𝜽-𝝁)}^{\top }\mathrm{\Lambda }\sin (𝜽-𝝁)]$$ (1)

where $𝜽\in {[0,2\pi )}^p$ is a $p-$variate circular random variable, $𝝁\in {[0,2\pi )}^p$ is the location vector, $𝜿={(k_1,\mathrm{\dots },k_p)}^{\top }$, with $k_j>0$, is the concentration vector, ${\mathrm{\Lambda }}_{ij}={\lambda }_{ij}={\lambda }_{ji}$ and ${\mathrm{\Lambda }}_{jj}=0$ for $i,j=1,\mathrm{\dots },p$ [Mardia et al., 2008]. The normalising constant cannot be expressed in closed form for $p>2$, so we make the usual assumption that the concentration values $k_j$ are sufficiently large, and that the matrix $\mathrm{\Sigma }$, defined as ${\mathrm{\Sigma }}_{ii}^{-1}=k_i,{\mathrm{\Sigma }}_{ij}^{-1}=-{\lambda }_{ij}$, is positive definite. In this case, the density in eq. (1) can be approximated by the concentrated multivariate sine distribution (CMS) of Mardia et al. [2012]

$$m^c(𝜽;𝝁,𝜿,\mathrm{\Lambda })={\left(2\pi \right)}^{-p/2}|\mathrm{\Sigma }{|}^{-1/2}\exp [-{𝜿}^{\top }\cos (𝜽-𝝁)-\frac{1}{2}\sin {(𝜽-𝝁)}^{\top }\mathrm{\Lambda }\sin (𝜽-𝝁)].$$

Figure 1: Sample of size $n=250$ from a bivariate von Mises distribution with $𝝁=(0,0)$, ${𝜿}_0=(10,20)$ and ${\lambda }_0=15$ with the addition of 50 outliers. (Left) Data are displayed as points on ${𝕋}^2$. (Right) Fitted density contours obtained using approximate maximum likelihood estimates superimposed.

Of course, circular data are not exempt from contaminations, which can very badly affect MLE-based inference. Figure 1 shows, in the left panel, a sample of 250 observations from a bivariate von Mises distribution, which has been contaminated with 50 outliers. The contours superimposed on the scatter plot in the right plane highlight how the contaminated estimates differ from the original parameter values. The weighted likelihood methodology [Lindsay, 1994, Markatou et al., 1998] has been proved both efficient and effective for robust estimation, and has been recently applied also to univariate [Agostinelli, 2007] and multivariate [Greco et al., 2021, Saraceno et al., 2021, Agostinelli et al., 2024] circular data. However, the study of the $p-$variate von Mises distribution is still missing. In [Bertagnolli et al., 2024], we fill this gap.

1.1 Weighted likelihood estimating equations

Intuitively, robustness is achieved by down-weighting, in the score estimating equation, those observations that are unlikely w.r.t. the assumed model. Let us assume $𝝁=\mathrm{𝟎}$, and denote by $𝝉=(𝜿,\mathrm{\Lambda })$ the parameters of interest, by $m^c(\theta ;𝝉)$ the assumed CMS model, by $u(𝜽;𝝉)={\nabla }_{𝝉}\log m^c(𝜽;𝝉)$ the corresponding score, by ${\widehat{F}}_n$ the empirical distribution function of a random sample, and by $f^{*}$ a non-parametric kernel estimate of the true unknown density $f$ of $𝜽$. The agreement of the $i-$th observation with the assumed model is quantified through its Pearson residual:

$$\delta ({𝜽}_{𝒊})=\delta ({𝜽}_i;𝝉,{\widehat{F}}_n)=\frac{f^{*}({𝜽}_i)}{m^{c*}({𝜽}_i;𝝉)}-1\in [-1,\mathrm{\infty }).$$ (2)

$\delta ({𝜽}_i)=0$ indicates perfect agreement with the model, while values away from zero label ${𝜽}_i$ as an out-lier, $\delta (𝒚)\to \mathrm{\infty }$, or in-lier, $\delta (𝒚)\to -1$. To guarantee that the Pearson residuals converge to 0 with probability 1 when the model is correct, the model is also smoothed through the same kernel of $f^{*}$, that is $m^{c*}$. The WL estimating equation (WLEE) then reads

$$\sum _{i=1}^{n}w(\delta ({𝜽}_i);𝝉,{\widehat{F}}_n)u({𝜽}_i;𝝉)=0,\text{ with }w(\delta )=\frac{A(\delta )+1}{\delta +1},$$ (3)

where $A(\delta )$ is a so-called residual adjustment function (RAF), whose aim is to down-weight observations with large or negative residuals, but which also provides a direct link with disparity minimisation and the robustness properties of the derived estimator [Lindsay, 1994]. We use a symmetric chi-square (SCHI) RAF.

2 Results

Figure 2: 8TIM protein data. Left: fitted density contours from WLE. Right: fitted density contours from MLE. The point estimate is denoted by the symbol X.

The proposed technique proves satisfactory both on synthetic and real data examples. Figure 2 shows results of estimation based on WL (left) and ML (right) for the well-known protein 8TIM bivariate data, concerning $n=490$ backbone torsion angle pairs $(\varphi ,\psi )$ [Chakraborty and Wong, 2021]. This dataset has multiple clusters, which are not detected by MLE. The WL methodology, instead, gives strong indications of the presence of such sub-structures, which can be unveiled by further inspection of the multiple roots of the WLEE.