14
$\begingroup$

See this Wikipedia page: Binomial proportion confidence interval.

To get the Agresti-Coull Interval, one needs to calculate a percentile of the normal distribution, called $z$. How do I calculate the percentile? Is there a ready-made function that does this in Wolfram Mathematica and/or Python/NumPy/SciPy?

Improve this question
$\endgroup$
17
  • 1
    $\begingroup$ The integral expression in the "normal cdf I got exactly from Wiki" is unfortunately off by a factor of $1/\sqrt{\pi}$. There is no known exact formula for the normal cdf or its inverse using a finite number of terms involving standard functions ($\exp, \log, \sin \cos$ etc) but both the normal cdf and its inverse have been studied a lot and approximate formulas for both are programmed into many calculator, spreadsheets, not to mention statistical packages. I am not familiar with R but I would be astounded if it did not have what you are looking for built in already. $\endgroup$
    – Dilip Sarwate
    Commented Feb 14, 2012 at 20:52
  • $\begingroup$ @DilipSarwate, it's fixed! I am doing this using inverse tranformation, also "not allowed" to use too much built in. It's for the sake of learning I suppose. $\endgroup$
    – user1061210
    Commented Feb 14, 2012 at 21:26
  • 1
    $\begingroup$ @Dilip: Not only is there no known exact formula, better yet, it is known that no such formula can exist! $\endgroup$
    – cardinal
    Commented Feb 14, 2012 at 21:30
  • 1
    $\begingroup$ The Box-Muller method generates samples from a joint distribution of independent standard normal random variables. So histograms of the values generated will resemble standard normal distributions. But the Box-Muller method is not a method for computing values of $\Phi(x)$ except incidentally as in "I generated $10^4$ standard normal samples of which $8401$ has value $1$ or less, and so $\Phi(1) \approx 0.8401$, and $\Phi^{-1}(0.8401) \approx 1$. $\endgroup$
    – Dilip Sarwate
    Commented Feb 14, 2012 at 22:29
  • 1
    $\begingroup$ I just chose $8401$ as an example of the kinds of numbers you might expect. $\Phi(1) = 0.8413\ldots$ and so if you generate $10^4$ samples of a standard normal distribution, you should expect close to $8413$ of the $10000$ samples to have value $\leq 1$. You are implementing the Box-Muller method correctly, but are not understanding the results that you are getting and are not relating them to the cdf etc. $\endgroup$
    – Dilip Sarwate
    Commented Feb 15, 2012 at 0:27

5 Answers 5

Reset to default
3
$\begingroup$

For Mathematica $VersionNumber > 5 you can use

Quantile[NormalDistribution[μ, σ], 100 q]

for the q-th percentile.

Otherwise, you have to load the appropriate Statistics package first.

Improve this answer
$\endgroup$
2
  • $\begingroup$ (I have version 7.) I have no problem loading the Statistics package. But what's the function in there called? Because I get the impression that this Quantile line will do the calculation manually instead of using a formula. $\endgroup$
    – Ram Rachum
    Commented Oct 9, 2010 at 14:13
  • $\begingroup$ Evaluate it with symbolic parameters (i.e. don't assign values to mu, sigma, and q); you should get an expression involving the inverse error function. $\endgroup$
    – J. M. is not a statistician
    Commented Oct 9, 2010 at 14:24
23
$\begingroup$

John Cook's page, Distributions in Scipy, is a good reference for this type of stuff:

In [15]: import scipy.stats

In [16]: scipy.stats.norm.ppf(0.975)
Out[16]: 1.959963984540054
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ This is what I was looking for. BTW, if you want to calculate the number of standard deviations around the mean which capture p% of probability, you can do it as follows: sigmas_around_mean = lambda p: (norm.ppf(0.5 + p/2) - norm.ppf(0.5 - p/2))/2 E.g. You should get sigmas_around_mean(0.95) = 1.96. This is useful when you want to know how many stds around the mean capture, say, 50% of the probability mass. $\endgroup$
    – Abhishek Divekar
    Commented Aug 26, 2023 at 14:18
4
$\begingroup$

Well, you didn't ask about R, but in R you do it using ?qnorm

(It's actually the quantile, not the percentile, or so I believe)

> qnorm(.5)
[1] 0
> qnorm(.95)
[1] 1.644854
Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Quantile vs. percentile (it's merely a matter of terminology), j.mp/dsYz9z. $\endgroup$
    – chl
    Commented Oct 9, 2010 at 14:22
  • 1
    $\begingroup$ While we are in, in R Wald-adjusted CIs (e.g. Agresti-Coull) are available in the PropCIs package. Wilson's method is the default in Hmisc::binconf (as suggested by Agresti and Coull). $\endgroup$
    – chl
    Commented Oct 9, 2010 at 14:36
4
$\begingroup$

In Python, you can use the stats module from the scipy package (look for cdf(), as in the following example).

(It seems the transcendantal package also includes usual cumulative distributions).

Improve this answer
$\endgroup$
1
$\begingroup$

You can use the inverse erf function, which is available in MatLab and Mathematica, for instance.

For the normal CDF, starting from

$$y=\Phi\left(x\right)=\frac{1}{2}\left[1+\text{erf}\left(\frac{x}{\sqrt{2}}\right)\right]$$

We get

$$x=\sqrt{2}\ \text{erf}^{-1}\left(2y-1\right)$$

For the log-normal CDF, starting from

$$y=F_{x}(x;\mu,\sigma)=\frac{1}{2}\text{erfc}\left(\frac{-\log x-\mu}{\sigma\sqrt{2}}\right)$$

We get

$$-\log \left(x\right)=\mu+\sigma\sqrt{2}\ \text{erfc}^{-1}\left(2y\right)$$

Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ isn't this more of a comment than an answer? $\endgroup$
    – Macro
    Commented Feb 15, 2012 at 22:18
  • $\begingroup$ My idea was that if you have inverses for the erf and erfc functions, then the problem is solved. MatLab, for instance, has such preprogrammed functions. $\endgroup$
    – Jean-Victor Côté
    Commented Mar 5, 2012 at 19:31
  • $\begingroup$ @Jean-VictorCôté Please, develop your ideas in your reply. Otherwise, it merely looks like a comment as suggested above. $\endgroup$
    – chl
    Commented Mar 5, 2012 at 22:31
  • 1
    $\begingroup$ The lognormal calculation doesn't look right. After all, its inverse CDF should be identical to the inverse CDF for the normal apart for the use of $\log(x)$ instead of $x$. $\endgroup$
    – whuber ♦
    Commented Mar 6, 2012 at 19:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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