Prediction interval: but instead, the probability that the next datum is above a fixed threshold?

Question

I've been struggling with this problem, and I think I must be missing some important conceptual step. Imagine we observe $\theta_1 \sim N(\mu, \sigma^2)$, with unknown $\mu$ and known $\sigma^2$ (for simplicity). We will draw a second datum, $\theta_2 \sim N(\mu, \sigma^2)$ (for now, unobserved). (Point of clarification—$\theta$ is just a random variable here, not a stand-in for the parameter $\mu$, it might have been more clear if I wrote $X_1,X_2 \sim N(\mu, \sigma^2)$.)

It is straightforward to compute a "prediction interval" for $\theta_2$. By definition, this is a frequentist procedure that guarantees that $\theta_2$ is captured by the interval the correct proportion of the time. Note that $\theta_2 - \theta_1 \sim N(0, 2\sigma^2)$. Thus, $P(\theta_2 > \theta_1 - z_\alpha(\sqrt{2} \sigma) )= \alpha $ (where $z_\alpha$ is the z-score, e.g. $z_{.10} = 1.28$ gives a 90% one-sided prediction interval). (See Wikipedia for details).

However, I am interested in a slightly different property. For some fixed threshold $T$, I want to estimate the probability that $\theta_2 > T$. That is, rather than a procedure that creates an interval with "90% coverage", I want a procedure that estimates the probability that the next datum will fall within a fixed interval, $\theta_2 \in [T, \infty]$. Presumably, the "guarantee" I want is that the average probability of this estimate matches the true probability, $P(\theta_2 > T)$.

I'm stumped! Given that $\theta_2 - \theta_1 \sim N(0, 2\sigma^2)$, my instinct is just to compute $P(Z > (T - \theta_1)/(\sqrt{2}\sigma))$, where $Z \sim N(0,1)$ (i.e. just use the normal CDF). But under simulation, that isn't right. I assume the issue is that $\theta_1$ and $\theta_2 - \theta_1$ aren't independent? But it feels like this problem is doable, and I can't figure it out.

Answer 1

It is not very difficult to find an unbiased estimator for $P(\theta > T)$ : you want to find some function $g(\theta ; T)$ that has the property that $\mathbb E[ g(\theta;T)] = P(\theta > T)$. But $P(\theta > T) = \mathbb E[ 1(\theta > T)]$, so clearly $g(\theta;T) = 1(\theta > T)$ does the job.

Furthermore this is the only unbiased estimator that can be constructed with one sample: you can write the above unbiasedness condition as the equality

$$ \int d\theta \left[ g(\theta;T) - 1(\theta > T) \right ]f(\theta|\mu,\sigma^2) = 0$$

which need to hold for all values of $\mu$ (where $f(\theta|\mu,\sigma^2)$ is the pdf of the normal distribution). Clearly this can happen only if the integrand itself is zero.

Of course this is not a very good estimator : the estimate is one if $\theta > T$ and zero otherwise. But this just demonstrate the general potential difficulties when dealing with frequentist quantities - unbiasedness doesn't guarantee anything other than unbiasedness.

Note that if you want instead a confidence interval for $P(\theta > T)$, you can immediately get it from the standard confidence interval for $\mu$, since $P(\theta > T) = \Phi(\frac{T-\mu}{\sigma})$ is a monotonically decreasing function of $\mu$. So just mapping the interval using the normal CDF will maintain the coverage property.

Answer 2

A maximum likelihood estimate of the probability $p_T = P(Y\geq T)$ can be derived using the estimate of the mean $\hat{\mu}$

$$\hat p_T = 1-\Phi\left( \frac{T-\hat\mu}{\sigma}\right) $$

(You have a factor $\sqrt{2}$ which I find not intuitive, with this formula you just reparameterize the distribution and describe it in terms of $\sigma$ and $p_T$ instead of $\sigma$ and $\mu$)

The expectation value of this is

$$\begin{array}{}E\left[\hat p_T)\right] &=& E\left[1-\Phi\left( \frac{T-X_1}{\sigma}\right)\right]\\ &=&1- \int_{-\infty}^{\infty} \Phi\left(\frac{T-\mu}{\sigma} - x\right)\phi(x) dx \\&= &1 - \Phi\left(\frac{T-\mu}{\sqrt{2}\sigma}\right) \end{array}$$

(For that final equality see https://en.m.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions )

and it will thus be biased.

This bias depends on the value of $\mu$ and can not be corrected.