Frequentist Predictive Distribution for a Cauchy variable

Question

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%.

Answer 1

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.

Answer 2

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="")

Answer 3

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.