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?

$\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$ 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$ 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
    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$ 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$ Feb 15, 2012 at 0:27

5 Answers 5

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.

$\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
    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$ Oct 9, 2010 at 14:24
22
$\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
$\endgroup$
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$ 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
$\endgroup$
2
  • 1
    $\begingroup$ Quantile vs. percentile (it's merely a matter of terminology), j.mp/dsYz9z. $\endgroup$
    – chl
    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
    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).

$\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)$$

$\endgroup$
4
  • 2
    $\begingroup$ isn't this more of a comment than an answer? $\endgroup$
    – Macro
    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$ 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
    Mar 5, 2012 at 22:31
  • $\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
    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.