Tell me more ×
Cross Validated is a question and answer site for statisticians, data analysts, data miners and data visualization experts. It's 100% free, no registration required.
up vote 6 down vote favorite

See this Wikipedia page:

http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Agresti-Coull_Interval

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

share|improve this question
1  
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. – Dilip Sarwate Feb 14 '12 at 20:52
@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. – user1061210 Feb 14 '12 at 21:26
1  
@Dilip: Not only is there no known exact formula, better yet, it is known that no such formula can exist! – cardinal Feb 14 '12 at 21:30
1  
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$. – Dilip Sarwate Feb 14 '12 at 22:29
1  
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. – Dilip Sarwate Feb 15 '12 at 0:27
show 12 more comments

5 Answers

up vote 5 down vote accepted

For Mathematica, $VersionNumber > 5 you can use

Quantile[NormalDistribution[mu,sigma], 100 q]

for the q-th percentile.

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

share|improve this answer
(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. – Ram Rachum Oct 9 '10 at 14:13
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. – J. M. Oct 9 '10 at 14:24
up vote 5 down vote

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
share|improve this answer
Thanks for adding a real illustration with ppf. – chl♦ Oct 9 '10 at 16:22
up vote 4 down vote

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
share|improve this answer
1  
Quantile vs. percentile (it's merely a matter of terminology), j.mp/dsYz9z. – chl♦ Oct 9 '10 at 14:22
1  
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). – chl♦ Oct 9 '10 at 14:36
Thanks for the comments chl – Tal Galili Oct 10 '10 at 4:11
up vote 2 down vote

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).

share|improve this answer
up vote 0 down vote

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

share|improve this answer
1  
isn't this more of a comment than an answer? – Macro Feb 15 '12 at 22:18
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. – Jean-Victor Côté Mar 5 '12 at 19:31
@Jean-VictorCôté Please, develop your ideas in your reply. Otherwise, it merely looks like a comment as suggested above. – chl♦ Mar 5 '12 at 22:31
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$. – whuber♦ Mar 6 '12 at 19:15

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.