一种有效自由度的python实现与双尾t检验测试

用户11172986

发布于 2024-06-20 16:38:10

36900

代码可运行

文章被收录于专栏：气python风雨气python风雨

运行总次数：0

代码可运行

版本：python3.7 数据：随机生成的正态分布数组

�eff≈1�+2�∑�=1��−��(�)��(�

这里 N 是样本大小，ρXX （j）和 ρYY （j）分别是两个采样时间序列 X 和 Y 在时间滞后 j 处的自相关。

NAO implicated as a predictor of the surface air temperature multidecadal variability over East Asia

Tiejun Xie, Jianping Li

这个公式是用于估计两个时间序列X和Y之间的有效样本量（Effective Sample Size）或者是有效自由度，并且有以下各项含义：公式中的每一项都有特定的含义：

第一项

表示基础样本量的倒数，反映了样本大小对估计精度的贡献。第二项

�∑�=1��−��(�)

是关于自相关系数的求和部分，其中

是一种加权系数，用于调整每个时间滞后的贡献，\rho_{xx}(j)\rho_{yy}(j) 表示X和Y在时间滞后j处的自相关系数的乘积。通过估计有效样本量，我们可以更准确地考虑到自相关性的影响，并进行统计推断和假设检验等分析。

函数实现

np.correlate

def compute_equation(N, X, Y):
    Neff = 1 / N  # 初始化Neff的值为1 / N
    ρ_xx = np.correlate(X, X, mode='full')[N-1:]  # 计算X的自相关系数
    ρ_yy = np.correlate(Y, Y, mode='full')[N-1:]  # 计算Y的自相关系数
    
    for j in range(1, N+1):
        Neff += 2/N  *  ((N - j) / N) * ρ_xx[j-1] * ρ_yy[j-1]
        
    return 1 / Neff


# 调用函数并打印结果
N = 100
xx = np.random.normal(0, 1, 100)
yy = np.random.normal(0, 1, 100)

result = compute_equation(N, xx, yy)
print(result)

0.006419097817664573

sm.tsa.acf

import numpy as np
import statsmodels.api as sm

def compute_equation2(N, X, Y):
    Neff = 1 / N  # 初始化Neff的值为1 / N
    ρ_xx = sm.tsa.acf(X, nlags=N-1)  # 计算X的自相关系数
    ρ_yy = sm.tsa.acf(Y, nlags=N-1)  # 计算Y的自相关系数
    
    for j in range(1, N+1):
        Neff += 2/N  *  ((N - j) / N) * ρ_xx[j-1] * ρ_yy[j-1]
        
    return 1 / Neff


# 调用函数并打印结果
N = 100
xx = np.random.normal(0, 1, 100)
yy = np.random.normal(0, 1, 100)

result = compute_equation2(N, xx, yy)
print(result)

33.852749934880265

数据测试

import numpy as np
# 输入两组样本数据（假设为sample1和sample2）
sample1 = np.random.normal(0, 1, 100)
sample2 = np.random.normal(0, 1, 100)

import scipy.stats as stats
import matplotlib.pyplot as plt



# 执行双尾t检验
statistic, p_value = stats.ttest_ind(sample1, sample2)

# 输出结果
print("T统计量：", statistic)
print("P值：", p_value)

# 计算有效自由度数
N = len(sample1)
df = 1 / compute_equation(N, sample1, sample2)

# 计算t分布的临界值
alpha = 0.90
t_critical = stats.t.ppf(1 - alpha / 2, df)
print('tc',t_critical)

T统计量：-1.7797481108940942
P值：0.076650975597045
tc 0.1258061608005306

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# # 输入两组样本数据（假设为sample1和sample2）
# sample1 = np.random.normal(0, 1, 100)
# sample2 = np.random.normal(0, 1, 100)

# 执行双尾t检验
statistic, p_value = stats.ttest_ind(sample1, sample2)

# 输出结果
print("T统计量：", statistic)
print("P值：", p_value)

# 计算有效自由度数
N2 = len(sample1)
df2 = 1 / compute_equation2(N, sample1, sample2)

# 计算t分布的临界值
alpha = 0.90
t_critical2 = stats.t.ppf(1 - alpha / 2, df2)
print('tc',t_critical2)

T统计量：-1.7797481108940942
P值：0.076650975597045
tc 2.591595063118441

绘图

numpy计算结果

# 设置图形大小
plt.figure(figsize=(10, 6))  # 设置宽度为10，高度为6

# 绘制折线图
x = np.arange(len(sample1))
plt.plot(x, sample1, label='Sample 1')
plt.plot(x, sample2, label='Sample 2')
plt.xlabel('Index')
plt.ylabel('Value')
plt.title('Line Chart of Samples')
plt.legend(loc='upper left')

# 绘制虚线区域
n1 = len(sample1)
n2 = len(sample2)
fill_x = np.array([min(x), max(x)])
fill_y = t_critical * np.sqrt((np.var(sample1) / n1) + (np.var(sample2) / n2))
plt.fill_between(fill_x, -fill_y, fill_y, color='gray', alpha=0.3)

# 显示图形
plt.show()

statsmodels计算结果

# 设置图形大小
plt.figure(figsize=(10, 6))  # 设置宽度为10，高度为6

# 绘制折线图
x = np.arange(len(sample1))
plt.plot(x, sample1, label='Sample 1')
plt.plot(x, sample2, label='Sample 2')
plt.xlabel('Index')
plt.ylabel('Value')
plt.title('Line Chart of Samples')
plt.legend(loc='upper left')

# 绘制虚线区域
n1 = len(sample1)
n2 = len(sample2)
fill_x = np.array([min(x), max(x)])
fill_y = t_critical2 * np.sqrt((np.var(sample1) / n1) + (np.var(sample2) / n2))
plt.fill_between(fill_x, -fill_y, fill_y, color='gray', alpha=0.3)

# 显示图形
plt.show()

小结

画出来和文献不一样，可能数据处理不当或者数据差异导致。
随机数组是随机，f需要多试几次，也可自己搓一个符合的数组。
同一数组计算对比来说，sm.tsa.acf函数计算结果置信区间更大。

封面生成：ai

/imagine prompt: color photo of a mathematician surrounded by equations and calculations, busy calculating degrees of freedom in a complex system —c 10 —ar 2:3

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2023-11-30，如有侵权请联系 cloudcommunity@tencent.com 删除

数组