Model fitting when errors take a Cauchy distribution

Question

It's my understanding that the sum of squared errors (SSE) serves as a maximum likelihood estimator when a model's errors are normally distributed. (That is, if you find model parameters that minimize the SSE, they also maximize the likelihood.) However, the error distribution of my model looks much more like a Cauchy distribution. Would minimizing the SSE still result in the maximum likelihood parameter set for my model? If not, what statistic should I look at?

Forgive me if this doesn't make any sense, or I'm missing something simple. Please feel free to link to sources that might help me understand the basics a bit better. Thanks!

Answer 1

The least squares estimates for the regression coefficients are only equal to the maximum-likelihood estimates when the errors have a normal distribution (see here for the proof).

If you really wanted maximum likelihood estimates for regression parameters with Cauchy errors, just look at that likelihood: $$L(\beta,\sigma)=\prod_{i=1}^n {\frac{1}{\pi\sigma\left(1+\left(\frac{y_i-\beta^\mathrm{T}x_i}{\sigma}\right)^2\right)}}$$

($y_i$ is the $i$th observation, $x_i$ the vector of predictors, $\sigma$ the scale parameter, & $\beta$ the vector of coefficients.) There's no sufficient statistic of lower dimensionality than the entire dataset, so it's not so easy to maximize, though there's probably a better method than brute force. But without some theoretical motivation for assuming Cauchy errors, you can just say they have some fat-tailed distribution. In this situation some form or other of robust regression would be worth considering.

Note that the least squares approach isn't the worst thing you could use even so. Provided the variance is constant (& finite, which it isn't for the Cauchy) it still gives consistent estimates, even the best linear unbiased estimates, though you'd have to take confidence intervals with a pinch of salt.

Answer 2

The Jeffreys (posterior) distribution is quite nice to do inference in a linear regression model with location-scale errors. Inference based on the Jeffreys distribution achieves very good frequentist properties: it provides confidence intervals whose coverage is close to the nominal coverage, even for very small sample sizes.

Let $\Phi$ be any derivable cdf and $\phi=\Phi'$ the corresponding pdf. Consider the linear regression model $Y=X\beta+\sigma\epsilon$ with $\epsilon_i \sim_{\text{iid}} \mathrm{d}\Phi$.

The Jeffreys posterior distribution is given, up to a proportionality constant, by: $$ \pi(\beta, \sigma \mid y) \propto \frac{1}{\sigma^{n+1}} \prod_{i=1}^n \phi\left(\frac{y_i-\mu_i}{\sigma}\right) =: f(\beta, \sigma \mid y) $$ where $\mu_i=X_i \beta$ is the expected value of $y_i$.

The problem, given a set $A$ in the $q$-dimensional parameter space $\Theta$ (with $q=p+1$ since parameters are $\beta_1$, $\ldots$, $\beta_p$ and $\sigma$), is to evaluate $$ \int_A \pi(\beta, \sigma \mid y) \mathrm{d}\beta \mathrm{d}\sigma = \frac{\int_A f(\beta, \sigma \mid y)\mathrm{d}\beta \mathrm{d}\sigma}{\int_\Theta f(\beta, \sigma \mid y)\mathrm{d}\beta\mathrm{d}\sigma} \approx ? $$

Change of variables

It is possible to transform this integral to an integral on ${[0,1]}^q$ as follows. The key point is the fact that $$ \frac{1}{\sigma^{q+1}}\prod_{i=1}^q\phi\left(\frac{y_i-\mu_i}{\sigma}\right) $$ is, up to a proportionality constant not depending on $(\beta,\sigma)$, the Jacobian of the function $$ F\colon(\beta,\sigma)\mapsto \left(\Phi\left(\frac{y_1-\mu_1}{\sigma}\right), \ldots, \Phi\left(\frac{y_q-\mu_q}{\sigma}\right)\right) \in {[0,1]}^q. $$ I have not tried to prove this point, but one can numerically check it, for a simple linear regression for example:

x <- 1:10; y <- rcauchy(length(x), 2+x)
X <- model.matrix(~x)
Phi <- pcauchy; phi <- dcauchy # use any cdf and pdf you want
F <- function(betasigma){
  sapply(1:3, function(i) Phi((y[i]-X[i,]%*%betasigma[1:2])/betasigma[3]))
}
library(pracma) # provides the jacobian() function
f <- function(betasigma){
  prod(sapply(1:3, function(i){ phi((y[i]-X[i,]%*%betasigma[1:2])/betasigma[3]) })) / betasigma[3]^4
}
# look, the ratio is always the same:
det(jacobian(F, c(1,1,1)))/f(c(1,1,1))
## [1] -19.01263
det(jacobian(F, c(1,2,1)))/f(c(1,2,1))
## [1] -19.01263
det(jacobian(F, c(2,2,2)))/f(c(2,2,2))
## [1] -19.01263

