<jupyter_start><jupyter_text>[Python for Probability, Statistics, and Machine Learning](https://www.springer.com/fr/book/9783319307152)<jupyter_code>from __future__ import division
%pylab inline<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>We considered Maximum Likelihood Estimation (MLE) and Maximum A-Posteriori(MAP) estimation and in each case we started out with a probability densityfunction of some kind and we further assumed that the samples were identicallydistributed and independent (iid). The idea behind robust statistics[[maronna2006robust]](maronna2006robust) is to construct estimators that can survive theweakening of either or both of these assumptions. More concretely, suppose youhave a model that works great except for a few outliers. The temptation is tojust ignore the outliers and proceed. Robust estimation methods provide adisciplined way to handle outliers without cherry-picking data that works foryour favored model. The Notion of LocationThe first notion we need is *location*, which is  a generalization of the ideaof *central value*. Typically, we just use an estimate of the mean for this,but we will see later why this could be a bad idea.  The general idea oflocation satisfies the following requirements Let $X$ be a random variable withdistribution $F$, and let $\theta(X)$ be some descriptive measure of $F$. Then$\theta(X)$ is said to be a measure of *location* if for any constants *a* and*b*, we have the following:$$\begin{equation}\theta(X+b)  = \theta(X) +b \label{_auto1} \tag{1}\end{equation}$$ $$\begin{equation} \\theta(-X)   = -\theta(X)  \label{_auto2} \tag{2}\end{equation}$$ $$\begin{equation} \X \ge 0 \Rightarrow \theta(X)  \ge 0  \label{_auto3} \tag{3}\end{equation}$$ $$\begin{equation} \\theta(a X) = a\theta(X)\label{_auto4} \tag{4}\end{equation}$$  The first condition is called *location equivariance* (or *shift-invariance* insignal processing lingo). The fourth condition is called *scale equivariance*,which means that the units that $X$ is measured in should not effect the valueof the  location estimator.  These requirements capture the intuition of*centrality* of a distribution, or where most of theprobability mass is located.For example, the sample mean estimator is $ \hat{\mu}=\frac{1}{n}\sum X_i $. The firstrequirement is obviously satisfied as $ \hat{\mu}=\frac{1}{n}\sum (X_i+b) = b +\frac{1}{n}\sum X_i =b+\hat{\mu}$. Let us consider the second requirement:$\hat{\mu}=\frac{1}{n}\sum -X_i = -\hat{\mu}$. Finally, the last requirement issatisfied with $ \hat{\mu}=\frac{1}{n}\sum a X_i =a \hat{\mu}$. Robust Estimation and ContaminationNow that we have the generalized location of centrality embodied in the*location* parameter, what can we do with it?  Previously, we assumed that our sampleswere all identically distributed. The key idea is that the samples might beactually coming from a *single* distribution that is contaminated by another nearbydistribution, as in the following: $$F(X) = \epsilon G(X) + (1-\epsilon)H(X)$$  where $ \epsilon $ randomly toggles between zero and one. This meansthat our data samples $\lbrace X_i \rbrace$ actually derived from two separatedistributions, $ G(X) $ and $ H(X) $. We just don't know how they are mixedtogether. What we really want  is an estimator  that captures the location of $G(X) $ in the face of random intermittent contamination by $ H(X)$.  Forexample, it may be that this contamination is responsible for the outliers in amodel that otherwise works well with the dominant $F$ distribution. It can geteven worse than that because we don't know that there is only one contaminating$H(X)$ distribution out there. There may be a whole family of distributionsthat are contaminating $G(X)$. This means that whatever estimators we constructhave to be derived from a more generalized family of distributions instead offrom a single distribution,  as the maximum-likelihood method assumes.  This iswhat makes robust estimation so difficult --- it has to deal with *spaces* offunction distributions instead of parameters from a particular probabilitydistribution. Generalized Maximum Likelihood EstimatorsM-estimators are generalized maximum likelihood estimators. Recall that formaximum likelihood, we want to maximize the likelihood function as in thefollowing: $$L_{\mu}(x_i) = \prod f_0(x_i-\mu)$$  and then to find the estimator $\hat{\mu}$ so that $$\hat{\mu} = \arg \max_{\mu} L_{\mu}(x_i)$$   So far, everything is the same as our usual maximum-likelihoodderivation except for the fact that we don't assume a specific $f_0$ as thedistribution of the $\lbrace X_i\rbrace$. Making the definition of $$\rho = -\log f_0$$  we obtain the more convenient form of the likelihood product and theoptimal $\hat{\mu}$ as $$\hat{\mu} = \arg \min_{\mu} \sum \rho(x_i-\mu)$$  If $\rho$ is differentiable, then differentiating  this with respectto $\mu$ gives $$\begin{equation}\sum \psi(x_i-\hat{\mu}) = 0 \label{eq:muhat} \tag{5}\end{equation}$$  with $\psi = \rho^\prime$, the first derivative of $\rho$ , and for technical reasons we will assume that$\psi$ is increasing. So far, it looks like we just pushed some definitionsaround, but the key idea is we want to consider general $\rho$ functions thatmay not be maximum likelihood estimators for *any* distribution. Thus, ourfocus is now on uncovering the nature of $\hat{\mu}$. Distribution of M-estimatesFor a given distribution $F$, we define $\mu_0=\mu(F)$ as the solution to thefollowing $$\mathbb{E}_F(\psi(x-\mu_0))= 0$$  It is technical to show, but it turns out that $\hat{\mu} \sim\mathcal{N}(\mu_0,\frac{v}{n})$ with $$v = \frac{\mathbb{E}_F(\psi(x-\mu_0)^2)}{(\mathbb{E}_F(\psi^\prime(x-\mu_0)))^2}$$  Thus, we can say that $\hat{\mu}$ is asymptotically normal with asymptoticvalue $\mu_0$ and asymptotic variance $v$. This leads to the efficiency ratiowhich is defined as  the following: $$\texttt{Eff}(\hat{\mu})= \frac{v_0}{v}$$  where $v_0$ is the asymptotic variance of the MLE and measures hownear $\hat{\mu}$ is to the optimum. In other words, this provides a sense ofhow much outlier contamination costs in terms of samples. For example, if fortwo estimates with asymptotic variances $v_1$ and $v_2$, we have $v_1=3v_2$,then first estimate requires three times as many observations to obtain thesame variance as the second. Furthermore, for the sample mean (i.e.,$\hat{\mu}=\frac{1}{n} \sum X_i$) with $F=\mathcal{N}$, we have $\rho=x^2/2$and $\psi=x$ and also $\psi'=1$. Thus, we have $v=\mathbb{V}(x)$.Alternatively, using the sample median as the estimator for the location, wehave $v=1/(4 f(\mu_0)^2)$.  Thus, if we have $F=\mathcal{N}(0,1)$, for thesample median, we obtain $v={2\pi}/{4} \approx 1.571$. This means that thesample median takes approximately 1.6 times as many samples to obtain the samevariance for the location as the sample mean. The sample median is far more immune to the effects of outliers than the sample mean, so this gives a sense of how much this robustness costs in samples.** M-Estimates as Weighted Means.** One way to think about M-estimates is aweighted means. Operationally, thismeans that we want weight functions that can circumscribe theinfluence of the individual data points, but, when taken as a whole,still provide good estimated parameters. Most of the time, we have $\psi(0)=0$ and $\psi'(0)$ exists sothat $\psi$ is approximately linear at the origin. Using the followingdefinition: $$W(x)  =  \begin{cases}                \psi(x)/x & \text{if} \: x \neq 0 \\\                \psi'(x)  & \text{if} \: x =0             \end{cases}$$  We can write our Equation ref{eq:muhat} as follows: $$\begin{equation}\sum W(x_i-\hat{\mu})(x_i-\hat{\mu}) = 0 \label{eq:Wmuhat} \tag{6}\end{equation}$$  Solving this for $\hat{\mu} $ yields the following, $$\hat{\mu} = \frac{\sum w_{i} x_i}{\sum w_{i}}$$  where $w_{i}=W(x_i-\hat{\mu})$. This is not practically usefulbecause the $w_i$ contains $\hat{\mu}$, which is what we are trying to solvefor. The question that remains is how to pick the $\psi$ functions. This isstill an open question, but the Huber functions are a well-studied choice. Huber FunctionsThe family of Huber functions is defined by the following: $$\rho_k(x ) = \begin{cases}                x^2         & \mbox{if }  |x|\leq k \\\                2 k |x|-k^2 & \mbox{if }  |x| >  k                \end{cases}$$  with corresponding derivatives $2\psi_k(x)$ with $$\psi_k(x ) = \begin{cases}                x              & \mbox{if } \: |x| \leq k \\\                \text{sgn}(x)k & \mbox{if } \: |x| >  k                \end{cases}$$  where the limiting cases $k \rightarrow \infty$ and $k \rightarrow 0$correspond to the mean and median, respectively. To see this, take$\psi_{\infty} = x$ and therefore $W(x) = 1$ and thus the defining Equationref{eq:Wmuhat} results in $$\sum_{i=1}^{n} (x_i-\hat{\mu}) = 0$$  and then solving this leads to $\hat{\mu} = \frac{1}{n}\sum x_i$.Note that choosing $k=0$ leads to  the sample median, but that is not sostraightforward to solve for. Nonetheless, Huber functions provide a wayto move between two extremes of estimators for location (namely, the mean vs. the median) with a tunable parameter $k$. The $W$ function corresponding to Huber's $\psi$ is the following: $$W_k(x) = \min\Big{\lbrace} 1, \frac{k}{|x|} \Big{\rbrace}$$  [Figure](fig:Robust_Statistics_0001) shows the Huber weightfunction for $k=2$ with some sample points. The idea is that the computedlocation, $\hat{\mu}$ is computed from Equation ref{eq:Wmuhat} to lie somewherein the middle of the weight function so that those terms (i.e., *insiders*)have their values fully reflected in the location estimate. The black circlesare the *outliers* that have their values attenuated by the weight function sothat only a fraction of their presence is represented in the location estimate. -->This shows the Huber weight function, $W_2(x)$ and some cartoon data points that are insiders or outsiders as far as the robust location estimate is concerned. Breakdown PointSo far, our discussion of robustness has been very abstract.  A more concreteconcept of robustness comes from the breakdown point.  In the simplest terms,the breakdown point describes what happens when a single data point in anestimator is changed in the most damaging way possible. For example, suppose wehave the sample mean, $\hat{\mu}=\sum x_i/n$, and we take one of the $x_i$points to be infinite. What happens to this estimator? It also goes infinite.This means that the breakdown point of the estimator is 0%. On the other hand,the median has a breakdown point of 50%, meaning that half of the data forcomputing the median could go infinite without affecting the median value. The medianis a *rank* statistic that cares more about the relative ranking of the datathan the values of the data, which explains its robustness.The simpliest but still formal way to express the breakdown point is totake $n$ data points, $\mathcal{D} = \lbrace (x_i,y_i) \rbrace$. Suppose $T$is a regression estimator that yields a vector of regression coefficients,$\boldsymbol{\theta}$, $$T(\mathcal{D}) = \boldsymbol{\theta}$$  Likewise, consider all possible corrupted samples of the data$\mathcal{D}^\prime$. The maximum *bias* caused by this contamination isthe following: $$\texttt{bias}_{m} = \sup_{\mathcal{D}^\prime} \Vert T(\mathcal{D^\prime})-T(\mathcal{D}) \Vert$$  where the $\sup$ sweeps over all possible sets of $m$ contaminated samples.Using this, the breakdown point is defined as the following: $$\epsilon_m = \min \Big\lbrace \frac{m}{n} \colon \texttt{bias}_{m} \rightarrow \infty \Big\rbrace$$  For example, in our least-squares regression, even one point atinfinity causes an infinite $T$. Thus, for least-squares regression,$\epsilon_m=1/n$. In the limit $n \rightarrow \infty$, we have $\epsilon_m\rightarrow 0$. Estimating ScaleIn robust statistics, the concept of *scale* refers to a measure of thedispersion of the data. Usually, we use theestimated standard deviation for this, but this has a terrible breakdown point.Even more troubling, in order to get a good estimate of location, we have toeither somehow know the scale ahead of time, or jointly estimate it. None ofthese methods have easy-to-compute closed form solutions and must be computednumerically.The most popular method for estimating scale is the *median absolute deviation* $$\texttt{MAD} = \texttt{Med} (\vert \mathbf{x} - \texttt{Med}(\mathbf{x})\vert)$$  In words, take the median of the data $\mathbf{x}$ andthen subtract that median from the data itself, and then take the median of theabsolute value of the result. Another good dispersion estimate is the *interquartile range*, $$\texttt{IQR} = x_{(n-m+1)} - x_{(n)}$$  where $m= [n/4]$. The $x_{(n)}$ notation means the $n^{th}$ dataelement after the data have been sorted. Thus, in this notation,$\texttt{max}(\mathbf{x})=x_{(n)}$. In the case where $x \sim\mathcal{N}(\mu,\sigma^2)$, then $\texttt{MAD}$ and $\texttt{IQR}$ are constantmultiples of $\sigma$ such that the normalized $\texttt{MAD}$ is the following, $$\texttt{MADN}(x) = \frac{\texttt{MAD} }{0.675}$$   The number comes from the inverse CDF of the normal distributioncorresponding to the $0.75$ level. Given the complexity of thecalculations, *jointly* estimating both location and scale is a purelynumerical matter. Fortunately, the Statsmodels module has many of theseready to use. Let's create some contaminated data in the following code,<jupyter_code>import statsmodels.api as sm
from scipy import stats
data=np.hstack([stats.norm(10,1).rvs(10),stats.norm(0,1).rvs(100)])<jupyter_output><empty_output><jupyter_text>These data correspond to our model of contamination that we startedthis section with. As shown in the  histogram in [Figure](fig:Robust_Statistics_0002), there are two normal distributions, onecentered neatly at zero, representing the majority of the samples, and anothercoming less regularly from the normal distribution on the right. Notice thatthe group of infrequent samples on the right separates the mean and medianestimates (vertical dotted and dashed lines).  In the absence of thecontaminating distribution on the right, the standard deviation for this datashould be close to one. However, the usual non-robust estimate for standarddeviation (`np.std`) comes out to approximately three.  Using the$\texttt{MADN}$ estimator (`sm.robust.scale.mad(data)`) we obtain approximately1.25. Thus, the robust estimate of dispersion is less moved by the presence ofthe  contaminating distribution. -->Histogram of sample data. Notice that the group of infrequent samples on the right separates the mean and median estimates indicated by the vertical lines.The generalized maximum likelihood M-estimation extends to jointscale and location estimation using Huber functions. For example,<jupyter_code>huber = sm.robust.scale.Huber()
loc,scl=huber(data)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10532647196417831, 'fraction_numerical': 0.008821760342177371, 'mean_word_length': 4.598204264870931, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 16, 'https://': 1, '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': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '556501', 'license': 'MIT', 'n_tokens_mistral': 4263, 'n_tokens_neox': 3983, 'n_words': 1935, 'path': 'chapters/statistics/notebooks/Robust_Statistics.ipynb', 'repo_name': 'leopoldodonnell/Python-for-Probability-Statistics-and-Machine-Learning', 'size': '14963'}