<jupyter_start><jupyter_text>Introduction and examples of robust directional and multivariate median estimators<NAME> and <NAME>Astrophysics Group, Cavendish Laboratory, University of Cambridge<jupyter_code>import numpy as np
import scipy
import seaborn as sns
from scipy import stats
from sklearn.covariance import MinCovDet

from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import mark_inset, zoomed_inset_axes
from scipy.stats import shapiro
from scipy.stats.mstats import gmean as geometric_mean

from robstat.robstat import circ_mean_dev, geometric_median, mardia_median, \
mv_median, mv_normality, mv_outlier, tukey_median
from robstat.utils import round_down, round_up
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Directional data  Mardia MedianMardia median given by the angle that the circular mean deviation:$$ d(\tilde\theta) = \pi - \frac{1}{n} \sum_{i=1}^{n} \left| \pi - \left| \theta_i - \tilde\theta \right| \right| $$where $\tilde\theta$ is the estimate of the preferred direction, and it is used as a measure of dispersion.The Mardia median occasionally leads to a non-unique estimate of the circular median since there can sometimes be two or more diameters that divide the data equally and have the same circular mean deviation.Weighted Mardia median:$$ d(\tilde\theta) = \pi - \frac{1}{\sum \eta_i} \sum_{i=1}^{n} \eta_i \left| \pi - \left| \theta_i - \tilde\theta \right| \right| $$<jupyter_code># create sample angular data
angle_center = 0.5
np.random.seed(0)
angles = angle_center + np.random.normal(loc=0, scale=0.2, size=200)

x_coords = np.cos(angles)
y_coords = np.sin(angles)
fig, ax = plt.subplots(figsize=(6, 6), dpi=100)

circ_rad = 1
lim_rng = circ_rad * 1.25

ax.set(xlim=(-lim_rng, lim_rng), ylim = (-lim_rng, lim_rng))

a_circle = plt.Circle((0, 0), 1, fill=False, color='blue', alpha=0.5)
ax.add_artist(a_circle)

ax.axhline(0, color='grey', ls='--', lw=0.5)
ax.axvline(0, color='grey', ls='--', lw=0.5)
ax.plot(x_coords, y_coords, 'o', color='orange', markersize=4, alpha=0.5)

plt.show()
# circular mean deviation
np.array([circ_mean_dev(angles, angle_center)]).item()
# checking uniqueness of Mardia Median

x = mardia_median(angles, init_guess=angle_center-0.3)
y = mardia_median(angles, init_guess=angle_center+0.3)

print('Mardia Medians found:')
print(x)
print(y, '\n')

print('Circular mean: \n{}\n'.format(stats.circmean(angles, nan_policy='omit')))

cmd_x = circ_mean_dev(angles, x).item()
cmd_y = circ_mean_dev(angles, y).item()
print('Circular mean deviations for the Mardia Medians (check):')
print(cmd_x)
print(cmd_y, '\n')

if ~np.isclose(x, y):
    print('The Mardia Median in this example is not unique')
    take_mean = 'Mean'
else:
    take_mean = ''
fig, ax = plt.subplots(figsize=(6, 6), dpi=100, subplot_kw={'projection': 'polar'})

ax.set(rlim=(0, lim_rng))
ax.set_rlabel_position(-90)

# ax.plot(np.linspace(0, 2*np.pi, 1000), np.ones(1000), color='blue', linestyle='-', alpha=0.5)
ax.plot(angles, np.ones_like(angles), 'o', color='orange', markersize=4, alpha=0.5, zorder=0)
ax.scatter(np.mean([x, y]), 1, zorder=1, label='{} Mardia median'.format(take_mean))

ax.set_xlim(round_down(angles.min(), 1), round_up(angles.max(), 1))

