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!

• Try writing down the Normal & Cauchy likelihoods to start with & see how far you can simplify each. – Scortchi - Reinstate Monica Aug 28 '13 at 19:24
• Though "looks much more like a Cauchy" isn't very convincing & you might be better off looking at robust analysis methods. – Scortchi - Reinstate Monica Aug 28 '13 at 19:37
• +1 to @Scortchi. I don't know your data, clearly, but jumping to the conclusion of a Cauchy I will still bet as being very pessimistic. (Harold Jeffreys suggested that real errors of measurement were like t with roughly 7 d.f. rather than Gaussian.) – Nick Cox Aug 28 '13 at 20:41
• @NickCox: I recently read that they Cauchy is like a t with only 1 d.f. – Wayne Aug 28 '13 at 21:33
• Indeed, that's what it (or a standardised version) is, but it's way pessimistic for errors. You use Stata: look at some samples using rt(1) as a function. – Nick Cox Aug 28 '13 at 21:39

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.

• I'm not using likelihood for regression but rather for inter-model comparisons, though, so technically finding the MAXIMUM likelihood is less important than being able to say "Model B is significantly more likely than model A." I suppose in this case I can just compare the models' respective likelihoods instead of finding some proxy statistic like SSE, huh? – Sam R Aug 29 '13 at 17:30
• (1) You still need to fit the models by maximum likelihood, if they're not completely specified, & (2) you want to penalize the likelihood for the model that has more free parameters (see AIC), if one does. – Scortchi - Reinstate Monica Aug 29 '13 at 17:42

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  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 • As can R (the rlm method in the MASS package), & SAS (PROC ROBUSTREG). – Scortchi - Reinstate Monica Aug 28 '13 at 22:00 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. 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. • @DJE The Cauchy distribution is a particular case of the Student-$t$distribution for the case$\nu=1$. I did not say preferred, I said more robust. The robustness comes from assuming that$\nu$is an unknown parameter. This is, the tail behaviour of the errors is dictated by the data. In addition when$\nu\rightarrow \infty$, the Student-$t$converges to the Normal distribution. – Drake Aug 29 '13 at 13:54 • Although the OP does refer to a "Cauchy" distribution, it's reasonable--given a context in which the typical residual size is a priori unknown--to understand that they mean a scale family of Cauchy distributions. Also, experience shows that it will be difficult in most circumstances to estimate$\nu$if you allow that as a parameter. If the OP has good reason to believe the distribution is scaled Cauchy (or$t_1$), then they may be better off not varying$\nu$. It is therefore unclear what the claim of "robustness" really means in this answer. – whuber Aug 29 '13 at 14:22 • @whuber Indeed, the estimation of$\nu$can be difficult in some cases. However, it is a good(better) practice to, at least, select a grid of values of$\nu$(including$\nu=1$) to check if the value$\nu=1$is indeed as good as our believes point out. The data may help to check this. I think the use of the Student-$t$distribution as a means of including models with different tail behaviour has been laargely discussed in the literature. I have included a quote to clarify what "robustness" means. – Drake Aug 29 '13 at 14:36 • That is a well-known definition of robustness but it does not apply directly here: it refers to statistics, not estimators. One could argue that somehow the likelihood function you wrote might have some robust qualities, but whether it does or not is not immediately apparent. Regardless, it is not at all clear that your procedure has the ability to perform well for a wide range of probability distributions. – whuber Aug 29 '13 at 14:56 • In Venables & Ripley's MASS they caution strongly against estimating$\nu$from data when using$t\$-regression; that destroys all robustness properties, they also give references. – kjetil b halvorsen May 7 '15 at 14:19