Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Given a small sample from a normally-distributed population, how do I calculate the confidence that a specified percentage of the population is within some bounds [A,B]?

To make it concrete, if I get a sample of [50,51,52] from a normally-distributed set, I should be able to calculate a fairly high confidence that 50% of the population lies within the range of 0-100, even with such a small sample.

This is certainly related to the "tolerance interval", but differs in an important way. In all of the examples I can find for tolerance intervals, the required percentile and confidence is given, and the interval is found. In my problem, the interval and percentile are given, and I need to find the confidence.

The relevant equation is this one: (Guttman 1970)

$$1 - \gamma = P\left[P(X \geqq t_0) \geqq 1 - p\right] = P\left[T_{n-1}^*(\sqrt n z_p) \leqq \sqrt n K\right]$$

With definitions:

  • $1 - \gamma$ is the confidence
  • $n$ is the number of samples
  • $100p$ is the percentage of the population required to be within the interval, as estimated from the sample mean and sample variance.
  • $t_0 = x - K_{S, z_p}$ is the $(1 - p) 100$th percentile of the standard normal distribution
  • $T_v^*(\delta)$ is the noncentral Student’s t distribution with $v$ degrees of freedom and noncentrality parameter $\delta$.

This solves the one-sided problem, but I'm having trouble extending this to the two-sided problem. In confidence-interval land, I'd use the fact that $P(t_1 \leqq X \leqq t_2) = 1 - P(t_1 \gt X) - P(X \gt t_2)$, to break this into two one-sided problems, but in tolerance-interval land I need to relate these back to the confidence ($1-\gamma$), and I don't see how.

$$1 - \gamma = P\left[P(t_1 \geqq X \geqq t_2) \geqq 1 - p\right] = ??? $$

If I attempt to turn this into two one-sided problems:

$$1 - \gamma = P\left[1 - P(t_1 \lt X) - P(X \lt t_2) \geqq 1 - p\right] = ??? $$

And I'm utterly stuck there. I don't see how to relate this back to the one-sided tolerance interval solution.


I'm not certain this is useful for people to understand the question, but it might, so I'm putting it in this addenda.

In scipy, I'm able to pretty easily calculate $K$ given $p$ $\gamma$ and $n$ as:

def K(p, gamma, n):
    from scipy import stats
    return stats.nct.ppf(1-gamma, n-1, sqrt(n) * stats.norm.ppf(1-p)) / sqrt(n)

I'm also able to find $\gamma$ given $K$ $p$ and $n$ as:

def gamma(p, n, K):
    from scipy import stats                                                                                              
    z_p = stats.norm.ppf(1-p)
    return 1 - stats.nct.cdf(sqrt(n) * K, n-1, sqrt(n) * z_p)

Much less important, but is this a valid simplification of the Guttman's formula?

$$1 - \gamma = P\left[P(X \geqq t_0) \geqq 1 - p\right] = P\left[T_{n-1}^*(\sqrt n z_p) \leqq \sqrt n K\right]$$ $$\gamma = P\left[P(X \geqq t_0) \lt 1 - p\right] = P\left[T_{n-1}^*(\sqrt n z_p) \gt \sqrt n K\right]$$ $$\gamma = P\left[P(X \lt t_0) \lt p\right] = P\left[T_{n-1}^*(\sqrt n z_p) \gt \sqrt n K\right]$$

If so, this form seems way easier to understand, to me.

share|improve this question
    
Could you rephrase your "concrete question" in more mathematical terms ? Do you want to estimate the probability that the unknown normal distribution falls within a given interval ? –  Stéphane Laurent Feb 13 '13 at 22:48
    
The probability that any distribution "falls within" a given range will be 1, since all normal distributions have the reals as their domain. I'm trying to estimate the probability that any one sample in the population falls within a specified acceptable range. This is what I try to express by "percentage of the population". –  bukzor Feb 14 '13 at 2:57
    
I try to rephrase my question: do you want to estimate $\Pr(X \in [a,b])$ for a given interval $[a,b]$ ? (this probability is not $1$...) –  Stéphane Laurent Feb 14 '13 at 7:14
    
@StéphaneLaurent Yes, I believe so. I take that to mean the probibility that the random variable X will be between a and b. Seem right. Secondarily, I believe I need the confidence of this estimate (I think). –  bukzor Feb 14 '13 at 7:19
2  
I'll try to help you when I'll be less busy. But you can also take a look (with the help of Google) at the notion of tolerance intervals. I think this is related to your aim. –  Stéphane Laurent Feb 14 '13 at 16:34
show 6 more comments

3 Answers

One possibility: you have a certain number of successes and failures, and you want to guess the percentage of the underlying population. You calculate the Agresti-Coull interval (see Brown, Cai, & DasGupta 2001) with the required parameters.

On the bright side this is easy and makes no assumptions about the type of distribution. On the downside it doesn't take advantage of the extra information about where the samples in the interval fall -- if you know something about how they are distributed this might get you a more precise estimate based on the data you have.

