Я попробовал тест хи-квадрат, используя scipy, как показано ниже
import numpy as np
import scipy
vals = np.array([[70, 20], [50, 60]])
x2, p, dof, expected = scipy.stats.chi2_contingency(vals)
print('x2 = {:.5f}'.format(x2))
print('p-value = {}'.format(p))
print(expected)
a = scipy.stats.chisquare(f_obs= vals, # Array of observed counts
f_exp= expected)
Я получил
- x2 = 20.22306
- p-значение = 6.8917007718498866e-06
- [[54. 36.] [66. 44.]]
Однако этот результат отличается от моей реализации.
def Chi2_test(vals, k=1):
r, c = vals.shape
a_sum = vals.sum(axis=0)
b_sum = vals.sum(axis=1)
S = vals.sum()
Pa= a_sum / S
Pb = b_sum / S
Pa = np.tile(Pa, c).reshape(r, c)
Pb = np.repeat(Pb, c).reshape(r, c)
Pab = Pa * Pb
E = Pab * S
x2 = np.sum(((vals - E) ** 2) / E)
# chi square -> p value
# Gamma function
def Gamma(x):
if x == 1:
return 1
elif x == 0.5:
return np.sqrt(np.pi)
else:
return (x - 1) * Gamma(x - 1)
# chi square
def Chi2(x, k):
return (x ** (k / 2 - 1)) * (np.exp(- x / 2)) / ((2 ** (k / 2)) * Gamma(k / 2))
p_value = integrate.quad(lambda x: Chi2(x, k=k), x2, np.inf)[0]
return x2, p_value
vals = np.array([[70, 20],
[50, 60]])
x2, p_value = Chi2_test(vals)
print('x2 :', x2)
print('p-value :', p_value)
- x2: 21.548821548821547
- p-значение: 3.449345362777984e -06
Не знаю, что не так.