8

我正在尝试编写自己的 Python 代码来计算一个和两个有尾独立 t 检验的 t 统计量和 p 值。我可以使用正态近似值，但目前我试图只使用 t 分布。我未能将 SciPy 统计库的结果与我的测试数据相匹配。我可以用一双新的眼睛来看看我是否只是在某个地方犯了一个愚蠢的错误。

请注意，这是从 Cross-Validated 交叉发布的，因为它在那里已经有一段时间没有回复了，所以我认为获得一些软件开发人员的意见也无妨。我试图了解我正在使用的算法是否存在错误，它应该重现 SciPy 的结果。这是一个简单的算法，所以很奇怪为什么我找不到错误。

我的代码：

import numpy as np
import scipy.stats as st

def compute_t_stat(pop1,pop2):

    num1 = pop1.shape[0]; num2 = pop2.shape[0];

    # The formula for t-stat when population variances differ.
    t_stat = (np.mean(pop1) - np.mean(pop2))/np.sqrt( np.var(pop1)/num1 + np.var(pop2)/num2 )

    # ADDED: The Welch-Satterthwaite degrees of freedom.
    df = ((np.var(pop1)/num1 + np.var(pop2)/num2)**(2.0))/(   (np.var(pop1)/num1)**(2.0)/(num1-1) +  (np.var(pop2)/num2)**(2.0)/(num2-1) ) 

    # Am I computing this wrong?
    # It should just come from the CDF like this, right?
    # The extra parameter is the degrees of freedom.

    one_tailed_p_value = 1.0 - st.t.cdf(t_stat,df)
    two_tailed_p_value = 1.0 - ( st.t.cdf(np.abs(t_stat),df) - st.t.cdf(-np.abs(t_stat),df) )    


    # Computing with SciPy's built-ins
    # My results don't match theirs.
    t_ind, p_ind = st.ttest_ind(pop1, pop2)

    return t_stat, one_tailed_p_value, two_tailed_p_value, t_ind, p_ind

更新：

在阅读了更多关于 Welch t 检验的内容后，我看到我应该使用 Welch-Satterthwaite 公式来计算自由度。我更新了上面的代码以反映这一点。

有了新的自由度，我得到了更接近的结果。我的双边 p 值与 SciPy 版本相差约 0.008 ......但这仍然是一个太大的错误，所以我仍然必须做一些不正确的事情（或者 SciPy 分布函数非常糟糕，但很难相信它们只精确到小数点后 2 位）。

第二次更新：

在继续尝试的同时，我认为当自由度足够高（大约 > 30）时，SciPy 的版本可能会自动计算 t 分布的正态近似值。所以我改用正态分布重新运行我的代码，计算结果实际上比我使用 t 分布时更远离 SciPy。

奖金问题: )（更多统计理论相关；随意忽略）

此外，t 统计量为负。我只是想知道这对单边 t 检验意味着什么。这是否通常意味着我应该在负轴方向进行测试？在我的测试数据中，人群 1 是没有接受过某种就业培训计划的对照组。人口 2 确实收到了，测量的数据是治疗前后的工资差异。

所以我有理由认为人口 2 的平均值会更大。但从统计理论的角度来看，以这种方式编造一个测试似乎并不正确。我怎么知道在不依赖对数据的主观知识的情况下检查（对于单方面的测试）？或者这只是那些虽然在哲学上不严谨但需要在实践中完成的常客事情之一？

4

3 回答 3

9

通过使用 SciPy 内置函数source()，我可以看到该函数源代码的打印输出ttest_ind()。根据源代码，内置 SciPy 正在执行 t 检验，假设两个样本的方差相等。它没有使用 Welch-Satterthwaite 自由度。SciPy 假设方差相等，但没有说明这个假设。

我只想指出，至关重要的是，这就是为什么您不应该只信任库函数的原因。在我的情况下，我实际上确实需要对不等方差的总体进行 t 检验，并且自由度调整对于我将运行的一些较小的数据集可能很重要。

正如我在一些评论中提到的那样，对于 30 到 400 之间的样本量，我的代码和 SciPy 之间的差异约为 0.008，然后对于更大的样本量会慢慢变为零。这是等方差 t 统计量分母中额外 (1/n1 + 1/n2) 项的影响。准确性方面，这非常重要，尤其是对于小样本量。它肯定向我证实了我需要编写自己的函数。（可能还有其他更好的 Python 库，但至少应该知道这一点。坦率地说，令人惊讶的是，这并不是 SciPy 文档中的任何位置ttest_ind()）。

于 2012-04-12T21:07:40.920 回答
2

您不是在计算样本方差，而是在使用总体方差。样本方差除以n-1，而不是n。np.var有一个可选参数调用ddof，原因与此类似。

这应该会给您预期的结果：

import numpy as np
import scipy.stats as st

def compute_t_stat(pop1,pop2):

    num1 = pop1.shape[0]
    num2 = pop2.shape[0];
    var1 = np.var(pop1, ddof=1)
    var2 = np.var(pop2, ddof=1)

    # The formula for t-stat when population variances differ.
    t_stat = (np.mean(pop1) - np.mean(pop2)) / np.sqrt(var1/num1 + var2/num2)

    # ADDED: The Welch-Satterthwaite degrees of freedom.
    df = ((var1/num1 + var2/num2)**(2.0))/((var1/num1)**(2.0)/(num1-1) + (var2/num2)**(2.0)/(num2-1)) 

    # Am I computing this wrong?
    # It should just come from the CDF like this, right?
    # The extra parameter is the degrees of freedom.

    one_tailed_p_value = 1.0 - st.t.cdf(t_stat,df)
    two_tailed_p_value = 1.0 - ( st.t.cdf(np.abs(t_stat),df) - st.t.cdf(-np.abs(t_stat),df) )    


    # Computing with SciPy's built-ins
    # My results don't match theirs.
    t_ind, p_ind = st.ttest_ind(pop1, pop2)

    return t_stat, one_tailed_p_value, two_tailed_p_value, t_ind, p_ind

PS：SciPy 是开源的，主要使用 Python 实现。您可以自己检查源代码ttest_ind并找出错误。

对于奖励方面：您不会通过查看 t 值来决定单尾检验的一方。你事先用你的假设来决定它。如果您的零假设是均值相等，而您的备择假设是第二个均值较大，那么您的尾巴应该在左侧（负）侧。因为您的 t 值足够小（负）值将表明备择假设更有可能为真，而不是零假设。

于 2012-04-06T04:33:10.787 回答
0

看起来您忘记了 df 的分子 **2。Welch-Satterthwaite 自由度。

df = (np.var(pop1)/num1 + np.var(pop2)/num2)/(   (np.var(pop1)/num1)**(2.0)/(num1-1) +  (np.var(pop2)/num2)**(2.0)/(num2-1) )

应该：

df = (np.var(pop1)/num1 + np.var(pop2)/num2)**2/(   (np.var(pop1)/num1)**(2.0)/(num1-1) +  (np.var(pop2)/num2)**(2.0)/(num2-1) )
于 2012-04-06T03:43:35.123 回答