我正在尝试编写自己的 Python 代码来计算一尾独立 t 检验和二尾独立 t 检验的 t 统计量和 p 值。我可以使用正态近似,但目前我尝试仅使用 t 分布。我未能成功地将 SciPy 统计库的结果与我的测试数据进行匹配。我可以用一双新的眼睛来看看我是否只是在某个地方犯了一个愚蠢的错误。
注意,这是从交叉验证交叉发布因为它已经存在了一段时间而没有任何回应,所以我认为获得一些软件开发人员的意见也没什么坏处。我试图了解我正在使用的算法是否存在错误,这应该会重现 SciPy 的结果。这是一个简单的算法,所以令人困惑的是为什么我无法定位错误。
My code:
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
Update:
在阅读了更多关于韦尔奇 t 检验的内容后,我发现我应该使用韦尔奇-萨特思韦特公式来计算自由度。我更新了上面的代码以反映这一点。
有了新的自由度,我得到了更接近的结果。我的两侧 p 值与 SciPy 版本相差约 0.008...但这仍然是一个太大的错误,所以我一定仍然做了一些不正确的事情(或者 SciPy 分布函数非常糟糕,但很难相信它们仅精确到小数点后两位)。
第二次更新:
在继续尝试的同时,我认为当自由度足够高(大约> 30)时,SciPy 的版本可能会自动计算 t 分布的正态近似。因此,我使用正态分布重新运行代码,计算结果实际上比使用 t 分布时更远离 SciPy。
奖金问题:)(更多统计理论相关;随意忽略)
此外,t 统计量为负。我只是想知道这对于单方面 t 检验意味着什么。这通常是否意味着我应该在负轴方向上进行测试?在我的测试数据中,人口1是没有接受过某种就业培训计划的对照组。人群2确实接受了治疗,测量的数据是治疗前后的工资差异。
所以我有理由认为人口 2 的平均值会更大。但从统计理论的角度来看,以这种方式炮制测试似乎并不正确。我怎么知道在不依赖于数据的主观知识的情况下要在负方向上进行检查(对于单方面的测试)?或者这只是那些虽然在哲学上不严谨但需要在实践中完成的常客主义事情之一?