RE: st: Numerical accuracy of standard normal distribution

[Georges Casamatta <georges.casamatta@ecopa.com>]

Yes this helps! 

Thanks Marteen, 

Georges



-----Message d'origine-----
De : owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] De la part de Maarten Buis
Envoyé : mercredi 11 septembre 2013 09:55
À : statalist@hsphsun2.harvard.edu
Objet : Re: st: Numerical accuracy of standard normal distribution

On Wed, Sep 11, 2013 at 8:49 AM, Georges Casamatta wrote:
> I am coding a program for estimating a probit model using the IRLS algorithm. For this, I use the cumulative distribution function of the standard normal distribution in mata. I am surprised that it is equal to 1 for a value as small as 7 (the same kind of comment of course also applies to the density function). Isn't there a way to improve the numerical accuracy of this function ?

As you can see by typing -di %21x normal(7)-, the -normal()- function does not return an exact 1, but something very close. See:
<http://blog.stata.com/2011/02/02/how-to-read-the-percent-21x-format/>
and <http://blog.stata.com/2011/02/10/how-to-read-the-percent-21x-format-part-2/>
This should not come as a surprise, because normal(7) is very very very close to 1, See:

twoway function y= normal(x), range(-10 10) xline(7)

When using -normal()- in a program the precision problem that often comes up occurs when one computes:

1-normal(x)

Instead you should code that as the mathametically equivalent but numerically more stable form:

normal(-x)

For more on this see: William Gould (2006) "Mata matters: Precision".
The Stata Journal, 6(4): 550-560.
<http://www.stata-journal.com/article.html?article=pr0025>

As an aside: In that article Bill starts with an example of what can go wrong by showing a graph of the function normalden(x)/(1-normal(x)), which after approximately x = 7 starts to oscillate. Some time ago I had a discussion with some R users and I experimented with that same function in R and R showed exactly the same oscillation pattern as Stata. Not only that, the numbers were exactly the same. That to me suggest that Stata and R use exactly the same algorithm. The solution for both programs is ofcourse to code that function as normalden(x)/normal(-x).

Hope this helps,
Maarten

---------------------------------
Maarten L. Buis
WZB
Reichpietschufer 50
10785 Berlin
Germany

http://www.maartenbuis.nl
---------------------------------

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/