вычисление квадрата ци с использованием scipy - PullRequest
1 голос
/ 10 апреля 2020

Я попробовал тест хи-квадрат, используя 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

Не знаю, что не так.

1 Ответ

1 голос
/ 10 апреля 2020

Нет ничего плохого! Разница, которую вы видите, состоит в том, что scipy.stats.chi2_contingency применяет «коррекцию непрерывности», когда входной массив равен 2x2. Вы можете отключить это исправление, передав аргумент correction=False. При этом выходные данные соответствуют вашим расчетам:

In [12]: vals = np.array([[70, 20], [50, 60]])

In [13]: x2, p, dof, expected = scipy.stats.chi2_contingency(vals, correction=False)

In [14]: x2
Out[14]: 21.54882154882155

In [15]: p
Out[15]: 3.449345750127958e-06
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...