Критерий Пирсона: проверка гипотез и анализ данных на Python

pip install numpy scipy pandas matplotlib seaborn

import numpy as np
import scipy.stats as stats
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Для красивого отображения графиков
plt.style.use('ggplot')
sns.set(style="whitegrid")

# Создаем данные: опыт (строки) × язык программирования (столбцы)
observed = np.array([
[30, 14, 6], # Junior (Python, Java, C++)
[45, 30, 25], # Middle (Python, Java, C++)
[25, 36, 49] # Senior (Python, Java, C++)
])

# Создаем датафрейм для лучшей визуализации
df_observed = pd.DataFrame(
observed, 
index=['Junior', 'Middle', 'Senior'],
columns=['Python', 'Java', 'C++']
)
print("Наблюдаемые частоты:")
print(df_observed)

# Суммы по строкам и столбцам
row_totals = observed.sum(axis=1)
col_totals = observed.sum(axis=0)
total = observed.sum()

# Расчет ожидаемых частот
expected = np.outer(row_totals, col_totals) / total

# Создаем датафрейм ожидаемых частот
df_expected = pd.DataFrame(
expected, 
index=['Junior', 'Middle', 'Senior'],
columns=['Python', 'Java', 'C++']
)
print("\nОжидаемые частоты:")
print(df_expected)

# Ручной расчет статистики хи-квадрат
chi2_stat = np.sum((observed – expected) ** 2 / expected)

# Степени свободы = (строки – 1) * (столбцы – 1)
dof = (observed.shape[0] – 1) * (observed.shape[1] – 1)

# p-значение
p_value = 1 – stats.chi2.cdf(chi2_stat, dof)

print(f"\nСтатистика χ²: {chi2_stat:.4f}")
print(f"Степени свободы: {dof}")
print(f"p-значение: {p_value:.6f}")

# Использование встроенной функции SciPy
chi2_result = stats.chi2_contingency(observed)
print("\nРезультаты с использованием stats.chi2_contingency:")
print(f"Статистика χ²: {chi2_result[0]:.4f}")
print(f"p-значение: {chi2_result[1]:.6f}")
print(f"Степени свободы: {chi2_result[2]}")

# Сравнение результатов
print("\nРазница между ручным расчетом и функцией SciPy:")
print(f"Разница в статистике χ²: {abs(chi2_stat – chi2_result[0]):.10f}")
print(f"Разница в p-значении: {abs(p_value – chi2_result[1]):.10f}")

# Создаем данные: образование × ОС
education_os = np.array([
[150, 100, 50], # Среднее (Windows, macOS, Linux)
[120, 130, 50], # Бакалавр (Windows, macOS, Linux)
[80, 150, 120] # Магистр/PhD (Windows, macOS, Linux)
])

# Проведем тест хи-квадрат
chi2, p, dof, expected = stats.chi2_contingency(education_os)

print(f"Статистика χ²: {chi2:.4f}")
print(f"p-значение: {p:.6f}")
print(f"Степени свободы: {dof}")

# Интерпретация результата
alpha = 0.05
print("\nИнтерпретация:")
if p < alpha:
print(f"p-значение ({p:.6f}) < {alpha}: отклоняем нулевую гипотезу")
print("Существует зависимость между уровнем образования и предпочтениями ОС")
else:
print(f"p-значение ({p:.6f}) ≥ {alpha}: не отклоняем нулевую гипотезу")
print("Нет достаточных доказательств зависимости между образованием и предпочтениями ОС")

# Наблюдаемые частоты оценок
observed_grades = np.array([10, 25, 30, 25, 10]) # Оценки от 1 до 5

# Ожидаемое равномерное распределение
expected_grades = np.ones(5) * np.sum(observed_grades) / 5

# Проведем тест хи-квадрат
chi2_stat = np.sum((observed_grades – expected_grades) ** 2 / expected_grades)
p_value = 1 – stats.chi2.cdf(chi2_stat, len(observed_grades) – 1)

