# Frequentist Predictive Distribution for a Cauchy variable

I have not been able to find this in the literature, but that probably means I am looking in the wrong spot. I am looking to find the Frequentist predictive distribution, assuming it exists, for a one dimensional and an n-dimensional Cauchy variate.

The issue with the n-dimensional version is that there is nothing like a covariate matrix, instead, there is but one scale parameter making the errors hyper-circular. I could see this interfering with the existence of a pivotal value.

EDIT

I am either looking to predict $x_{i+1}$ from a set of observations $x_1\dots{x_i}$ drawn from a Cauchy distribution with center $\mu$ and scale $\sigma,$ or to predict $y_{i+1}$ from some equation $y=mx+b,$ where $x$ is drawn from a Cauchy distribution as above. It could be a vector or multidimensional, but I am trying to determine the relative properties of the Bayesian versus the Frequentist prediction. My data is drawn from either a truncated Cauchy or a Cauchy depending on which set.

A prediction interval will work as I will just set the interval to 100%.

• When you say Cauchy variate, do you mean a regressor $X_{j,i}$ in the regression equation $Y_i = \sum_{j=1}^JX_{j,i}\beta_j + \varepsilon_i$, and are you looking for the confidence interval of the OLS estimate for $\beta_j$? I am not sure I understand the question. Dec 18 '17 at 13:42
• There are papers, like this one from 2008, on estimation with multivariate t distributions. The multivariate Cauchy is a special case of the multivariate t. As an aside, this allows for a fully flexible correlation structure. Does that help you, or does it not match your question? Dec 18 '17 at 16:44
• Are you trying to do this stats.stackexchange.com/questions/16349 for a multivariate distribution with zero means? Dec 18 '17 at 16:49
• @eric_kernfeld I have to read it carefully, but yes that is sort of it, except that I want to know how to find its predictive density using Frequentist methods. Dec 18 '17 at 18:50
• it seems that you're trying to estimate Cauchy distribution parameters from $x_i$. Is that right? Dec 21 '17 at 17:13

The general solution to your problem is Maximum Likelihood Estimation (MLE) of your parameters $\theta$. Once they are obtained as $\hat{\theta}$, you substitute them into your pdf for the unknown parameters, i.e. you estimate the pdf of your random variable as $\hat{f}(x_i) = f(x_i|\hat{\theta})$. This allows you to construct the the predictive distribution of your Cauchy Random Variable.

For the univariate case, this paper is an excellent resource. For the univariate Cauchy with center $\mu$ and scale $\sigma$, one has a closed form if you have $3-4$ observations. If you have $n>4$ observations, the MLE exists$^{\ast}$. If you have $n$ observations, you will have to solve two equations that are easily derived by setting the first derivative of the log-likelihood to zero, see here for their exact form. (In their notation, $x_0 = \mu$ and $\sigma = \gamma$.) Solving this problem numerically has an implementation in the R language, see here.

For the multivariate case, all you need to note is that the multivariate Cauchy distribution is simply a multivariate $t$-distribution where the degree of freedom parameter is set to $1$, as was already pointed out in the comments. For the multivarate-$t$, you can do MLE inference as explained excellently in this answer, which is based on the paper that eric_kernfeld has pointed out. I did not find ready-to-roll implementation for this algorithm, but as you will see when you take a look at the supplied answer in the post, it should really easy to implement it yourself.

Difference to Bayesian prediction: In the Bayesian setting, you would put a prior on the parameters $\mu$ and $\sigma$, modelling your uncertainty about them as a random variable. Thus, you will get posterior distributions for both parameters, which indicate the relative certainty you have about them given your data. If you have the posterior $q(\mu, \sigma|x_1,\dots,x_n)$, you then obtain your predictive distribution as $\int f(x|\mu, \sigma)q(\mu, \sigma|x_1,\dots,x_n)d\mu d\sigma$, integrating out your uncertainty. In contrast, the MLE-setting will give you point estimates of $\mu$ and $\sigma$ that you plug into your pdf's functional form. Equivalently, you could say that MLE leads to a posterior with point mass $1$ at the tuple $(\hat{\mu}, \hat{\sigma})$ and $0$ probability at any other value. Thus, you ignore all parameter uncertainty in this case, and you rely on the fact that $\hat{\theta}$ is asymptotically equivalent to $\theta$, meaning that $\hat{f}(x) \to f(x)$ (uniformly over $x$).

