How do I know whether future samples will remain below a threshold?

Question

I have N samples of quantity x (N=40 or 50 or so; I illustrate 7 of them here).

0.529578449
0.24137483
0.715602119
0.541831981
0.426401788
0.426117433
0.678243369
....

To these I can reasonably fit e.g. a normal distribution with, say, mean 0.51 and standard deviation 0.16.

The question I have to address is: How confidently can we assert that >99% of future samples of x will remain below a given threshold, e.g. 1.0?

If I just read the value of the fitted distribution mentioned above, it says that 99.9% will be below 1.0. So far so good.

But how well is that fitted distribution likely to represent the underlying distribution? And how do I translate this into a measure of confidence?

What would be a relevant test/method to address this?

Answer 1

You want to assess whether the probability $\Pr(X>1)$ is significantly higher than $99\%$.

To do so, you can derive a confidence interval about $\Pr(X>1)$. For example you can get such a confidence interval using the Bayesian approach with the Jeffreys prior.

Another way is to use a lower tolerance limit. If the lower $(1-\alpha, p=99\%)$-tolerance limit is higher than $1$, then the probability $\Pr(X>1)$ is significantly higher than $99\%$ at the $\alpha$-level of significance.

Example:

> # simulated sample
> set.seed(666)
> y <- rnorm(40, mean=5, sd=1)
> # tolerance limit
> library(tolerance)
> normtol.int(y, alpha=0.05, P=0.99, side=1)
  alpha    P    x.bar 1-sided.lower 1-sided.upper
1  0.05 0.99 4.874011      1.353383      8.394639

The lower tolerance limit is $\approx 1.35 > 1$, then $\Pr(X>1)$ is significantly higher than $99\%$ at the significance level $\alpha=5\%$.

Using the Jeffreys Bayesian approach:

> Jeffreys <- function(y, nsims=100000){
+   n <- length(y)
+   sigma <- sqrt(c(crossprod(y-mean(y)))/rchisq(nsims,n))
+   mu <- rnorm(nsims, mean(y), sigma/sqrt(n))
+   list(mu=mu, sigma=sigma)
+ }
> # posterior sampling of Pr(Y>1)
> nsims <- 100000
> sims_musigma <- Jeffreys(y, nsims)
> sims_pr <- numeric(nsims)
> for(i in 1:nsims){
+   sims_pr[i] <- 1 - pnorm(1, mean=sims_musigma$mu[i], sd=sims_musigma$sigma[i])
+ }
> # lower confidence bound of Pr(Y>1)
> quantile(sims_pr, 0.05)
       5% 
0.9954999 

The lower $95\%$-confidence bound of $\Pr(X>1)$ is $\approx 99.5\%$, then $\Pr(X>1)$ is significantly higher than $99\%$ at the significance level $\alpha=5\%$.

If you don't like the Jeffreys approach, you can use these approximate confidence bounds of $\Pr(X>q)$:

• lower bound: $1 - \Phi\left[\frac{q-\hat\mu}{\hat\sigma}\left(1-\Phi^{-1}(1-\alpha)\sqrt{\dfrac{1}{n{\left(\frac{q-\hat\mu}{\hat\sigma}\right)}^2}+\dfrac{1}{2(n-1)}}\right) \right]$

• upper bound: $1 - \Phi\left[\frac{q-\hat\mu}{\hat\sigma}\left(1+\Phi^{-1}(1-\alpha)\sqrt{\dfrac{1}{n{\left(\frac{q-\hat\mu}{\hat\sigma}\right)}^2}+\dfrac{1}{2(n-1)}}\right) \right]$

sources:

• Bissell, A. F. (1990), "How Reliable Is Your Capability Index?" Applied Statistics, 30, 331 - 340.

• Kushler, R. H. and Hurley, P. (1992), "Confidence Bounds for Capability Indices," Journal of Quality Technology, 24, 188 - 195.

The lower bound is similar to the previous one:

> alpha <- 5/100
> n <- length(y)
> q <- 1
> 1 - pnorm((q-mean(y))/sd(y) * (1-qnorm(1-alpha)*sqrt(1/n/((q-mean(y))/sd(y))^2 + 1/2/(n-1))))
[1] 0.9950559