print(f"Наблюдаемые частоты: {observed_grades}")
print(f"Ожидаемые частоты: {expected_grades}")
print(f"Статистика χ²: {chi2_stat:.4f}")
print(f"p-значение: {p_value:.6f}")

# Интерпретация
if p_value < 0.05:
print("Распределение оценок не соответствует равномерному")
else:
print("Распределение оценок соответствует равномерному")

# Создаем данные
observed_data = np.array([
[30, 15, 5],
[25, 25, 10],
[20, 30, 40]
])

# Рассчитываем ожидаемые частоты
chi2, p, dof, expected = stats.chi2_contingency(observed_data)

# Создаем датафреймы для визуализации
df_observed = pd.DataFrame(observed_data, 
index=['Группа A', 'Группа B', 'Группа C'],
columns=['Категория 1', 'Категория 2', 'Категория 3'])

df_expected = pd.DataFrame(expected, 
index=['Группа A', 'Группа B', 'Группа C'],
columns=['Категория 1', 'Категория 2', 'Категория 3'])

# Настраиваем график
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Визуализация наблюдаемых частот
sns.heatmap(df_observed, annot=True, cmap="YlGnBu", fmt=".1f", ax=axes[0])
axes[0].set_title('Наблюдаемые частоты')

# Визуализация ожидаемых частот
sns.heatmap(df_expected, annot=True, cmap="YlOrRd", fmt=".1f", ax=axes[1])
axes[1].set_title('Ожидаемые частоты')

plt.suptitle(f'Анализ частот (p-значение: {p:.6f})', fontsize=16)
plt.tight_layout()
plt.show()

# Расчет остатков (разница между наблюдаемыми и ожидаемыми частотами)
residuals = (observed_data – expected) / np.sqrt(expected)
df_residuals = pd.DataFrame(residuals,
index=['Группа A', 'Группа B', 'Группа C'],
columns=['Категория 1', 'Категория 2', 'Категория 3'])

# Визуализация остатков
plt.figure(figsize=(10, 8))
sns.heatmap(df_residuals, annot=True, cmap="coolwarm", center=0, fmt=".2f")
plt.title('Стандартизированные остатки')
plt.tight_layout()
plt.show()

# Создаем диапазон значений для распределения хи-квадрат
x = np.linspace(0, 30, 1000)
y = stats.chi2.pdf(x, dof) # PDF для распределения хи-квадрат с df степенями свободы

# Рисуем распределение и область p-значения
plt.figure(figsize=(10, 6))
plt.plot(x, y, 'b-', lw=2, label=f'χ² распределение (df={dof})')

# Заштрихованная область для критического значения
critical_value = stats.chi2.ppf(0.95, dof)
plt.fill_between(x, y, where=(x >= critical_value), color='red', alpha=0.3, 
label=f'Критическая область (α=0.05)')

# Вертикальная линия для наблюдаемого значения χ²
plt.axvline(chi2, color='green', linestyle='--', lw=2, 
label=f'Наблюдаемое χ²={chi2:.4f}')

# Заштрихованная область для p-значения
plt.fill_between(x, y, where=(x >= chi2), color='green', alpha=0.2, 
label=f'p-значение={p:.4f}')

plt.xlabel('χ² значение')
plt.ylabel('Плотность вероятности')
plt.title('Распределение χ² и p-значение')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Мозаичный график для визуализации взаимосвязи
from statsmodels.graphics.mosaicplot import mosaic

plt.figure(figsize=(10, 8))
props = {}
labelizer = lambda k: {('Группа A', 'Категория 1'): 'A1', 
('Группа A', 'Категория 2'): 'A2',
('Группа A', 'Категория 3'): 'A3',
('Группа B', 'Категория 1'): 'B1', 
('Группа B', 'Категория 2'): 'B2',
('Группа B', 'Категория 3'): 'B3',
('Группа C', 'Категория 1'): 'C1', 
('Группа C', 'Категория 2'): 'C2',
('Группа C', 'Категория 3'): 'C3'}.get(k, '')