Thus, $$ \begin{align} \int_A f(\beta, \sigma \mid y)d\beta d\sigma & \propto \int_A \bigl|\det J_F(\mu,\sigma)\bigr| \frac{1}{\sigma^{n-q}} \prod_{i=q+1}^n \phi\left(\frac{y_i-\mu_i}{\sigma}\right) \mathrm{d}\beta\mathrm{d}\sigma \\ & = \int_{F(A)} g\bigl(F^{-1}(u_1, \ldots, u_q)\bigr)\mathrm{d}u_1\ldots\mathrm{d}u_q \end{align} $$ where $g(\beta,\sigma)=\frac{1}{\sigma^{q+1}} \prod_{i=1}^q \phi\left(\frac{y_i-\mu_i}{\sigma}\right)$.

It is not difficult to get the inverse of $F$: $$ F^{-1}(u_1, \ldots, u_q) = {(\beta,\sigma)}' = {(H'H)}^{-1}H'y_{1:q} $$ where the matrix $H$ is $H=\left[ X_{1:q}, {\bigl(\Phi^{-1}(u_i)\bigr)}_{i\in(1:q)}\right]$.

Note that $F^{-1}(u_1, \ldots, u_q)$ yields $\sigma<0$ for some values of the $u_i$. In fact, if $F^{-1}\bigl(\Phi(z_1), \ldots, \Phi(z_q))={(\beta,\sigma)}'$, then $F^{-1}\bigl(\Phi(-z_1), \ldots, \Phi(-z_q))={(\beta,-\sigma)}'$, therefore the set of $u_i$'s for which $\sigma>0$ has Lesbegue measure $1/2$.

In fact, the Jeffreys distribution for a location-scale linear regression is the same as the fiducial distribution. The method I present is a particular case of the general method given in the paper Computational issues of generalized fiducial inference by Hannig & al.. But there is a high simplification of the general method in the case of a location-scale linear regression (we can take $K=1$ with the notations of the paper, but I will not develop this point).

Algorithm

The Jeffreys function below returns an approximation of the Jeffreys distribution for the linear regression model when errors follow a Student distribution with degrees of freedom df, to be set by the user. For df=Inf (default), this is the Gaussian linear regression; for df=1 this is the Cauchy linear regression. In the Gaussian case df=Inf, we can compare the results to the exact Jeffreys distribution which is known and elementary. Moreover the inference based on the Jeffreys distribution in the Gaussian case is the same as the usual least-squares inference (as we will see on examples).

By default, the X matrix is the matrix of the intercept-only model y~1.

The approximation is obtained by a Riemann-like integration on ${[0,1]}^q$ using a uniform partition into hypercubes. The partition is controlled by the argument L, giving the number of centers of the hypercubes on each coordinate (hence there are $L^q$ hypercubes).

#' parameters: y (sample), X (model matrix), L (number of points per coordinate)
Jeffreys <- function(y, X=as.matrix(rep(1,length(y))), L=10, df=Inf){
  qdistr <- function(x, ...) qt(x, df=df, ...)
  ddistr <- function(x, ...) dt(x, df=df, ...)
  n <- nrow(X)
  q <- ncol(X)+1
  # centers of hypercubes (volume 1/L^p)
  centers <- as.matrix(do.call(expand.grid, rep(list(seq(0, 1, length.out=L+1)[-1] - 1/(2*L)), q)))
  # remove centers having equal coordinates (H'H is not invertible)
  centers <- centers[apply(centers, 1, function(row) length(unique(row))>1),]
  # outputs
  M <- (L^q-L)/2 # number of centers yielding sigma>0
  J <-  numeric(M)
  Theta <- array(0, c(M, q))
  # algorithm
  I <- 1:q
  yI <- y[I]; ymI <- y[-I]
  XI <- X[I,]; XmI <- X[-I,]
  counter <- 0
  for(m in 1:nrow(centers)){
    H <- unname(cbind(XI, qdistr(centers[m,])))
    theta <- solve(crossprod(H))%*%t(H)%*%yI
    if(theta[q]>0){ # sigma>0
      counter <- counter+1
      J[counter] <- sum(ddistr((ymI-XmI%*%head(theta,-1))/theta[q], log=TRUE)) - (n-q)*log(theta[q])  
      Theta[counter,] <- theta
    }
  }
  J <- exp(J)
  return(list(Beta=Theta[,-q], sigma=Theta[,q], W=J/sum(J)))
}