Edit: for a brief overview see http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Agresti-Coull_Interval

share|improve this answer
    
Yes but I think he wants to take the advantage of the Gaussian distributional assumption. –  Stéphane Laurent May 23 '13 at 18:35
    
I don't believe this will ever accomplish the scenario in my second paragraph: "if I get a sample of [50,51,52] from a normally-distributed set, I should be able to calculate a fairly high confidence that 50% of the population lies within the range of 0-100, even with such a small sample." –  bukzor May 23 '13 at 18:48
add comment

The following R code simulates the posterior distribution of the Bayesian Gaussian model $y_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2)$ starting with the Jeffreys prior and returns the $100(1-\alpha)\%$ equi-tailed two-sided confidence interval about the probability $\Pr({\cal N}(\mu, \sigma^2) \in [a,b])$.

# confidence inteval about  Pr(a < N(mu,sigma^2) < b)
bounds <- function(y, a, b, alpha=.05, nsims=10000){
        n <- length(y)
    ssq <- (n-1)*var(y)
    rho <- rchisq(nsims,n)/ssq
    mu <- rnorm(nsims, mean(y), 1/sqrt(n*rho))
    p <- pnorm(b,mu,1/sqrt(rho))-pnorm(a,mu,1/sqrt(rho))
return(quantile(p, c(alpha/2,1-alpha/2)))
} 

Example:

> n <- 100
> y <- rnorm(n,0,1)
> bounds(y,a=-1.96,b=1.96)
     2.5%     97.5% 
0.9146264 0.9773575 

There's no theoretical guarantee that the confidence level is well controlled (that is, for example, the effective confidence level of the confidence interval with $95\%$ nominal confidence level is not exactly $95\%$). You should use simulations to check, but I think the effective confidence level is rather close to the nominal one.

share|improve this answer
    
Thanks! Where is a good place to find the definition of those functions? Also I don't quite understand the final output. What exactly are those four numbers? –  bukzor May 25 '13 at 0:54
    
I'll come back later (or tomorrow) to provide more details. Actually there's only 2 numbers in the output : $0.9146$ is the lower $2.5\%$-confidence bound and $0.977$ is the upper $97.5\%$-confidence bound, thus $[0.9146,0.977]$ is a two-sided $95\%$-confidence interval (called equi-tailed). –  Stéphane Laurent May 25 '13 at 5:48
    
Could you please define these? yi ~ iid ssq rho –  bukzor Jun 3 '13 at 17:08
    
If I understand this solution correctly, you're using simulated random outcomes to avoid integration across the estimated parameters. This adds another layer of uncertainty which doesn't seem entirely necessary. Is there a similar solution which uses integration rather than simulation? –  bukzor Jun 3 '13 at 17:10
    
@bukzor Type something like "conjugate bayesian analysis gaussian distribution" in Google. The posterior distribution of $\sigma^2$ is an inverse-Gamma and the conditional posterior distribution of $\mu$ given $\sigma$ is Gaussian. Thus the joint distribution of $\mu$ and $\sigma$ is not complicated but we need the inverse cumulative distribution of the probability, which is quite more complicated. –  Stéphane Laurent Jun 3 '13 at 17:34
show 2 more comments

I would suggest using bootstraping. It is especially very useful for small sample size case. You can get more information about it here: http://en.wikipedia.org/wiki/Bootstrapping_(statistics)

share|improve this answer
    
Will bootstrapping be able to accomplish the scenario in my second paragraph? I've looked at bootstrapping before, but all the descriptions were too vague for me to create a concrete implementation. Do you know of any simple implementations I could look at? –  bukzor Jun 3 '13 at 17:06
    
Boostraping is only asymptotically valid (i.e. for large sample sizes). –  Stéphane Laurent Jun 3 '13 at 17:35
    
@bukzor you can check here: photon.isy.liu.se:81/~fredrik/spcourse/Zoubir.pdf and eprints.qut.edu.au/14131/1/14131.pdf it has many applications including the composite hypothesis testing. It is used primarily when the sample size is small. It is computationally expensive for real time applications. Please have a look at the references. If you are talking about small sample size and lack of samples to build you statistical model, then atually you are talking about bootstrap. –  Seyhmus Güngören Jun 3 '13 at 20:48
    
@StéphaneLaurent what do you mean? –  Seyhmus Güngören Jun 3 '13 at 20:50
    
@SeyhmusGüngören Bootstrap is not valid for small sample sizes. For instance the effective coverage of the bootstrap $95\%$-confidence interval is not $95\%$, and generally it is really smaller. Mathematical results about bootstrap are asymptotic results ($n \to \infty$). en.wikipedia.org/wiki/Bootstrapping_(statistics)#Disadvantages –  Stéphane Laurent Jun 4 '13 at 5:33
show 11 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.