What would a robust Bayesian model for estimating the scale of a roughly normal distribution be?

Question

There exists a number of robust estimators of scale. A notable example is the median absolute deviation which relates to the standard deviation as $\sigma = \mathrm{MAD}\cdot1.4826$. In a Bayesian framework there exist a number of ways to robustly estimate the location of a roughly normal distribution (say a Normal contaminated by outliers), for example, one could assume the data is distributed as a t distribution or Laplace distribution. Now my question:

What would a Bayesian model for measuring the scale of a roughly normal distribution in a robust way be, robust in the same sense as the MAD or similar robust estimators?

As is the case with MAD, it would be neat if the Bayesian model could approach the SD of a normal distribution in the case when the distribution of the data actually is normally distributed.

edit 1:

A typical example of a model that is robust against contamination/outliers when assuming the data $y_i$ is roughly normal is using a t distribution like:

$$y_i \sim \mathrm{t}(m, s,\nu)$$

Where $m$ is the mean, $s$ is the scale, and $\nu$ is the degree-of-freedom. With suitable priors on $m, s$ and $\nu$, $m$ will be an estimate of the mean of $y_i$ that will be robust against outliers. However, $s$ will not be a consistent estimate of the SD of $y_i$ as $s$ depends on $\nu$. For example, if $\nu$ would be fixed to 4.0 and the model above would be fitted to a huge number of samples from a $\mathrm{Norm}(\mu=0,\sigma=1)$ distribution then $s$ would be around 0.82. What I'm looking for is a model which is robust, like the t model, but for the SD instead of (or in addition to) the mean.

edit 2:

Here follows a coded example in R and JAGS of how the t-model mentioned above is more robust with respect to the mean.

# generating some contaminated data
y <- c( rnorm(100, mean=10, sd=10), 
        rnorm(10, mean=100, sd= 100))

#### A "standard" normal model ####
model_string <- "model{
  for(i in 1:length(y)) {
    y[i] ~ dnorm(mu, inv_sigma2)
  }

  mu ~ dnorm(0, 0.00001)
  inv_sigma2 ~ dgamma(0.0001, 0.0001)
  sigma <- 1 / sqrt(inv_sigma2)
}"

model <- jags.model(textConnection(model_string), list(y = y))
mcmc_samples <- coda.samples(model, "mu", n.iter=10000)
summary(mcmc_samples)

### The quantiles of the posterior of mu
##  2.5%   25%   50%   75% 97.5% 
##   9.8  14.3  16.8  19.2  24.1 

#### A (more) robust t-model ####
library(rjags)
model_string <- "model{
  for(i in 1:length(y)) {
    y[i] ~ dt(mu, inv_s2, nu)
  }

  mu ~ dnorm(0, 0.00001)
  inv_s2 ~ dgamma(0.0001,0.0001)
  s <- 1 / sqrt(inv_s2)
  nu ~ dexp(1/30) 
}"

model <- jags.model(textConnection(model_string), list(y = y))
mcmc_samples <- coda.samples(model, "mu", n.iter=1000)
summary(mcmc_samples)

### The quantiles of the posterior of mu
## 2.5%   25%   50%   75% 97.5% 
##8.03  9.35  9.99 10.71 12.14 

Answer 1

Bayesian inference in a T noise model with an appropriate prior will give a robust estimate of location and scale. The precise conditions that the likelihood and prior need to satisfy are given in the paper Bayesian robustness modelling of location and scale parameters by Andrade and O'Hagan (2011). The estimates are robust in the sense that a single observation cannot make the estimates arbitrarily large, as demonstrated in figure 2 of the paper.

When the data is normally distributed, the SD of the fitted T distribution (for fixed $\nu$) does not match the SD of the generating distribution. But this is easy to fix. Let $\sigma$ be the standard deviation of the generating distribution and let $s$ be the standard deviation of the fitted T distribution. If the data is scaled by 2, then from the form of the likelihood we know that $s$ must scale by 2. This implies that $s = \sigma f(\nu)$ for some fixed function $f$. This function can be computed numerically by simulation from a standard normal. Here is the code to do this:

library(stats)
library(stats4)
y = rnorm(100000, mean=0,sd=1)
nu = 4
nLL = function(s) -sum(stats::dt(y/s,nu,log=TRUE)-log(s))
fit = mle(nLL, start=list(s=1), method="Brent", lower=0.5, upper=2)
# the variance of a standard T is nu/(nu-2)
print(coef(fit)*sqrt(nu/(nu-2)))

For example, at $\nu=4$ I get $f(\nu)=1.18$. The desired estimator is then $\hat{\sigma} = s/f(\nu)$.

Answer 2