The function returns the values of $(\beta,\sigma)$ corresponding to every hypercube center in the partition of ${[0,1]}^q$. It also computes the values of the integrand evaluated at every center in the vector J, and returns the normalized vector of weights W=J/sum(J). We will see how to deal with these outputs on some examples.

First example: Gaussian sample

Let's try it for an i.i.d. Gaussian sample $y_i \sim_{\text{i.i.d.}} {\cal N}(\mu, \sigma^2)$:

set.seed(666)
n <- 4
y <- rnorm(n)
results <- Jeffreys(y, L=100)
Mu <- results$Beta; Sigma <- results$sigma; W <- results$W

Now we can treat Mu and Sigma as if they were weighted samples of the Jeffreys distribution, with weights W.

The theoretical mean is the sample mean, and our approximation is quite good:

sum(W*Mu); mean(y)
## [1] 1.109794
## [1] 1.110175

We can get the approximate Jeffreys cdf with the ewcdf function (weighted empirical cdf) of the spatstat package, and compare with the theoretical one. Our approximation is quite perfect:

### approximate Jeffreys distribution of µ ###
F_mu <- spatstat::ewcdf(Mu, weights=W)
curve(F_mu, from=0, to=2.5, xlab="mu", ylim=c(0,1), col="blue", lwd=2)
### exact Jeffreys distribution ###
mean_y <- mean(y); sd_y <- sd(y)
curve(pt((x-mean_y)/(sd_y/sqrt(n)), df=n-1), add=TRUE, col="red", lwd=4, lty="dashed")

We can get confidence intervals by applying the quantile function to the weighted cdf F_mu. They are theoretically the same as the ususal confidence intervals in Gaussian linear regression, and we indeed get very close results:

quantile(F_mu, c(2.5,97.5)/100)
##       2.5%      97.5% 
## -0.7143989  2.9309891
confint(lm(y~1)) # theoretically the same
##                  2.5 %  97.5 %
## (Intercept) -0.7121603 2.93251

The same for $\sigma$ (knowing the inverse-Gamma distribution of $\sigma^2$):

F_sigma <- spatstat::ewcdf(Sigma,W)
curve(F_sigma, from=0, to=2.5, xlab="sigma", ylim=c(0,1), col="blue", lwd=2)
curve(1-pgamma(1/x^2, (n - 1)/2, (n - 1) * sd_y^2/2), add=TRUE, col="red", lwd=4, lty="dashed")

Second example: Cauchy sample

Now let's try a i.i.d. Cauchy sample with sample size $n=200$.

set.seed(666)
n <- 200
y <- rcauchy(n)
results <- Jeffreys(y, L=100, df=1)
Mu <- results$Beta; Sigma <- results$sigma; W <- results$W

Since $n=200$ is not a small sample size, the Jeffreys means are close to the maximum-likelihood estimates:

sum(W*Mu); sum(W*Sigma)
## [1] -0.01490355
## [1] 0.9081371
MASS::fitdistr(y, "cauchy")
##     location        scale   
##   -0.01345121    0.89958785 
##  ( 0.09185580) ( 0.08874509)

The MASS::rlm estimates are not so close:

rlmfit <- MASS::rlm(y~1)
rlmfit$coefficients
## (Intercept) 
##  -0.1160915
rlmfit$s # rlm estimate of sigma
## [1] 1.338744

Jeffreys confidence intervals are close to the ML asymptotic confidence intervals:

F_mu <- spatstat::ewcdf(Mu, weights=W); F_sigma <- spatstat::ewcdf(Sigma,W)
quantile(F_mu, c(2.5,97.5)/100)
##       2.5%      97.5% 
## -0.1971883  0.1707172
quantile(F_sigma, c(2.5,50,97.5)/100)
##      2.5%       50%     97.5% 
## 0.7471966 0.9055118 1.1027395
confint(MASS::fitdistr(y, "cauchy"))
##               2.5 %    97.5 %
## location -0.1934853 0.1665829
## scale     0.7256507 1.0735250

Third example : Gaussian simple linear regression

Nice:

