65.9K
CodeProject 正在变化。 阅读更多。
Home

如何计算卡方 P 值

Jacob F. W.
starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.56/5 (19投票s)

2012年8月1日

公共领域

7分钟阅读

viewsIcon

128821

downloadIcon

3610

学习宇宙的奥秘

引言

您是否曾经想过那些卡方检验的 P 值对照表是如何计算出来的？

您是否曾经想过自己计算这些值？

您是否曾在网上求助计算这些值，却得到一个恼人的链接，上面只有对照表或一个可以帮您计算的表单小程序，但丝毫没有提及实际的公式！？！？！？！？ 

那么您来对地方了。请继续阅读，这篇文章就是为您准备的。

背景 

卡方分布（Chi-Squared Distribution）可能是当今统计学中最广泛使用的分布。它最为人熟知的是在皮尔逊卡方检验（Pearson's Chi-Squared Test）中的应用，该检验用于衡量拟合优度。

int Observed[256];
int Expected[256]; 
double CriticalValue = 0.0;
double XSqr;
 
for(int I = 0; I < 256; I++)
{
    XSqr = Observed[I] - Expected[I];
    CriticalValue += ((XSqr * XSqr) / Expected[I]);
}
 
/*
   Degrees of Freedom = 255 (The number of measure values - 1)
   Test Statistic = CriticalValue
*/

上面是皮尔逊卡方检验第一部分的简要示例。这演示了我们如何获得计算 P 值所需的两个值：自由度和临界值。 

什么是“P”值，它为何如此有价值？

对于那些对此不熟悉的朋友，P 值是指观测值与期望值之间的差异纯粹由偶然发生的概率。例如，一枚公平硬币的正面和反面出现的概率应该相等。

正面 = 50.0%
反面 = 50.0%

所以理论上，如果您抛掷一枚硬币 1000 次，正面和反面应该各出现 500 次。但实际上，这种情况很少发生。您更有可能出现正面朝上 507 次，反面 493 次，或者反之亦然。所以如果有些许差异，那是可以理解的。

但是，什么时候这种差异就太大了呢？一枚公平的硬币很可能出现正面朝上 507/1000 次。但如果出现 515 次呢？如果出现 600 次呢？  

这就是 P 值的作用。它允许您通过与期望值进行比较，客观地为观测值随机发生的概率赋值。

使用代码

上面我简要介绍了如何计算卡方检验的临界值和自由度。这仅用于演示，不包含在提供的代码中，因为您无疑需要根据自己的情况进行定制。如果您在这一方面有任何问题，网上有许多优秀的教程可以演示如何计算这些值。我强烈推荐其中的任何一个。

但这不是我在这里要讲的。我要向您展示如何计算实际的 P 值。所以，让我们来看一些代码。

double chisqr(int Dof, double Cv)
{
    if(Cv < 0 || Dof < 1)
    {
        return 0.0;
    }
    double K = ((double)Dof) * 0.5;
    double X = Cv * 0.5;
    if(Dof == 2)
    {
	return exp(-1.0 * X);
    }
 
    double PValue = igf(K, X);
    if(isnan(PValue) || isinf(PValue) || PValue <= 1e-8)
    {
        return 1e-14;
    } 

    PValue /= gamma(K);
    //PValue /= tgamma(K); 
	
    return (1.0 - PValue);
}

这是您将要使用的主要函数，位于 chisqr.c 中。它非常简单；您输入您的值，它就会输出 P 值。您会注意到，如果您只有 2 个自由度，那您就走运了。有一个快捷方法可以使用一次 exp() 调用来计算 P 值，这比 chisqr() 正常运行时要快得多。

您还需要另外两个函数来完成 P 值的计算。

static double igf(double S, double Z)
{
    if(Z < 0.0)
    {
	return 0.0;
    }
    double Sc = (1.0 / S);
    Sc *= pow(Z, S);
    Sc *= exp(-Z);
 
    double Sum = 1.0;
    double Nom = 1.0;
    double Denom = 1.0;
 
    for(int I = 0; I < 200; I++)
    {
	Nom *= Z;
	S++;
	Denom *= S;
	Sum += (Nom / Denom);
    }
 
    return Sum * Sc;
}

接下来是 不完全伽玛函数。它用于计算 P 值的第一部分。如果您想了解它与计算卡方值之间的关系，可以在 维基百科卡方分布页面的此部分找到说明。

顺便说一下，在 chisqr.c 中，我简要解释了为什么我的 for 循环会迭代 200 次。简而言之，这是我通过反复试验得出的结果。如果您在这方面有经验，并且对此有更多了解，请在下方的评论中告诉我。

最后但同样重要…… 