contingency_data = {(i, j): observed_data[i, j] 
for i in range(3) 
for j in range(3)}

mosaic(contingency_data, properties=props, labelizer=labelizer, 
title=f'Мозаичный график (p={p:.4f})', axes_label=True)

plt.tight_layout()
plt.show()

# Данные о конверсии для трех типов маркетинговых кампаний
marketing_data = np.array([
[120, 880], # Кампания A: [конверсия, нет конверсии]
[150, 850], # Кампания B: [конверсия, нет конверсии]
[200, 800] # Кампания C: [конверсия, нет конверсии]
])

# Проведение теста хи-квадрат
chi2, p, dof, expected = stats.chi2_contingency(marketing_data)

print(f"Кейс 1: Эффективность маркетинговых кампаний")
print(f"Статистика χ²: {chi2:.4f}")
print(f"p-значение: {p:.6f}")

# Интерпретация результата
if p < 0.05:
print("Существуют статистически значимые различия в эффективности кампаний")

# Анализ остатков для определения конкретных различий
residuals = (marketing_data – expected) / np.sqrt(expected)
for i, campaign in enumerate(['A', 'B', 'C']):
for j, outcome in enumerate(['Конверсия', 'Нет конверсии']):
print(f"Кампания {campaign}, {outcome}: остаток = {residuals[i, j]:.2f}")

if abs(residuals[i, j]) > 1.96:
direction = "выше" if residuals[i, j] > 0 else "ниже"
print(f" Значимо {direction} ожидаемого (p < 0.05)")
else:
print("Нет статистически значимых различий в эффективности кампаний")

# Данные о предпочтениях продуктов по возрастным группам
age_product_data = np.array([
[150, 100, 50], # 18-24: [продукт A, продукт B, продукт C]
[120, 150, 130], # 25-34: [продукт A, продукт B, продукт C]
[100, 120, 180], # 35-44: [продукт A, продукт B, продукт C]
[80, 100, 220] # 45+: [продукт A, продукт B, продукт C]
])

# Проведение теста
chi2, p, dof, expected = stats.chi2_contingency(age_product_data)

print(f"\nКейс 2: Предпочтения продуктов в разных возрастных группах")
print(f"Статистика χ²: {chi2:.4f}")
print(f"p-значение: {p:.10f}") # Очень малое p-значение
print(f"Степени свободы: {dof}")

# Интерпретация p-значения
if p < 0.001:
significance = "очень сильная"
elif p < 0.01:
significance = "сильная"
elif p < 0.05:
significance = "умеренная"
else:
significance = "отсутствует"

print(f"Статистическая значимость: {significance}")

# Расчет коэффициента V Крамера для оценки силы связи
def cramers_v(confusion_matrix):
chi2 = stats.chi2_contingency(confusion_matrix)[0]
n = confusion_matrix.sum()
phi2 = chi2/n
r, k = confusion_matrix.shape
phi2corr = max(0, phi2 – ((k-1)*(r-1))/(n-1))
rcorr = r – ((r-1)**2)/(n-1)
kcorr = k – ((k-1)**2)/(n-1)
return np.sqrt(phi2corr / min((kcorr-1), (rcorr-1)))

# Рассчитываем V Крамера для обоих примеров
v_marketing = cramers_v(marketing_data)
v_age_product = cramers_v(age_product_data)

print("\nСила связи (V Крамера):")
print(f"Кейс 1 (маркетинг): {v_marketing:.4f}")
print(f"Кейс 2 (возраст-продукт): {v_age_product:.4f}")

# Интерпретация силы связи
def interpret_cramers_v(v):
if v < 0.1:
return "незначительная связь"
elif v < 0.3:
return "слабая связь"
elif v < 0.5:
return "умеренная связь"
elif v < 0.8:
return "сильная связь"
else:
return "очень сильная связь"

print(f"Интерпретация кейс 1: {interpret_cramers_v(v_marketing)}")
print(f"Интерпретация кейс 2: {interpret_cramers_v(v_age_product)}")