标准正态分布z值函数在C#
我一直在看杰夫阿特伍德最近的博客文章Alternate Sorting Orders。我试图将帖子中的代码转换为C#,但我遇到了一个问题。 .NET中没有函数,我知道这将返回z值,给定标准正态曲线下面积的百分比。用于算法的建议值分别为95%和97.5%,您可以在任何统计手册的z值表中查找。标准正态分布z值函数在C#
有没有人知道如何实现这样一个函数的所有z值或至少6个标准偏差的平均值。一种方法是将这些值硬编码成字典并使用查找,但必须有一种计算确切值的方法。 我试图解决这个问题是对标准正态曲线函数进行一个确定的积分。
Y =(1 /(SQRT(2 * PI)))* E ^( - (1/2)* X^2)
这使我的区域中的两个x值之间的曲线下但然后我卡住了...也许我是基地的方式,这不是你会怎么做?
谢谢。
下面是在统计程序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;
}
在你的情况,你想亩= 0.0,标准差= 1.0,lower_tail = TRUE,log_p =假。
您缺少https://svn.r-project.org/R/trunk/src/nmath/dpq.h中定义的“R_D_Lval”和“R_D_Cval”的定义。在C#中它们如下:'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; }' – 2012-04-03 22:58:51
谢谢;我已经添加了这些答案。 – 2012-04-04 14:55:50
'-Base.ExpM1()'指的是什么? – 2013-01-14 23:29:12
对于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);
}
的新版本我很喜欢你的博客,谢谢! – Lukasz 2009-11-03 01:52:42
+1“#A&S公式7.1.26”。 Abramowitz和Stegun非常棒 - 每个从事数值工作的人都应该了解它。 – duffymo 2009-11-03 02:17:45
...而不是将它翻译成C#你自己,你可以点击“C#”链接。 – Contango 2011-05-27 13:16:13