plt.legend(loc=0)
plt.show()<jupyter_output><empty_output><jupyter_text>Multivariate data  Geometric MedianThe geometric median is defined as the value of the argument $y$ that minimizes the sum of Euclidian distances between $y$ and all points $x_i$:$$ \underset{y \in \mathbb{R}^d}{\mathrm{arg\,min}} \sum_{i=1}^n || x_i - y ||  $$Properties & asides: - The geometric median has a breakdown point of 0.5: up to half of the sample data may be arbitrarily corrupted, and the median of the samples will still provide a robust estimator for the location of the uncorrupted data - The geometric median is unique whenever the points are not collinear - Weiszfeld's algorithm for faster computation We can also weight the geometric median (c.f. the Weber problem):$$ \underset{y \in \mathbb{R}^d}{\mathrm{arg\,min}} \sum_{i=1}^n \eta_i || x_i - y ||  $$  Tukey MedianTukey (1975) proposed the halfspace depth as a tool to visually describe multivariate datasets.For a finite set of data points $\mathcal{X}_n = \{x_1, ..., x_n\}$ in $\mathbb{R}^d$, the Tukey depth, or halfspace depth of any point $y \in \mathbb{R}^d$ determines how central the point is inside the data cloud; it is defined as the minimal number of data points in any closed halfspace determined by a hyperplane through $y$:$$ d_T(y; \mathcal{X}_n) = \underset{|| u || = 1}{\min} \ \{i \in \{1, \dots, n\} : u^{\intercal} x_i \geq u^{\intercal} y \} $$The set of all points with depth $\geq \kappa$ is called the Tukey $\kappa^{\mathrm{th}}$ depth region $\mathcal{D}(\kappa)$:$$\mathcal{D}(\kappa) = \{ x \in \mathbb{R}: d_T(x) \geq \kappa \} $$The halfspace depth regions form a sequence of closed convex polyhedra, with each polyhedron included in $\mathrm{conv} (x1, \dots, x_n)$, and making it compact. Tukey regions are also nested: they shrink with increasing $\kappa$. An empirical distribution is fully characterized by its Tukey regions.Donoho \& Gasko (1992) noted that this could be the basis of a definition for a multivariate median: the points that maximizes the Tukey depth, known as the Tukey median, is given by:$$ \mathcal{T} = \underset{y \in \mathbb{R}^d}{\arg \max} \; d_T (y; \mathcal{X}_n) $$<jupyter_code>points = np.random.random(500).reshape(-1, 2)*2
points = np.concatenate((points, np.random.random(50).reshape(-1, 2)+5))
points[:, 0] += 2
points[:, 1] += 2
points_c = points[:, 0] + points[:, 1]*1j
sample_gmean = geometric_mean(points_c)
sample_gmed = geometric_median(points_c, weights=None)
sample_tmed = tukey_median(points_c)['barycenter']
bad_med = lambda x : np.nanmedian(x.real) + np.nanmedian(x.imag)*1j
sample_bmed = bad_med(points_c)
med_ests = list(zip([sample_gmean, sample_gmed, sample_tmed, sample_bmed], \
               ['Geometric Mean', 'Geometric Median', 'Tukey Median', 'Separate Median'], \
               ['ro', 'co', 'yo', 'bo']))
fig, ax = plt.subplots(figsize=(6, 6), dpi=100)

ax.scatter(points[:, 0], points[:, 1], alpha=0.5)
for i, med_est in enumerate(med_ests):
    ax.plot(med_est[0].real, med_est[0].imag, med_est[2], label=med_est[1])
ax.set_xlim(1.5, 9)
ax.set_ylim(1.5, 9)

# zoomed in sub region of the original image
axins = zoomed_inset_axes(ax, zoom=6, loc=4)
axins.scatter(points[:, 0], points[:, 1], alpha=0.5)
for i, med_est in enumerate(med_ests):
    axins.plot(med_est[0].real, med_est[0].imag, med_est[2])

x1 = round_down(np.min([med_est[0].real for med_est in med_ests]), 1)
x2 = round_up(np.max([med_est[0].real for med_est in med_ests]), 1)
y1 = round_down(np.min([med_est[0].imag for med_est in med_ests]), 1)
y2 = round_up(np.max([med_est[0].imag for med_est in med_ests]), 1)
axins.set_xlim(x1, x2)
axins.set_ylim(y1, y2)

axins.tick_params(axis='x', direction='in', pad=-15)
mark_inset(ax, axins, loc1=1, loc2=3, fc='none', ec='0.5')

ax.legend(loc='upper left')
plt.show()
# effect of outliers in data

print('Geometric median')
print(geometric_median(points[:250], weights=None))
print(geometric_median(points, weights=None), '\n')

print('Tukey median')
print(tukey_median(points[:250], weights=None)['barycenter'])
print(tukey_median(points, weights=None)['barycenter'], '\n')

