C# 中的标准正态分布 z 值函数

Question

我一直在查看 Jeff Atwood 最近在 Alternate Sorting Orders 上发表的博客文章。我试图将帖子中的代码转换为 C#，但遇到了问题。给定标准正态曲线下面积的百分比，我知道.NET 中没有函数会返回 z 值。该算法的推荐值是 95% 和 97.5%，您可以在任何统计书籍的 z 值表中查找。

有谁知道如何为所有 z 值或至少 6 个标准差的平均值实现这样的函数。一种方法是将值硬编码到字典中并使用查找，但必须有一种计算确切值的方法。 我试图解决这个问题是对标准正态曲线函数进行定积分。

y = (1 / (sqrt(2 * PI))) * e^(-(1/2) * x^2)

这给了我两个 x 值之间曲线下的面积，但是我被卡住了……也许我是基础的方式，这不是你会做的吗？

谢谢。

Answer 1

这里有一些用 Python 编写的正态分布的 code，但可以通过添加一些标点符号轻松地将其转换为 C#。大约只有 15 行代码。

Answer 2

这是统计程序 R 中使用的普通分位数 C code 的 C# 翻译。

/// <summary>
/// Quantile function (Inverse CDF) for the normal distribution.
/// </summary>
/// <param name="p">Probability.</param>
/// <param name="mu">Mean of normal distribution.</param>
/// <param name="sigma">Standard deviation of normal distribution.</param>
/// <param name="lower_tail">If true, probability is P[X <= x], otherwise P[X > x].</param>
/// <param name="log_p">If true, probabilities are given as log(p).</param>
/// <returns>P[X <= x] where x ~ N(mu,sigma^2)</returns>
/// <remarks>See https://svn.r-project.org/R/trunk/src/nmath/qnorm.c</remarks>
public static double QNorm(double p, double mu, double sigma, bool lower_tail, bool log_p)
{
  if (double.IsNaN(p) || double.IsNaN(mu) || double.IsNaN(sigma)) return (p + mu + sigma);
  double ans;
  bool isBoundaryCase = R_Q_P01_boundaries(p, double.NegativeInfinity, double.PositiveInfinity, lower_tail, log_p, out ans);
  if (isBoundaryCase) return (ans);
  if (sigma < 0) return (double.NaN);
  if (sigma == 0) return (mu);

  double p_ = R_DT_qIv(p, lower_tail, log_p);
  double q = p_ - 0.5;
  double r, val;

  if (Math.Abs(q) <= 0.425)  // 0.075 <= p <= 0.925
  {
    r = .180625 - q * q;
    val = q * (((((((r * 2509.0809287301226727 +
               33430.575583588128105) * r + 67265.770927008700853) * r +
             45921.953931549871457) * r + 13731.693765509461125) * r +
           1971.5909503065514427) * r + 133.14166789178437745) * r +
         3.387132872796366608)
    / (((((((r * 5226.495278852854561 +
             28729.085735721942674) * r + 39307.89580009271061) * r +
           21213.794301586595867) * r + 5394.1960214247511077) * r +
         687.1870074920579083) * r + 42.313330701600911252) * r + 1.0);
  }
  else
  {
    r = q > 0 ? R_DT_CIv(p, lower_tail, log_p) : p_;
    r = Math.Sqrt(-((log_p && ((lower_tail && q <= 0) || (!lower_tail && q > 0))) ? p : Math.Log(r)));

    if (r <= 5)              // <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11
    {
      r -= 1.6;
      val = (((((((r * 7.7454501427834140764e-4 +
              .0227238449892691845833) * r + .24178072517745061177) *
            r + 1.27045825245236838258) * r +
           3.64784832476320460504) * r + 5.7694972214606914055) *
         r + 4.6303378461565452959) * r +
        1.42343711074968357734)
       / (((((((r *
                1.05075007164441684324e-9 + 5.475938084995344946e-4) *
               r + .0151986665636164571966) * r +
              .14810397642748007459) * r + .68976733498510000455) *
            r + 1.6763848301838038494) * r +
           2.05319162663775882187) * r + 1.0);
    }
    else                     // very close to  0 or 1 
    {
      r -= 5.0;
      val = (((((((r * 2.01033439929228813265e-7 +
              2.71155556874348757815e-5) * r +
             .0012426609473880784386) * r + .026532189526576123093) *
           r + .29656057182850489123) * r +
          1.7848265399172913358) * r + 5.4637849111641143699) *
        r + 6.6579046435011037772)
       / (((((((r *
                2.04426310338993978564e-15 + 1.4215117583164458887e-7) *
               r + 1.8463183175100546818e-5) * r +
              7.868691311456132591e-4) * r + .0148753612908506148525)
            * r + .13692988092273580531) * r +
           .59983220655588793769) * r + 1.0);
    }
    if (q < 0.0) val = -val;
  }

  return (mu + sigma * val);
}

一些辅助方法：

private static bool R_Q_P01_boundaries(double p, double _LEFT_, double _RIGHT_, bool lower_tail, bool log_p, out double ans)
{
  if (log_p)
  {
    if (p > 0.0)
    {
      ans = double.NaN;
      return (true);
    }
    if (p == 0.0)
    {
      ans = lower_tail ? _RIGHT_ : _LEFT_;
      return (true);
    }
    if (p == double.NegativeInfinity)
    {
      ans = lower_tail ? _LEFT_ : _RIGHT_;
      return (true);
    }
  }
  else
  {
    if (p < 0.0 || p > 1.0)
    {
      ans = double.NaN;
      return (true);
    }
    if (p == 0.0)
    {
      ans = lower_tail ? _LEFT_ : _RIGHT_;
      return (true);
    }
    if (p == 1.0)
    {
      ans = lower_tail ? _RIGHT_ : _LEFT_;
      return (true);
    }
  }
  ans = double.NaN;
  return (false);
}

private static double R_DT_qIv(double p, bool lower_tail, bool log_p)
{
  return (log_p ? (lower_tail ? Math.Exp(p) : -ExpM1(p)) : R_D_Lval(p, lower_tail));
}

private static double R_DT_CIv(double p, bool lower_tail, bool log_p)
{
  return (log_p ? (lower_tail ? -ExpM1(p) : Math.Exp(p)) : R_D_Cval(p, lower_tail));
}

private static double R_D_Lval(double p, bool lower_tail) 
{
  return lower_tail ? p : 0.5 - p + 0.5; 
} 

private static double R_D_Cval(double p, bool lower_tail) 
{ 
  return lower_tail ? 0.5 - p + 0.5 : p;
}
private static double ExpM1(double x) 
{
  if (Math.Abs(x) < 1e-5)
     return x + 0.5 * x * x;
  else
     return Math.Exp(x) - 1.0;
 }

在您的情况下，您需要 mu=0.0、sigma=1.0、lower_tail=true、log_p=false。

Answer 3

查找error function 的实现。 ... 书籍中的所有经典数字食谱中都有一个。

Answer 4

对于较新版本的 MathNet

    //standard normal cumulative distribution function
    static double F(double x)
    {
        MathNet.Numerics.Distributions.Normal result = new MathNet.Numerics.Distributions.Normal();
        return result.CumulativeDistribution(x);
    }

Answer 5

如果您使用的是 Math.Net，那么您可以使用 InverseCumulativeDistribution 函数。

所以如果你想要 Z 值，其中 80% 的标准正态曲线被覆盖，代码看起来像这样

var curve = new MathNet.Numerics.Distributions.Normal();
var z_value = curve.InverseCumulativeDistribution(0.8);
Console.WriteLine(z_value);

确保您已安装MathNet.Numerics NuGet 包。