High-accuracy Complementary normal distribution function. - CSharp System

CSharp examples for System:Math Number

  1. HOME
  2. CSharp
  3. System
  4. Math Number

Description

High-accuracy Complementary normal distribution function.

Demo Code

// Copyright (C) 2013, Alan Genz,  All rights reserved.               
using System;//from  w w  w .  ja v  a2 s  . c  om

public class Main{
        /// <summary>
        ///     High-accuracy Complementary normal distribution function.
        /// </summary>
        /// <remarks>
        ///     <para>
        ///         This function uses 9 tabled values to provide tail values of the
        ///         normal distribution, also known as complementary Phi, with an
        ///         absolute error of 1e-14 ~ 1e-16.
        ///     </para>
        ///     References:
        ///     - George Marsaglia, Evaluating the Normal Distribution, 2004.
        ///     Available in: http://www.jstatsoft.org/v11/a05/paper
        /// </remarks>
        /// <returns>
        ///     The area under the Gaussian p.d.f. integrated
        ///     from the given value to positive infinity.
        /// </returns>
        public static double HighAccuracyComplemented(double x)
        {
            double[] R =
            {
                1.25331413731550025, 0.421369229288054473, 0.236652382913560671,
                0.162377660896867462, 0.123131963257932296, 0.0990285964717319214,
                0.0827662865013691773, 0.0710695805388521071, 0.0622586659950261958
            };

            var j = (int) (0.5*(Math.Abs(x) + 1));

            if (j >= R.Length)
            {
                return x > 0 ? 0 : 1;
            }

            double a = R[j];
            double z = 2*j;
            double b = a*z - 1.0;

            double h = Math.Abs(x) - z;
            double q = h*h;
            double pwr = 1;

            double sum = a + h*b;
            double term = a;


            for (int i = 2; sum != term; i += 2)
            {
                term = sum;

                a = (a + z*b)/(i);
                b = (b + z*a)/(i + 1);
                pwr *= q;

                sum = term + pwr*(a + h*b);
            }

            sum *= Math.Exp(-0.5*x*x - 0.91893853320467274178);

            return (x >= 0) ? sum : (1.0 - sum);
        }
}

Related Tutorials