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
rt(1)
as a function. $\endgroup$ – Nick Cox Aug 28 '13 at 21:39