Model fitting when errors take a Cauchy distribution 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!
 A: 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.
A: 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
A: 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.
A: 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.