print('Bad median')
print(bad_med(points[:250]))
print(bad_med(points))<jupyter_output><empty_output><jupyter_text>Other location depth notions - Simplicial depth - Oja depth - Projection depth - Spatial depth  See e.g. https://www.csun.edu/~ctoth/Handbook/chap58.pdf, https://cran.r-project.org/web/packages/depth/depth.pdf<jupyter_code>for mvm in ['Tukey', 'Oja', 'Liu', 'Spatial', 'CWmed']:
    print('{:7s}: {}'.format(mvm, mv_median(points, method=mvm)))<jupyter_output><empty_output><jupyter_text>Normality tests  Shapiro-Wilk The Shapiro–Wilk test (1965) tests the null hypothesis that a sample $x_1, \dots, x_n$ comes from a normal distribution. The test statistic is given by$$ W = \frac{\left( \sum_{i=1}^n a_i x_{\left( i \right)} \right)^2}{\sum_{i=1}^n \left(x_i - \bar{x} \right)^2} $$where $x_{\left( i \right)}$ are the ordered sample values (i.e. the the $i$th-smallest number in the sample), $\bar{x}$ is the sample mean, and $a_i$ are constants generated from the means, with $\vec{a} = \left( a_1, \dots, a_n \right)$ given by$$ \vec{a} = \frac{m^{\intercal}V^{-1}}{C} $$with vector norm $C$$$ C = \lVert V^{-1} m \rVert = \left(m^{\intercal} V^{-1} V^{-1} m \right)^{1/2} $$and $\vec{m} = \left( m_1, \dots, m_n \right)^\intercal$ made of the expected values of the order statistics of *iid* random variables sampled from the standard normal distribution; finally, $V$ is the covariance matrix of those normal order statisticsSmall values of $W$ are evidence of departure from normality.The null-hypothesis of this test is that the sample is normally distributed. Thus, if the $p$ value is less than the chosen threshold level (usually 5%), then the null hypothesis is rejected and there is evidence that the sample data are not normally distributed.<jupyter_code>shapiro_x = shapiro(points[:, 0])
shapiro_y = shapiro(points[:, 1])
print(shapiro_x)
print(shapiro_y)<jupyter_output><empty_output><jupyter_text>Henze-Zirkler The Henze-Zirkler test (1990) statistic is based on a non-negative functional that measures the distance between two distribution functions, the hypothesized function (which is the multivariate normal) and the observed function.The HZ test statistic is given by$$ HZ = \frac{1}{n} \sum^n_{i=1} \sum^n_{j=1} e^{-\frac{\beta^2}{2} D_{ij}} - 2 (1 + \beta^2)^{-\frac{p}{2}} \sum^n_{j=1} e^{-\frac{\beta^2}{2 (1+\beta^2)} D_{i}} + n (1 + 2 \beta^2)^{-\frac{p}{2}} $$where$$ p = \ \; \mathrm{variables} $$$$ \beta = \frac{1}{\sqrt{2}} \left( \frac{n (2p + 1)}{4} \right)^{\frac{1}{p+4}} $$$$ D_{ij} = (x_i - x_j)' S^{-1} (x_i - x_j) $$$$ D_{i} = (x_i - \bar{x})' S^{-1} (x_i - \bar{x}) $$with $D_{i}$ giving the squared Mahalanobis distance of $i^{\mathrm{th}}$ observation to the centroid and $D_{ij}$ the Mahalanobis distance between the $i^{\mathrm{th}}$ and $j^{\mathrm{th}}$ observations, as $S$ is the covariance matrix. Mahalanobis distance is a metric that calculates how far each observation is to the center of joint distribution.If the data is multivariate normal, $HZ$ is approximately log-normally distributed.<jupyter_code>MVN_res = mv_normality(points)
print(MVN_res)<jupyter_output><empty_output><jupyter_text>Outlier detection Outliers are the common reason for violating the assumption of normality. A multivariate outlier method based on robust Mahalanobis distances is presented here.Mahalanobis distances are used in the Henze-Zirkler test above. They are given by:$$ \mathrm{MD}(x_i) = \left( (x_i - \hat{\mu}) \hat{\Sigma}^{-1} (x_i - \hat{\mu})^\intercal \right)^{1/2} $$ where $\hat{\mu}$ is the sample multivariate mean and $\hat{\Sigma}$ is the sample covariance matrix.Of course, the mean and covariance will be heavily influenced by the presence of outliers; obtaining good robust estimators of $\hat{\mu}$ and $\hat{\Sigma}$ are necessary to measure the "outlyingness" of data points and to have a proper distance-based outlier detection procedure. Therefore, in line with this package, we modify the above to get robust Mahalanobis distances:$$ \mathrm{rMD}(x_i) = \left( (x_i - \hat{\mu}_r) \hat{\Sigma}_r^{-1} (x_i - \hat{\mu}_r)^\intercal \right)^{1/2} $$where where $\hat{\mu}_r$ and $\hat{\Sigma}_r$ are robust estimators of centrality and covariance matrix.In practice, the most frequently used estimator is the minimum covariance determinant $\mathrm{MCD}$ estimator, which is based on the computation of the ellipsoid with the smallest volume or with the smallest covariance determinant that would encompass at least half of the data points.This method consists on determining the subset $J$ of observations of size $h$ that minimizes the determinant of the sample covariance matrix, computed from only these $h$ points. The choice of $h$ determines the robustness of the estimator; it is a compromise between robustness and efficiency. Once this subset of size $h$ is found, it is possible to estimate the centrality and the covariance matrix based only upon that subset.$$ J = \left\{ h : | \hat{\Sigma}_J | < | \hat{\Sigma}_K | \; \forall K \; \mathrm{s.t.} \; \K = h \right\}$$$$ \hat{\mu}_{\mathrm{MCD}} = \frac{1}{h} \sum_{i \in J} x_i $$$$ \hat{\Sigma}_{\mathrm{MCD}} = \frac{1}{h} \sum_{i \in J} (x_i - \hat{\mu}_{\mathrm{MCD}}) (x_i - \hat{\mu}_{\mathrm{MCD}})^\intercal $$where $||$ denotes the determinant of the matrix, and $\K$ denotes the cardinality of the subset $K$.The steps for finding outliers are as follows: - Compute robust Mahalanobis distances $\mathrm{rMD}(x_i)$ - Compute the 97.5% quantile $Q$ of the chi-square distribution $\chi^2_{d; \; 0.975}$ ($d = \ \mathrm{d.o.f.}$) - Declare $\mathrm{rMD}(x_i) > Q$ as possible outliers The quantile threshold (97.5% by default) can of course be changed as desired.We note that the threshold for this method can be modified such that the threshold is adjusted to the sample size; an adjusted quantiled $AQ$ is used (that replaces $Q$, Filzmoser 2005) in the the second step. Proceeding with the $AQ$ option generally improves the false classification rates, while maintaining the same correct classification rates (Cabana 2017). The adjusted threshold is computed by comparing the theoretical cumulative $\chi_d^2$ distribution function and the empirical cumulative distribution function of the squared robust distance samples,  and finding the the supremum of the difference between the two tails of these distributions.<jupyter_code>mvo_res = mv_outlier(points, method='quan')
mvo_res.head()
fig, axes = plt.subplots(ncols=2, figsize=(10, 6), sharey=True)
sns.scatterplot(x=points[:, 0], y=points[:, 1], hue=mvo_res['RS Mahalanobis Distance'], \
                ax=axes[0])