f <- function(x) 4+2*x
set.seed(666)
n <- 20
x <- seq_len(n) # covariates
y <- f(x)+rnorm(n)
# run algorithm
results <- Jeffreys(y, X=model.matrix(~x), L=60)
# outputs
W <- results$W; Beta0 <- results$Beta[,1]; Beta1 <- results$Beta[,2]
sum(W*Beta0); sum(W*Beta1) # Jeffreys means
## [1] 4.172721
## [1] 1.983503
coef(lm(y~x)) # theoretically the same
## (Intercept)           x 
##    4.179859    1.983008
F_Beta0 <- spatstat::ewcdf(Beta0, weights=W); F_Beta1 <- spatstat::ewcdf(Beta1, weights=W)
quantile(F_Beta0, c(2.5,97.5)/100); quantile(F_Beta1, c(2.5,97.5)/100)
##     2.5%    97.5% 
## 2.883869 5.499620
##     2.5%    97.5% 
## 1.872328 2.095764
confint(lm(y~x)) # theoretically the same
##                2.5 %   97.5 %
## (Intercept) 2.857903 5.501815
## x           1.872653 2.093362

Fourth example : Cauchy simple linear regression

set.seed(666)
y <- f(x)+rcauchy(n)
# run algorithm
results <- Jeffreys(y, X=model.matrix(~x), L=60, df=1)
# outputs
W <- results$W; Beta0 <- results$Beta[,1]; Beta1 <- results$Beta[,2]; Sigma <- results$sigma
# Jeffreys means
sum(W*Beta0); sum(W*Beta1); sum(W*Sigma)
## [1] 4.157664
## [1] 1.997121
## [1] 0.685825

While $n=20$ is not large, the ML estimates of the regression parameters are close to their Jeffreys means, but they are not so close for $\sigma$:

X <- model.matrix(~x)
likelihood <- function(y, beta0, beta1, sigma){
  prod(dcauchy((y-X%*%c(beta0,beta1))/sigma)/sigma)
}
(ML <- MASS::fitdistr(y, likelihood, list(beta0=sum(W*Beta0), beta1=sum(W*Beta1), sigma=1)))
##      beta0        beta1        sigma   
##   4.20188590   1.99239112   0.60087397 
##  (0.54295228) (0.04433536) (0.18660186)

The Jeffreys confidence intervals are close the ML confidence intervals, except for $\sigma$:

confint(ML)
##          2.5 %    97.5 %
## beta0 3.137719 5.2660528
## beta1 1.905495 2.0792868
## sigma 0.235141 0.9666069
F_Beta0 <- spatstat::ewcdf(Beta0, weights=W); F_Beta1 <- spatstat::ewcdf(Beta1, weights=W)
quantile(F_Beta0, c(2.5,97.5)/100); quantile(F_Beta1, c(2.5,97.5)/100)
##     2.5%    97.5% 
## 3.098442 5.167351
##     2.5%    97.5% 
## 1.913328 2.089146
F_sigma <- spatstat::ewcdf(Sigma,W); quantile(F_sigma, c(2.5,50,97.5)/100)
##      2.5%       50%     97.5% 
## 0.3418978 0.6491162 1.2464163

The MASS::rlm estimates of the regression parameters are rather close to their Jeffreys means too:

rlmfit <- MASS::rlm(y~x)
rlmfit$coefficients
## (Intercept)           x 
##    3.945603    2.042590
rlmfit$s # rlm estimate of sigma
## [1] 1.019974

Answer 3

GraphPad Prism can do nonlinear regression assuming a Cauchy distribution. That is our robust method. The mathematical details are explained in detail on pages 11-14 of BMC Bioinformatics 2006, 7:123 doi:10.1186/1471-2105-7-123, Detecting outliers when fitting data with nonlinear regression – a new method based on robust nonlinear regression and the false discovery rate

Answer 4

A bit too late, but it may be useful for others in the future.

Scortchi wrote the likelihood as if the errors where independent and there is no simple answer in that case. To complement what he said in the end about OLS, note that if you instead, specify the joint likelihood with dependent errors, you will end up with a Multivariate Cauchy Distribution (a Multivariate t Distribution with $\nu = 1$) and from here you can recover the the maximum likelihood estimates of the regression coefficients which will be exactly what you would get with OLS.

$$\mathcal{L} = \frac{\Gamma((N+1)/2)}{\pi^{(N+1)/2}\sigma^{N/2}}\bigg(1+\frac{\mathbf{\epsilon}^\intercal\mathbf{\epsilon}}{\sigma^{2}}\bigg)^{(N+1)/2}$$

$$\mathcal{L} = \frac{\Gamma((N+1)/2)}{\pi^{(N+1)/2}\sigma^{N/2}}\bigg[1+\frac{\sum_{i=1}^{N}(y_{i}-\beta x_{i})^{2}}{\sigma^{2}}\bigg]^{(N+1)/2}$$