$^\ast$Well, that is unless for the exotic case where $n$ is even and $n/2$ of your observations take value $x_1$ while the other half takes value $x_2$, which happens with probability zero because the Cauchy distribution is continuous.

• Jeremias. Do you think there are possibilities to incorporate the uncertainty about $\hat{\theta}$ into the prediction interval? And how do we construct a predictive distribution from the pdf of a multidimensional cauchy distribution? Dec 21 '17 at 12:31
• If you assume that the $x_i$ are random draws from a Cauchy random variable with unknown parameters, plugging in the estimated parameters into the functional form directly gives you the predictive distribution of the next draws of $x_i$. Dec 21 '17 at 13:33
• As for incorporating parameter uncertainty, if you want to do that you have to go the Bayesian way. Note that as a side product, the parameter posterior of Bayesian inference concentrates on the MLE via the Bernstein Mises theorem asymptotically. Dec 21 '17 at 13:36

One could use a Monte Carlo method to obtain empirical estimates for relationships between the $x_1....x_i$ and the prediction interval for $x_{i+n}$.

Motivation: If we estimate the prediction interval based on the quartiles/CDF of a distribution that follows from maximum likelihood estimates (or other type of parameter estimates), then we underestimate the size of the interval. Effectively, in practice, the point $x_{i+n}$ will fall out of the range more often than predicted.

The figure below demonstrates by how much we underestimate the size of the interval, by expressing how many more times a new measurement $x_i$ is outside the predictive range based on parameter estimates. (based on computations with 2000 repetitions for the prediction)

For instance, if we use a prediction interval of 99% (thus expecting 1% errors), then we get 5 times more errors if the sample size was 3.

These type of computations can be used to make empirical relationships for how we can correct the range, as well the computations show that for large $n$ the difference becomes smaller(and at some point one may consider it irrelevant).

set.seed(1)

# likelihood calculation
like<-function(par, x){
scale = abs(par[2])
pos   = par[1]
n <- length(x)
like <- -n*log(scale*pi) - sum(log(1+((x-pos)/scale)^2))
-like
}

# obtain effective predictive failure rate rate
tryf <- function(pos, scale, perc, n) {

# random distribution
draw <- rcauchy(n, pos, scale)

# estimating distribution parameters based on median and interquartile range
first_est <- c(median(draw), 0.5*IQR(draw))

# estimating distribution parameters based on likelihood
out <- optim(par=first_est, like, method='CG', x=draw)
# making scale parameter positive (we used an absolute valuer in the optim function)
out$par[2] <- abs(out$par[2])

# calculate predictive interval
ql <- qcauchy(perc/2, out$par[1], out$par[2])
qh <- qcauchy(1-perc/2, out$par[1], out$par[2])

# calculate effective percentage outside predicted predictive interval
pl <- pcauchy(ql, pos, scale)
ph <- pcauchy(qh, pos, scale)
error <- pl+1-ph
error
}

# obtain mean of predictive interval in 2000 runs
meanf <- function(pos,scale,perc,n) {
trueval <- sapply(1:2000,FUN <- function(x) tryf(pos,scale,perc,n))
mean(trueval)
}

#################### generate image

# x-axis chosen desired interval percentage
percentages <- 0.2/1.2^c(0:30)

# desired sample sizes n
ns <- c(3,4,5,6,7,8,9,10,20,30)

