A good code for a normal cumulative distribution?

Lucetios · June 10th, 2005, 7:53 pm

I have ued this code. Let me know what you think. Here is the link

JWD · June 11th, 2005, 5:22 pm

Hi Lucetios –The SunPro code at http://www.ks.uiuc.edu/Research/namd/do ... ource.html that you quote claims great accuracy. It uses different approximations for different values of the argument, which is certainly reasonable. However it’s never clear how to verify the coefficients of the various rational functions (Pade approximants) that appear. Just as a sanity check it might be an idea to run a comparison for a sequence of arguments against a standard reference, e.g. Abramowitz and Stegun (AS), equation 26.2.17 (for the related normal integral). The AS equation has a uniform error upper bound 7.5*10^(-8). The SunPro code is supposed to be much more accurate than that, so presumably the difference of the AS code and SunPro’s code (with appropriate change of variable etc.) should be within the same uniform error.On a slightly possibly more fundamental note, typical model errors due to discretization, parameter estimation, etc. can be orders of magnitude greater than even the AS approximation, so it’s hard for me to see the motivation for using a complicated code like SunPro’s. What is your application?Having said this, I am definitely in favor of using different numerical methods for evaluating the bivariate integral depending on the parameters in order to get improved accuracy.--------

Cuchulainn · June 11th, 2005, 7:06 pm

QuoteOriginally posted by: LucetiosI have ued this code. Let me know what you think. Here is the linkLucetios,What are the copyright restrictions on this code?If I make a s/w product can I redistrbute the code or is it for internal use only?Daniel

AVt · June 12th, 2005, 7:48 am

You find good codes at Graeme West accompaning his paper "Better Approximations to Cumulative Normal Functions", they are non-easy 'translations' from Fortran, see his references for 'Genz'.

JWD · June 12th, 2005, 10:50 am

Most of G. West’s attention is focused on the difficulties of numerically evaluating bivariate and higher variate integrals. There, a big problem is well known to be accuracy near the singular points corr = +-1.For this reason as I mentioned before I think it is a good idea to use different numerical methods for different values of the correlations. For the bivariate case, I use Owen’s method (including Borth’s modification), Curnow/Dunnett, and Drezner. The choice of the method for different values of the correlation is taken from the paper of Sowden and Ashford, and then modified further. A list of references I use is as follows:1. Borth, D., A Modification of Owen’s Method for Computing the Bivariate Normal Integral. Applied Statistics Vol 22, No 1, 1973. Pp. 82-85.2. Sowden, R. and Ashford, J., Computation of the Bivariate Normal Integral. Appl. Statistics, Vol 18, 1969, Pp. 169-180.3. Curnow, R. and Dunnett, C., The Numerical Evaluation of Certain Multivariate Normal Integrals. Ann. Math. Statist., 1962. Pp. 571-579.4. Owen, D.B., Tables for Computing Bivariate Normal Probabilities. Ann. Math. Statist., Vol 27, 1956. Pp 1075-1090.5. Drezner, Z., Computation of the Bivariate Normal Integral. Mathematics of Computation, Vol 32, 1978, Pp 277-279.--------

AVt · June 12th, 2005, 8:07 pm

The original Q was for the univariate case if i got it right ...But beyond:The linked sources are ok & fast for any practical situation on a PCwith exactness of about 15 digits, Genz already cares for approachingthe singular case. Graeme's preprint will come up at Wilmott.If one wants better solutions in dim <= 2 one can avoid integrationcompletely and switches to rapidly convergent recursions which arestopped by desired exactness (depending on your environment).

Errrb · June 12th, 2005, 8:56 pm

You can earn a PhD trying to reinvent the wheel and write your own good implementation of the normal CDF. I would suggest to use existing good libraries. Two choices come to mind 1. cephes library, can be found at netlib.org2. gnu scientific library (GSL)hope this helps

Lucetios · June 13th, 2005, 11:54 am

JanDash,I have not verified the accuracy. By, taking the two large and small approximation polynomials I wrote one that was good to 12 decimals. I simply verified it with the actual integral calculated in MatLab. I will do the same thing with this guys and post the result.Cuchulainn,If I remember correctly, I think the waiver is printed on the top of the code. I think you can redistribute it as long as you leave the waiver there.AVt,Are yo usaying that the code is not going to work on UNIX/LINUX?ThanksL.

Cuchulainn · June 13th, 2005, 6:11 pm

Cuchulainn,If I remember correctly, I think the waiver is printed on the top of the code. I think you can redistribute it as long as you leave the waiver there.thxDD

AVt · June 13th, 2005, 6:12 pm

I said it will run on a PC, that is not a question of the OS, just compileit. It is far less sophisticated than Moshier's cephes (which i never getto work with MSVC, but i am a lousy programmer) which Errrb mentions. Forthe GSL the easiest is to use a wrapper, so you just call the lib and getall the other cdf as well. Note that in a commercial sense both of these_ sources _ are not free.And as already mentioned there is a quite different way, a series solutionwith recursions. Find it enclosed for Excel, it also covers dim = 2 , evenfor approaching correlation ~ 1 and might meet JanDash's needs (and no, i donot have such a thing in dim=3, but would switch to integration for that).

JWD · June 20th, 2005, 1:10 pm

What methods do people use for N-variate-integral calculations, besides brute force integration with or without Monte Carlo? I have one reference that, as I remember, gives an explicit solution in terms of an integral over a product of factors; I’ll dig it up. In particular, the concerns for limiting cases of correlations would naturally exist in higher dimensions also. Does anyone know what is done in practice?--------

DavidJN · June 20th, 2005, 4:40 pm

For N-variate-integral calculations you can find a Fortran algorithm in the ACM Transactions on Mathematical Software (TOMS ) and translate as necessary. The one by Drezner (1992) seems robust. It may be algorithm number 725, if memory serves me correct.

quantie · June 20th, 2005, 5:51 pm

QuoteOriginally posted by: JanDashWhat methods do people use for N-variate-integral calculations, besides brute force integration with or without Monte Carlo? I have one reference that, as I remember, gives an explicit solution in terms of an integral over a product of factors; I’ll dig it up. In particular, the concerns for limiting cases of correlations would naturally exist in higher dimensions also. Does anyone know what is done in practice?--------For the trivial zero correlation case it is simply a product of error functions and for non zero correlation there are algorithms like mulnor And from what I gather mathematica uses brute-force numerical integration in the background. I get 8-12 digits precision between mathematica and the following referenceStuart, A. and Ord, J.K. (see pg 515) They provide an expression for orthant probabilities for example, P[X<=0,Y<=0] = (1/4 + ArcSin[rho]/2*pi) and simillar expressions for the trivariate case.QuoteOn a slightly possibly more fundamental note, typical model errors due to discretization, parameter estimation, etc. can be orders of magnitude greater than even the AS approximation, so it’s hard for me to see the motivation for using a complicated code like SunPro’s. What is your application?Agreed, perhaps deep OTM "wing" options.?