As you are asking a question about a very precise problem (robust estimation), I will offer you an equally precise answer. First, however, I will begin be trying to dispel an unwarranted assumption. It is not true that there is a robust bayesian estimate of location (there are bayesian estimators of locations but as I illustrate below they are not robust and, apparently, even the simplest robust estimator of location is not bayesian) . In my opinion, the reasons for the absence of overlap between the 'bayesian' and 'robust' paradigm in the location case goes a long way in explaining why there also are no estimators of scatter that are both robust and bayesian .

With suitable priors on $m, s$ and $\nu$, $m$ will be an estimate of the mean of $y_i$ that will be robust against outliers.

Actually, no. The resulting estimates will only be robust in a very weak sense of the word robust. However, when we say that the median is robust to outliers we mean the word robust in a much stronger sense. That is, in robust statistics, the robustness of the median refers to the property that if you compute the median on a data-set of observations drawn from a uni-modal, continuous model and then replace less than half of these observations by arbitrary values, the value of the median computed on the contaminated data is close to the value you would have had had you computed it on the original (uncontaminated) data-set. Then, it is easy to show that the estimation strategy you propose in the paragraph I quoted above is definitely not robust in the sense of how the word is typically understood for the median.

I'm wholly unfamiliar with Bayesian analysis. However, I was wondering what is wrong with the following strategy as it seems simple, effective and yet has not been considered in the other answers. The prior is that the good part of the data is drawn from a symmetric distribution $F$ and that the rate of contamination is less than half. Then, a simple strategy would be to:

1. compute the median/mad of your dataset. Then compute: $$z_i=\frac{|x_i-\mbox{med}(x)|}{\mbox{mad}(x)}$$
2. exclude the observations for which $z_i>q_{\alpha}(z|x\sim F)$ (this is the $\alpha$ quantile of the distribution of $z$ when $x\sim F$). This quantity is avalaible for many choice of $F$ and can be bootstrapped for the others.
3. Run a (usual, non-robust) Bayesian analysis on the non-rejected observations.

EDIT:

Thanks to the OP for providing a self contained R code to conduct a bonna fide bayesian analysis of the problem.

the code below compares the the bayesian approach suggested by the O.P. to it's alternative from the robust statistics literature (e.g. the fitting method proposed by Gauss for the case where the data may contain as much as $n/2-2$ outliers and the distribution of the good part of the data is Gaussian).

central part of the data is $\mathcal{N}(1000,1)$:

n<-100
set.seed(123)
y<-rnorm(n,1000,1)

Add some amount of contaminants:

y[1:30]<-y[1:30]/100-1000 
w<-rep(0,n)
w[1:30]<-1

the index w takes value 1 for the outliers. I begin with the approach suggested by the O.P.:

library("rjags")
model_string<-"model{
  for(i in 1:length(y)){
    y[i]~dt(mu,inv_s2,nu)
  }
  mu~dnorm(0,0.00001)
  inv_s2~dgamma(0.0001,0.0001)
  s<-1/sqrt(inv_s2)
  nu~dexp(1/30) 
}"

model<-jags.model(textConnection(model_string),list(y=y))
mcmc_samples<-coda.samples(model,"mu",n.iter=1000)
print(summary(mcmc_samples)$statistics[1:2])
summary(mcmc_samples)

I get:

     Mean        SD 
384.2283  97.0445 

and:

2. Quantiles for each variable:

 2.5%   25%   50%   75% 97.5% 
184.6 324.3 384.7 448.4 577.7 

(quiet far thus from the target values)

For the robust method,

z<-abs(y-median(y))/mad(y)
th<-max(abs(rnorm(length(y))))
print(c(mean(y[which(z<=th)]),sd(y[which(z<=th)])))

one gets:

 1000.149 0.8827613

(very close to the target values)

The second result is much closer to the real values. But it gets worst. If we classify as outliers those observations for which the estimated $z$-score is larger than th (remember that the prior is that $F$ is Gaussian) then the bayesian approach finds that all the observations are outliers (the robust procedure, in contrast, flags all and only the outliers as such). This also implies that if you were to run a usual (non-robust) bayesian analysis on the data not classified as outliers by the robust procedure, you should do fine (e.g. fulfil the objectives stated in your question).
This is just an example, but it's actually fairly straightforward to show that (and it can done formally, see for example, in chapter 2 of [1]) the parameters of a student $t$ distribution fitted to contaminated data cannot be depended upon to reveal the outliers.

• [1]Ricardo A. Maronna, Douglas R. Martin, Victor J. Yohai (2006). Robust Statistics: Theory and Methods (Wiley Series in Probability and Statistics).
• Huber, P. J. (1981). Robust Statistics. New York: John Wiley and Sons.

Answer 3

In bayesian analysis using the inverse Gamma distribution as a prior for the precision (the inverse of the variance) is a common choice. Or the inverse Wishart distribution for multivariate models. Adding a prior on the variance improves robustness against outliers.

There is a nice paper by Andrew Gelman: "Prior distributions for variance parameters in hierarchical models" where he discusses what good choices for the priors on the variances can be.

Answer 4

