Re: st: Numerical accuracy of standard normal distribution

[Maarten Buis <maartenlbuis@gmail.com>]

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/