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

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.

• Because $\Pr(\theta_2 \gt T)$ is a distributional property, you are asking how to estimate this property. This is not a prediction interval; moreover, an appropriate interval related to such an estimate would be a confidence interval. But since you characterize this as a prediction interval problem, I am concerned that you might not have communicated the situation as you intended to. Otherwise, aren't you just re-asking the question answered at stats.stackexchange.com/questions/511265 and stats.stackexchange.com/questions/492260?
– whuber
Apr 12, 2022 at 22:14
• But under simulation, that isn't right. How did you simulate this? Apr 12, 2022 at 22:16
• My apologies for the confusion—$\theta_1, \theta_2$ are just data, poor notation from me simplifying a different problem. I added a note in the question, it would be more typical if they were replaced by $X_1, X_2$. Apr 17, 2022 at 17:13

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.

• +1. I found the use of "$\theta$," although correct, to be a little confusing. In the question the random variable is "$\theta_1$" and it's just too conventional to read "$\theta$" as being a generic reference to the parameter, which here is $\mu$ ($\sigma$ is known).
– whuber
Apr 13, 2022 at 16:44
• Could you clarify what $\theta$ refers to in your post? Also, how did you arrive at the formula in the integral? Apr 13, 2022 at 20:55
• @Learningstatsbyexample $\theta$ is the random variable under consideration (as commented by whuber this letter usually denotes parameters, but I kept with the notation of the question). The integral simply follows from explicitly expressing expectation values as integrals with respect to the probability distribution Apr 14, 2022 at 12:35
• @J.Delaney I'm confused by this part: [𝑔(𝜃;𝑇)−1(𝜃>𝑇)]. Why would we need to subtract in this way? You've set 𝑔(𝜃;𝑇)=1(𝜃>𝑇). So doesn't this just work out to zero? Apr 14, 2022 at 13:14
• Thanks for this response! Your integral for the unbiased estimator is (I believe) what I’m trying to get at—couldn’t get the formulation right, but when you write it that way, seems clear that there’s no “good” estimator that really does what I have in mind (and taking the CDF of the theta_2-theta_1 normal distribution doesn’t work) Apr 17, 2022 at 17:12

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.