Chức năng giá trị phân phối chuẩn Z bình thường trong C#

Question

Tôi đã xem bài đăng blog gần đây của Jeff Atwood trên Alternate Sorting Orders. Tôi đã cố gắng để chuyển đổi mã trong bài viết để C# nhưng tôi chạy vào một vấn đề. Không có hàm nào trong .NET mà tôi biết sẽ trả về giá trị z, được cho phần trăm diện tích theo đường chuẩn bình thường. Các giá trị được đề xuất để sử dụng cho thuật toán là 95% và 97,5% mà bạn có thể tra cứu trên bảng giá trị z trong bất kỳ cuốn sách thống kê nào.

Có ai biết cách triển khai chức năng như vậy cho tất cả các giá trị z hoặc ít nhất đến 6 độ lệch chuẩn so với giá trị trung bình. Một cách sẽ là mã hóa cứng các giá trị vào từ điển và sử dụng tra cứu nhưng phải có cách tính giá trị chính xác. Nỗ lực của tôi trong việc giải quyết điều này là lấy một tích phân xác định của hàm đường chuẩn bình thường.

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

này mang lại cho tôi diện tích dưới đường cong giữa hai x giá trị nhưng sau đó tôi bị mắc kẹt ... Có lẽ tôi là cách của cơ sở và đây không phải là cách bạn sẽ làm điều đó?

Cảm ơn.

Answer 1

Dưới đây là một số code cho bản phân phối bình thường được viết bằng Python, nhưng nó có thể dễ dàng được dịch sang C# bằng cách thêm một số dấu chấm câu. Nó chỉ là khoảng 15 dòng mã.

Answer 2

Tra cứu các triển khai của error function. Có một trong tất cả các Bí quyết số học cổ điển trong ... sách.

Answer 3

Đây là bản dịch C# của số lượng bình thường C code được sử dụng trong chương trình thống kê R.

/// <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); 
} 

Một số phương pháp helper:

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; 
} 

Trong trường hợp của bạn, bạn muốn mu = 0.0, sigma = 1.0, lower_tail = true, log_p = false.

Answer 4

Đối với một phiên bản mới hơn của 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); 
    }