標準常態分佈函數和標準常態分佈分位元函數

標準常態分佈函數和分位元函數的數值演算法可參考高惠璿編著的《統計計算》

以下是C#版本的實現代碼：

/// <summary>
/// 標準常態分佈函數Phi(x)
/// </summary>
/// <param name="x">隨機變數</param>
/// <returns>標準常態分佈機率</returns>
static double NormDistFunc(double x)
{
  double x0 = (x >= 0 ? x : -x);
  double[] b = { 0.196854, 0.115194, 0.000344, 0.019527 };
  double erf = 0;
  for (int i = 1; i <= 4; i++)
  {
    erf += b[i - 1] * Math.Pow(x0, i);
  }
  erf = 1 - Math.Pow(1.0 + erf, -4);
  double phi = (x >= 0 ? 0.5 * (1 + erf) : 0.5 * (1 - erf));
  return phi;
}
/// <summary>
/// 標準常態分佈函數分位元函數
/// </summary>
/// <param name="p">機率</param>
/// <returns>分位元</returns>
static double NormDistributionQuantile(double p)
{
  Debug.Assert((0 < p) && (p < 1));
  if (p == 0.5)
    return 0;
  double[] b ={0.1570796288E1,   0.3706987906E-1,
              -0.8364353589E-3, -0.2250947176E-3,
                0.6841218299E-5,  0.5824238515E-5,
              -0.1045274970E-5,  0.8360937017E-7,
              -0.3231081277E-8,  0.3657763036E-10,
                0.6936233982E-12};
  double alpha = 0;
  if ((0 < p) && (p < 0.5))
    alpha = p;
  else if ((0.5 < p) && (p < 1))
    alpha = 1 - p;
  double y = -Math.Log(4 * alpha * (1 - alpha));
  double u = 0;
#if USE_TODA_FORMULA
  //Toda近似公式，最大誤差1.2e-8
  for (int i = 0; i < b.Length; i++)
  {
    u += b[i] * Math.Pow(y, i);
  }
  u = Math.Sqrt(y * u);
#else
  //山內近似公式，最大誤差4.9e-4
  u = Math.Sqrt(y * (2.0611786 - 5.7262204 / (y + 11.640595)));
#endif
  double up = 0;
  if ((0 < p) && (p < 0.5))
    up = -u;
  else if ((0.5 < p) && (p < 1))
    up = u;
  return up;
}
/// <summary>
/// 標準常態分佈機率密度函數
/// </summary>
/// <param name="x">隨機變數</param>
/// <returns>機率密度</returns>
static double NormDensityFunc(double x)
{
  return Math.Exp(-x * x * 0.5) / Math.Sqrt(2 * Math.PI);
}