double gamma(double N)
{
    const long double SQRT2PI = 2.5066282746310005024157652848110452530069867406099383;
 
    long double Z = (long double)N;
    long double Sc = powl((Z + A), (Z + 0.5));
    Sc *= expl(-1.0 * (Z + A));
    Sc /= Z;
 
    long double F = 1.0;
    long double Ck;
    long double Sum = SQRT2PI;
 

    for(int K = 1; K < A; K++)
    {
        Z++;
	Ck = powl(A - K, K - 0.5);
	Ck *= expl(A - K);
	Ck /= F;
 
	Sum += (Ck / Z);
 
	F *= (-1.0 * K);
    }
 
    return (double)(Sum * Sc);
}

向 伽玛函数问好。对于不熟悉它的人来说，伽玛函数允许您计算小数的阶乘（例如 5.5!）。

快速附注。如果您以前从未使用过伽玛函数，有一件事您应该知道。它实际上会计算您输入的数字减 1 的阶乘。Gamma(X) = (X - 1)! 所以如果您输入 Gamma(5)，您不会得到 120，而是 24。“为什么？”您会问。我不知道：数学家就是这么奇怪。我认为这是他们惩罚社会多年来称他们为书呆子的方式，但这只是我的理论。

取决于您使用的语言，您的标准数学库可能已经内置了伽玛函数。例如，在 C 的 math.h 中，您可以通过调用 tgamma() 来计算伽玛函数。我建议您检查一下您是否有此函数，因为标准库中的实现通常更有效、更快、更准确。

这里您看到的伽玛函数是 Spouge 逼近法的实现。我喜欢它的原因在于，只要您有足够大的数据类型来存储所需数据，它就能让您以任意精度计算伽玛函数。您可以使用任意精度库（事实上，我认为它已经实现了）来实现这一点，并计算您所需的精度。

但是，它并不是最快的伽玛函数版本。有几种逼近法可以给出对于卡方检验而言足够准确的值。

double approx_gamma(double Z)
{
    const double RECIP_E = 0.36787944117144232159552377016147;  // RECIP_E = (E^-1) = (1.0 / E)
    const double TWOPI = 6.283185307179586476925286766559;  // TWOPI = 2.0 * PI

    double D = 1.0 / (10.0 * Z);
    D = 1.0 / ((12 * Z) - D);
    D = (D + Z) * RECIP_E;
    D = pow(D, Z);
    D *= sqrt(TWOPI / Z);
 
    return D;
}

这是一个简单的 斯特林逼近法。这将为您提供大约 6-7 位的精度。我将使用此方法得到的 P 值与我使用其他实现和 math.h 中的值进行了比较。

// Degrees_of_Freedom = 255
// Critical_Value = 290.285192
// Output = P_Value
// chisqr(Degrees_of_Freedom, Critical_Value) = P_Value

chisqr(255, 290.285192) = 0.06364235  // Using tgamma() in math.h
chisqr(255, 290.285192) = 0.06364235  // Using gamma()
chisqr(255, 290.285192) = 0.06364235  // Using approx_gamma()

所以，是的，我尝试了数学库中的伽玛函数、我实现的 Spouge 逼近法以及斯特林逼近法，并将它们打印到屏幕上：没有区别。一些低位可能不同，但这已经比您在网上找到的任何表格或计算器都精确得多。所以，除非您追求绝对精确，否则我建议您坚持使用 approx_gamma()。

结论

天哪，代码量真不少。感谢大家一直坚持到现在。我发现这些算法非常重要，我很乐意与大家分享它们。如果您有任何问题，请随时提出，我会尽力回答。如果您想了解关于具体公式本身的任何信息，我建议您查阅我提供的链接。  如果您有任何评论、批评、建议或改进意见，请随时发表评论。

最后一点。在过去几年里，我偶尔会需要这些函数，并费力地寻找它们，却常常发现它们被置于非常严格的许可之下。这本身并不是一件坏事，但考虑到这些函数如此重要，并且在数学和测试设备中如此受欢迎，我惊讶于竟然很少有实现是免费提供的。

所以，为了绝对清楚。

这些代码已放入公共领域，任何人都可以出于任何原因、任何目的免费使用，无需获得许可，等等。 

我承认这些不是这些函数的最佳实现，但它们可以工作，而且总比什么都没有好。 

谢谢，

Jacob W.

更新 

2012/8/2：我在进行一些测试时，意识到我忘记考虑 PValue 为 nan 或 inf 的情况。我已经更新了示例，并附带了更正后的版本，以及 ZIP 文件。如果这让任何人感到困惑，我深表歉意。

2012/11/6：上传了 chisqr_v2.zip。ChiSqr() 已经进行了重大升级。我修改了代码，使其现在使用自然对数来帮助计算 ChiSqr P 值。这极大地增加了其范围和准确性。不过，速度变慢了。很遗憾，我还没有找到改进这一方面的方法。尽管如此，我还是强烈建议使用这个新版本而不是旧版本，因为我认为它非常值得。

 

© . All rights reserved.