sns.scatterplot(x=points[:, 0], y=points[:, 1], hue=mvo_res['Outlier'], ax=axes[1])
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Alternative Python-only method<jupyter_code>robust_cov = MinCovDet().fit(points)
chi2_quantile = 0.975
chi2_thresh = stats.chi2.ppf(chi2_quantile, points.shape[1])
outliers = np.where(robust_cov.mahalanobis(points) > chi2_thresh)[0]
fig, ax = plt.subplots(figsize=(6, 6))

inliers = np.delete(points, outliers, axis=0)
sns.scatterplot(x=inliers[:, 0], y=inliers[:, 1], ax=ax)
sns.scatterplot(x=points[outliers, 0], y=points[outliers, 1], color='orange', ax=ax, \
                label='Outliers')

# Create meshgrid of feature values
xx, yy = np.meshgrid(np.linspace(plt.xlim()[0], plt.xlim()[1], 100),
                     np.linspace(plt.ylim()[0], plt.ylim()[1], 100))
zz = np.c_[xx.ravel(), yy.ravel()]

# Calculate the MCD based Mahalanobis distances
mahal_robust_cov = robust_cov.mahalanobis(zz)
mahal_robust_cov = mahal_robust_cov.reshape(xx.shape)
robust_contour = ax.contour(xx, yy, np.sqrt(mahal_robust_cov), cmap=plt.cm.YlOrBr_r, \
                            linestyles='--')
ax.clabel(robust_contour, robust_contour.levels, inline=True, fontsize=10, fmt=r'%r')
ax.annotate('Robust Mahalanobis Distance', xy=(0.55, 0.05), xycoords='axes fraction', \
            bbox=dict(boxstyle='round', facecolor='white'), size=10, color='orange')

plt.legend()
plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14962512926577043, 'fraction_numerical': 0.019131334022750777, 'mean_word_length': 5.1425168717745136, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 37, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '58121', 'license': 'MIT', 'n_tokens_mistral': 5586, 'n_tokens_neox': 5229, 'n_words': 1660, 'path': 'notebooks/example_notebook.ipynb', 'repo_name': 'matyasmolnar/robstat', 'size': '15486'}