If $\ell = \log\mathcal{L}$, then

$$\ell = \log{\bigg[\frac{\Gamma((N+1)/2)}{\pi^{(N+1)/2}}}\bigg] -\frac{N}{2}\log{|\sigma|} - \frac{N+1}{2}\log\bigg[1+\frac{\sum_{i=1}^{N}(y_{i}-\beta x_{i})^{2}}{\sigma^{2}}\bigg]$$

$$\frac{\partial\ell}{\partial\beta} = -\frac{N+1}{2}\frac{1}{1+\frac{\sum_{i=1}^{N}(y_{i}-\beta x_{i})^{2}}{\sigma^{2}}}\sum_{i=1}^{N}(y_{i}-\beta x_{i})(-x_{i})\frac{1}{\sigma^{2}} = 0$$

$$\frac{\sum_{i=1}^{N}(y_{i}-\beta x_{i})(-x_{i})}{\sigma^{2}+\sum_{i=1}^{N}(y_{i}-\beta x_{i})^{2}} = 0$$

$$\sum_{i=1}^{N}(y_{i}-\beta x_{i})(x_{i}) = 0$$

$$\sum_{i=1}^{N}y_{i}x_{i}-\sum_{i=1}^{N}\beta x_{i}^{2} = 0$$

$$\beta \sum_{i=1}^{N}x_{i}^{2} = \sum_{i=1}^{N}y_{i}x_{i}$$

$$\beta = \frac{\sum_{i=1}^{N}y_{i}x_{i}}{\sum_{i=1}^{N}x_{i}^{2}}$$

And for $\sigma$

$$\sigma = \sqrt{\mathbf{\epsilon}^\intercal\mathbf{\epsilon}\bigg(\frac{N+2}{N}\bigg)}$$

The estimator for $\beta$ is the OLS estimate and the estimator for $\sigma$ is really similar to the OLS $s^2$.

Unfortunately these estimators don't have any moments due to the Cauchy density. A bit more complicated algebra when using a Multivariate t density instead of a Cauchy density, would yield the same estimator for $\beta$ and a similar for $\sigma$ and also the third parameter $\nu$, but with finite moments if ($\nu > 2$) and 'inflated' variance compared to standard OLS. You will note that the same inference holds true as in the case of normality: you can use t, and F statistics.

However, as you can see, there is not much of gain compared to standard OLS and maximum likelihood with a normal density. OLS does not require normality or independence.

Answer 5

The use of Cauchy errors IS NOT a robust method. It leads to a model that can capture outliers, but if there are no outliers, then the resulting model becomes very restrictive since it is being assumed that the distribution of the errors is heavy tailed with a specific tail behaviour. Because the Cauchy distribution is a special case of the t-distribution for $\nu=1$, this makes a very strong statement about how the errors are distributed. A more robust approach consists of using a Student $t$-distribution $t(0,\sigma,\nu)$, where $\nu$ are the degrees of freedom and $\sigma$ is the scale parameter, which are unknown and they are to be estimated using the data.

The model is

$$y_j = h(x_j^{\top}\beta) + e_j,$$

where $e_j\stackrel{ind}{\sim} t(0,\sigma,\nu)$, $j=1,\dots,n$, and $h$ is a real function.

The likelihood of $(\beta,\sigma,\nu)$ is then given by

$${\mathcal L}(\beta,\sigma,\nu) \propto \prod_{j=1}^n f(h(y_j- x_j^{\top}\beta);0,\sigma,\nu),$$

where

$$f(z;\mu,\sigma,\nu) = \dfrac{1}{\sigma}\dfrac{\Gamma\left(\dfrac{\nu+1}{2}\right)}{\sqrt{\pi\nu}\Gamma\left(\dfrac{\nu}{2}\right)} \left[1+\dfrac{1}{\nu}\left(\dfrac{z-\mu}{\sigma}\right)^2\right]^{-\frac{\nu+1}{2}}.$$

The maximum likelihood estimators can be obtained using numerical methods. Note that this structure covers data fitting, linear and nonlinear regression.

From Wikipedia:

Robust statistics are statistics with good performance for data drawn from a wide range of probability distributions, especially for distributions that are not normally distributed. Robust statistical methods have been developed for many common problems, such as estimating location, scale and regression parameters. One motivation is to produce statistical methods that are not unduly affected by outliers. Another motivation is to provide methods with good performance when there are small departures from parametric distributions. For example, robust methods work well for mixtures of two normal distributions with different standard-deviations, for example, one and three; under this model, non-robust methods like a t-test work badly.