累積常態分配反函數

前一篇 常態分配數值逼近法 曾經提到怎麼運用數值分析求得 CDF

但是常常我們不只是需要使用 CDF, 我們可能還需要取得標準常態分配 N(x) 中的 x 變數其值。

在這邊感謝同事提供他的資料給我，才得知原來已經有好幾個人針對反函數求解。

查找了一下文章後發現 Quantitative Finance Collector這篇提及了兩種作法，

分別是 Boris Moro 以及 Peter J Achlam所發表的理論。

在這邊我實作了 Boris Moro 的方法，其程式碼如下所示：

        const double c1 = 0.337475482272615;
        const double c2 = 0.976169019091719;
        const double c3 = 0.160797971491821;
        const double c4 = 2.76438810333863E-02;
        const double c5 = 3.8405729373609E-03;
        const double c6 = 3.951896511919E-04;
        const double c7 = 3.21767881768E-05;
        const double c8 = 2.888167364E-07;
        const double c9 = 3.960315187E-07;
        /// 

        /// Calculates the Normal Standard numbers given u, the associated uniform number (0, 1)
        /// 
        /// 
        /// 
        public static double Inverse(double u)
        {
            double[] a = new double[4] { 2.50662823884, -18.61500062529, 41.39119773534, -25.44106049637 };
            double[] b = new double[4] { -8.4735109309, 23.08336743743, -21.06224101826, 3.13082909833 };

            double x = u - 0.5;
            double r;
            if (Math.Abs(x) < 0.42)
            {
                r = x * x;
                r = x * (((a[3] * r + a[2]) * r + a[1]) * r + a[0]) / ((((b[3] * r + b[2]) * r + b[1]) * r + b[0]) * r + 1);
            }
            else
            {
                r = Math.Log(-Math.Log(1 - u));
                if (x <= 0)
                {
                    r = -r;
                }
                r = c1 + r * (c2 + r * (c3 + r * (c4 + r * (c5 + r * (c6 + r * (c7 + r * (c8 + r * c9)))))));
                if (x < 0)
                {
                    r = -r;
                }
            }
            return r;
        }