# computations
y <- matrix(rep(percentages, length(ns)), length(percentages))
for (i in which(ns>0)) {
y[,i] <- sapply(percentages, FUN <- function(x) meanf(0,1,x,ns[i]))
}

# plotting
plot(NULL,
xlim=c(0.0008,1), ylim=c(0,10),
log="x",
xlab="aimed error rate",
ylab="effective error rate / aimed error rate",
yaxt="n",xaxt="n",axes=FALSE)
axis(1,las=2,tck=-0.0,cex.axis=1,labels=rep("",2),at=c(0.0008,1),pos=0.0008)
axis(1,las=2,tck=-0.005,cex.axis=1,at=c(0.001*c(1:9),0.01*c(1:9),0.1*c(1:9)),labels=rep("",27),mgp=c(1.5,1,0),pos=0.0008)
axis(1,las=2,tck=-0.01,cex.axis=1,labels=c(0.001,0.01,0.1,1), at=c(0.001,0.01,0.1,1),mgp=c(1.5,1,0),pos=0.000)
#axis(2,las=1,tck=-0.0,cex.axis=1,labels=rep("",2),at=c(0.0008,1),pos=0.0008)
#axis(2,las=1,tck=-0.005,cex.axis=1,at=c(0.001*c(1:9),0.01*c(1:9),0.1*c(1:9)),labels=rep("",27),mgp=c(1.5,1,0),pos=0.0008)
#axis(2,las=1,tck=-0.01,cex.axis=1,labels=c(0.001,0.01,0.1,1), at=c(0.001,0.01,0.1,1),mgp=c(1.5,1,0),pos=0.0008)
axis(2,las=2,tck=-0.01,cex.axis=1,labels=0:15, at=0:15,mgp=c(1.5,1,0),pos=0.0008)

colours <- hsv(c(1:10)/20,1,1-c(1:10)/15)
for (i in which(ns>0)) {
points(percentages,y[,i]/percentages,pch=21,cex=0.5,col=colours[i],bg=colours[i])
}

legend(x=0.4,y=4.5,pch=21,legend=ns,col=colours,pt.bg=colours,title="sample size")

title("difference between confidence interval and effective confidence interval")

plot(ns,y[31,]/percentages[31],log="")

• What does the plot tell us other than using a small sample size will yield a bad estimate of your parameters when using mle? I fail to see how it invalidates using mle, since the error rates look excellent even for a very small sample size of 30. I'm also not sure I understand what the alternative is that you propose, would you mind expanding on the computational methods you mention at the beginning of your answer? Dec 21 '17 at 18:08
• @JeremiasK In practical applications, with small sample sizes, one could use these computations to as empirically determined correction factors. Dec 21 '17 at 18:13
• That makes sense! I don't think you mention it in the post, maybe you should edit it in so people don't have to read through the comments Dec 21 '17 at 18:56
• @MartijnWeterings so far you make the most sense. The pivot $\frac{\sqrt{n}(\hat{\mu}-\mu)}{\hat{\sigma}}$$follows the standard normal once the sample size gets to around 100, but I realized that I am beyond my skills to unwind this because instead of choosing a variable I am choosing a function for the minimization and I have not done that before. Dec 22 '17 at 0:43 • @DaveHarris I believe my method is in not so different from the case by Jeremia, except that I make an expression (and only by an experimental mathematics approach) for the underestimated range that occurs because the distribution$f(x,\hat{x}_0,\hat{\gamma})$is an over-dispersed version of$f(x,x_0,\gamma)$. Dec 22 '17 at 2:56 It seems that all you need is to estimate the parameters of Cauchy distribution from the dataset$x_i$. Here's what Stephens proposes, it's not MLE, and author claims this method is consistent and more stable than MLE though you have to take into account that this has been written in the last century. where Cauchy is parameterized as follows: Once you have the distribution, your point forecast will be$\hat\alpha\$. Note, that since it doesn't have moments, you won't be able to show that your forecast is optimal in usual sense such as minimizing expected square cost.