A robust estimator for the location parameter $\mu$ of some dataset of size $N$ is obtained when one assigns a Jeffreys prior to the variance $\sigma^2$ of the normal distribution, and computes the marginal for $\mu$, yielding a $t$ distribution with $N$ degrees of freedom.

Similarly, if you want a robust estimator for the standard deviation $\sigma$ of some data $D$, we can do the following:

First, we suppose that the data is normally distributed when its mean and standard deviation are known. Therefore, $$\left.D\right|_{\mu,\sigma} \sim \mathcal{N}(\mu,\sigma^2)$$ and if $D \equiv (d_1,\ldots,d_N)$ then $$p(D|\mu,\sigma^2) = \frac{1}{(\sqrt{2\pi}\sigma)^N} \exp\left(-\frac{N}{2\sigma^2}\left((m-\mu^2)+s^2\right)\right)$$ where the sufficient statistics $m$ and $s^2$ are $$m=\frac{1}{N}\sum_{i=1}^N d_i \quad s^2 = \frac{1}{N}\sum_{i=1}^N d_i^2 - m^2$$ In addition, using Bayes' theorem, we have $$p(\mu,\sigma^2|D) \propto p(D|\mu,\sigma^2) p(\mu,\sigma^2)$$ A convenient prior for $(\mu,\sigma^2)$ is the Normal-invese-gamma family, which covers a wide range of shapes and is conjugate to this likelihood. This means that the posterior distribution $p(\mu,\sigma^2|D)$ still belongs to the normal-inverse-gamma family, and its marginal $p(\sigma^2|D)$ is an inverse gamma distribution parameterized as $$\left.\sigma^2\right|_{D} \sim \mathcal{IG}\left(\alpha+N/2,2\beta+Ns^2\right) \qquad \alpha,\beta>0$$ From this distribution, we can take the mode, which will give us an estimator for $\sigma^2$. This estimator will be more or less tolerant to small excursions from misspecifications on the model by varying $\alpha$ and/or $\beta$. The variance of this distribution will then provide some indication on the fault-tolerance of the estimate. Since the tails of the inverse gamma are semi-heavy, you get the kind of behaviour you would expect from the $t$ distribution estimate for $\mu$ that you mention.

Answer 5

I have followed the discussion from the original question. Rasmus when you say robustness I am sure you mean in the data (outliers, not miss-specification of distributions). I will take the distribution of the data to be Laplace distribution instead of a t-distribution, then as in normal regression where we model the mean, here we will model the median (very robust) aka median regression (we all know). Let the model be:

$Y=\beta X+\epsilon$, $\epsilon$ has laplace$(0,$$\sigma^2)$.

Of course our goal is to estimate model parameters. We expect our priors to be vague to have an objective model. The model at hand has a posterior of the form $f(\beta,\sigma,Y,X)$. Giving $\beta$ a normal prior with large variance makes such a prior vague and a chis-squared prior with small degrees of freedom to mimic a jeffrey's prior(vague prior) is given to to $\sigma^2$. With a Gibbs sampler what happens? normal prior+laplace likehood=???? we do know. Also chi-square prior +laplace likelihood=??? we do not know the distribution. Fortunately for us there is a theorem in (Aslan,2010) that transforms a laplace likelihood to a scale mixture of normal distributions which then enable us to enjoy the conjugate properties of our priors. I think the whole process described is fully robust in terms of outliers. In a multivariate setting chi-square becomes a a wishart distribution, and we use multivariate laplace and normal distributions.

Answer 6

Suppose that you have $K$ groups and you want to model the distribution of their sample variances, perhaps in relation to some covariates $\bf{x}$. That is, suppose that your data point for group $k \in {1 \ldots K}$ is $\textrm{Var}(y_k) \in [0, \infty)$. The question here is, "What is a robust model for the likelihood of the sample variance?" One way to approach this is to model the transformed data $\textrm{ln}[\textrm{Var}(y_k)]$ as coming from a $t$ distribution, which as you have already mentioned is a robust version of the normal distribution. If you don't feel like assuming that the transformed variance is approximately normal as $n \rightarrow \infty$, then you could choose a probability distribution with positive real support that is known to have heavy tails compared to another distribution with the same location. For example, there is a recent answer to a question on Cross Validated about whether the lognormal or gamma distribution has heavier tails, and it turns out that the lognormal distribution does (thanks to @Glen_b for that contribution). In addition, you could explore the half-Cauchy family.

Similar reasoning applies if instead you are assigning a prior distribution over a scale parameter for a normal distribution. Tangentially, the lognormal and inverse-gamma distributions are not advisable if you want to form a boundary avoiding prior for the purposes of posterior mode approximation because they peak sharply if you parameterize them so that the mode is near zero. See BDA3 chapter 13 for discussion. So in addition to identifying a robust model in terms of tail thickness, keep in mind that kurtosis may matter to your inference, too.

I hope this helps you as much as your answer to one of my recent questions helped me.