diff --git a/idz3/ИДЗ 3_2 Артём.ipynb b/idz3/ИДЗ 3_2 Артём.ipynb index 9821edc..b6b3ddf 100644 --- a/idz3/ИДЗ 3_2 Артём.ipynb +++ b/idz3/ИДЗ 3_2 Артём.ipynb @@ -25,7 +25,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 4, "id": "57a523dd", "metadata": {}, "outputs": [], @@ -58,7 +58,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 5, "id": "db7e1a67", "metadata": {}, "outputs": [ @@ -95,7 +95,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 6, "id": "261ad18a", "metadata": {}, "outputs": [ @@ -159,7 +159,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 7, "id": "09541433", "metadata": {}, "outputs": [ @@ -214,7 +214,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 8, "id": "ead66cb6", "metadata": {}, "outputs": [ @@ -256,7 +256,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 9, "id": "a24ea7eb", "metadata": {}, "outputs": [ @@ -286,7 +286,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 10, "id": "e8490052", "metadata": {}, "outputs": [ @@ -316,7 +316,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 11, "id": "cc21a5b6", "metadata": {}, "outputs": [ @@ -347,7 +347,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 12, "id": "118d475e", "metadata": {}, "outputs": [ @@ -379,7 +379,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 13, "id": "08ea631c", "metadata": {}, "outputs": [ @@ -438,7 +438,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 14, "id": "7fa556a6", "metadata": {}, "outputs": [ @@ -470,7 +470,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 15, "id": "91b57523", "metadata": {}, "outputs": [ @@ -494,20 +494,24 @@ "metadata": {}, "source": [ "### 3. Смещение оценок\n", - "Для показательного распределения ОМП и ОММ совпадают. Найдём смещение:\n", + "Для показательного распределения ОМП и ОММ совпадают. Найдём смещение для $\\hat{\\lambda} = \\frac{1}{\\bar{X}}$:\n", "$$\n", "\\text{Смещение}(\\hat{\\lambda}) = E[\\hat{\\lambda}] - \\lambda\n", "$$\n", - "Для показательного распределения:\n", "$$\n", - "E[\\hat{\\lambda}_{\\text{ОМП}}] = E\\left[\\frac{1}{\\bar{X}}\\right] \\neq \\frac{1}{E[\\bar{X}]} = \\lambda\n", + "E[\\hat{\\lambda}_{\\text{ОМП}}] = E\\left[\\frac{1}{\\bar{X}}\\right] = \\frac{n\\lambda}{n-1} \\neq \\frac{1}{E[\\bar{X}]} = \\lambda \\space(n>1)\n", "$$\n", - "Оценка $\\hat{\\lambda}_{\\text{ОМП}}$ является смещённой, но асимптотически несмещённой." + "Смещение:\n", + "$$\n", + "\\text{Bias}(\\hat{\\lambda}) = \\frac{n\\lambda}{n-1} - \\lambda = \\frac{\\lambda}{n-1}.\n", + "$$\n", + "\n", + "Оценка $\\hat{\\lambda}_{\\text{ОМП}}$ является смещённой. Смещение положительно и уменьшается с ростом $n$. При $n→\\inf$, $Bias(\\hat{\\lambda})→0$, поэтому оценка *асимптотически несмещенная*." ] }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 16, "id": "1bb3bbda", "metadata": {}, "outputs": [ @@ -536,20 +540,54 @@ "metadata": {}, "source": [ "## Пункт d) Асимптотический доверительный интервал\n", - "Для построения асимптотического доверительного интервала используем тот факт, что ОМП асимптотически нормальна с дисперсией:\n", + "\n", + "### Оценка максимального правдоподобия (ОМП):\n", + "Для экспоненциального распределения ОМП параметра **λ** вычисляется как:\n", "$$\n", - "\\text{Var}(\\hat{\\lambda}) = \\frac{\\lambda^2}{n}\n", + "\\hat{\\lambda} = \\frac{1}{\\bar{X}},\n", "$$\n", - "Доверительный интервал уровня значимости $\\alpha_2$ имеет вид:\n", + "где $\\bar{X}$ — выборочное среднее.\n", + "\n", + "\n", + "### Асимптотическая нормальность ОМП:\n", + "При $n \\to \\infty$ распределение $\\hat{\\lambda}$ стремится к нормальному:\n", + "- **Математическое ожидание**:\n", + " $$\n", + " \\mathbb{E}[\\hat{\\lambda}] = \\lambda\n", + " $$\n", + "- **Дисперсия**:\n", + " $$\n", + " \\text{Var}(\\hat{\\lambda}) = \\frac{\\lambda^2}{n},\n", + " $$\n", + "\n", + "\n", + "### Стандартная ошибка оценки:\n", + "Оценивается через ОМП:\n", "$$\n", - "\\hat{\\lambda} \\pm z_{1-\\alpha_2/2} \\cdot \\frac{\\hat{\\lambda}}{\\sqrt{n}}\n", + "\\text{SE}(\\hat{\\lambda}) = \\sqrt{\\frac{\\hat{\\lambda}^2}{n}} = \\frac{\\hat{\\lambda}}{\\sqrt{n}}.\n", "$$\n", - "где $z_{1-\\alpha_2/2}$ — квантиль стандартного нормального распределения." + "\n", + "\n", + "### Квантиль нормального распределения:\n", + "Для уровня значимости $\\alpha/2$ (например, $\\alpha = 0.1$ для 90% доверительного интервала) квантиль:\n", + "$$\n", + "z_{1-\\alpha/2} \\approx 1.6449 \\quad (z_{0.95}).\n", + "$$\n", + "\n", + "\n", + "### Доверительный интервал:\n", + "Интервал для уровня доверия $1-\\alpha$:\n", + "$$\n", + "\\hat{\\lambda} \\pm z_{1-\\alpha/2} \\cdot \\text{SE}(\\hat{\\lambda}),\n", + "$$\n", + "что дает границы:\n", + "$$\n", + "\\left( \\hat{\\lambda} - 1.6449 \\cdot \\frac{\\hat{\\lambda}}{\\sqrt{n}}, \\quad \\hat{\\lambda} + 1.6449 \\cdot \\frac{\\hat{\\lambda}}{\\sqrt{n}} \\right)." ] }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 17, "id": "7f3db200", "metadata": {}, "outputs": [ @@ -596,7 +634,7 @@ }, { "cell_type": "code", - "execution_count": 43, + "execution_count": 18, "id": "d881725f", "metadata": {}, "outputs": [ @@ -649,20 +687,7 @@ "# # Критическое значение для alpha=0.10\n", "# critical_value = 1.22 / np.sqrt(n) \n", "\n", - "# print(f\"Статистика D: {D:.4f}, Критическое значение: {critical_value:.4f}\")\n", - "# print(f\" Гипотеза {'отвергается' if D > critical_value else 'не отвергается'}\")\n", - "# print(f\" p-значение: {p_value:.4f}\")\n", "\n", - "# # Визуализация\n", - "# plt.figure(figsize=(10, 6))\n", - "# plt.step(x_sorted, empirical_cdf, where='post', label='Эмпирическая ФР')\n", - "# plt.plot(x_sorted, theoretical_cdf, 'r-', label=f'Теоретическая ФР (λ={lambda0})')\n", - "# plt.xlabel('x')\n", - "# plt.ylabel('F(x)')\n", - "# plt.title('Сравнение эмпирической и теоретической функций распределения')\n", - "# plt.legend()\n", - "# plt.grid(True)\n", - "# plt.show()\n", "\n", "# e) Критерий Колмогорова-Смирнова с использованием scipy\n", "from scipy.stats import kstest\n", @@ -678,7 +703,13 @@ "print(f\"\\nКритерий Колмогорова-Смирнова:\")\n", "print(f\" Статистика Dn: {D_stat:.4f}, Критическое значение: {critical_value:.4f}\")\n", "print(f\" P-value: {p_value:.4f}\")\n", - "print(f\" Гипотеза {'отвергается' if D_stat > critical_value else 'не отвергается'}\")" + "print(f\" Гипотеза {'отвергается' if D_stat > critical_value else 'не отвергается'}\")\n", + "\n", + "\n", + "x_sorted = np.sort(data)\n", + "n = len(data)\n", + "empirical_cdf = np.arange(1, n + 1) / n\n", + "theoretical_cdf = 1 - np.exp(-lambda0 * x_sorted)\n" ] }, { @@ -696,7 +727,7 @@ }, { "cell_type": "code", - "execution_count": 53, + "execution_count": 19, "id": "4383629c", "metadata": {}, "outputs": [ @@ -704,17 +735,27 @@ "name": "stdout", "output_type": "stream", "text": [ - "χ² статистика: 14.4669\n", - "Критическое значение (α=0.1): 6.2514\n", - "p-значение: 0.002334\n", - "Степени свободы: 3\n", - "Гипотеза отвергается на уровне 0.1\n" + "\n", + "Таблица частот после объединения:\n", + "Интервал \tНаблюдаемо\tОжидаемо\tВклад в χ²\n", + "[0.00, 2.80)\t\t 34 \t30.15\t\t0.4906497\n", + "[2.80, 14.00)\t\t 7 \t11.97\t\t2.0627840\n", + "[14.00, 15.40)\t\t 9 \t7.68\t\t0.2261086\n", + "χ² статистика: 2.7795\n", + "Критическое значение (α=0.1): 4.6052\n", + "p-значение: 0.249132\n", + "Степени свободы: 2\n", + "Нет оснований отвергнуть гипотезу на уровне 0.1\n" ] } ], "source": [ + "# Теоретическое смещение для больших выборок\n", + "# Для показательного распределения смещение ОМП приближенно равно λ/n\n", "from scipy.stats import chi2, expon\n", "\n", + "n = len(data)\n", + "\n", "# Гистограмма и ожидаемые частоты\n", "observed, bin_edges = np.histogram(data, bins=bins)\n", "expected = []\n", @@ -725,17 +766,20 @@ " expected.append(prob * n)\n", "expected = np.array(expected)\n", "\n", - "# Объединение интервалов с expected < 5\n", - "while np.any(expected < 5):\n", - " min_idx = np.argmin(expected)\n", + "# Объединение интервалов\n", + "while np.any(observed < 5):\n", + " min_idx = np.argmin(observed)\n", " if min_idx == 0:\n", " # Объединяем с следующим интервалом\n", " expected[min_idx + 1] += expected[min_idx]\n", " observed[min_idx + 1] += observed[min_idx]\n", + " # bin_edges = np.delete(bin_edges, min_idx)\n", " else:\n", " # Объединяем с предыдущим интервалом\n", " expected[min_idx - 1] += expected[min_idx]\n", " observed[min_idx - 1] += observed[min_idx]\n", + " bin_edges = np.delete(bin_edges, min_idx-1)\n", + " # print(bin_edges)\n", " # Удаляем объединенный интервал\n", " expected = np.delete(expected, min_idx)\n", " observed = np.delete(observed, min_idx)\n", @@ -746,6 +790,14 @@ "chi2_crit = chi2.ppf(1 - alpha2, dof)\n", "p_value = 1 - chi2.cdf(chi2_stat, dof)\n", "\n", + "print(\"\\nТаблица частот после объединения:\")\n", + "print(\"Интервал \\tНаблюдаемо\\tОжидаемо\\tВклад в χ²\")\n", + "for i in range(len(observed)):\n", + " lower = bin_edges[i]\n", + " upper = bin_edges[i+1]\n", + " contrib = (observed[i]-expected[i])**2/expected[i]\n", + " print(f\"[{lower:.2f}, {upper:.2f})\\t\\t{observed[i]:^10}\\t{expected[i]:.2f}\\t\\t{contrib:.7f}\")\n", + "\n", "# Вывод результатов\n", "print(f\"χ² статистика: {chi2_stat:.4f}\")\n", "print(f\"Критическое значение (α={alpha2}): {chi2_crit:.4f}\")\n", @@ -754,7 +806,7 @@ "if chi2_stat > chi2_crit:\n", " print(f\"Гипотеза отвергается на уровне {alpha2}\")\n", "else:\n", - " print(f\"Нет оснований отвергнуть гипотезу на уровне {alpha2}\")" + " print(f\"Нет оснований отвергнуть гипотезу на уровне {alpha2}\")\n" ] }, { @@ -768,7 +820,7 @@ }, { "cell_type": "code", - "execution_count": 63, + "execution_count": 20, "id": "7456dc2c", "metadata": {}, "outputs": [], @@ -836,7 +888,7 @@ }, { "cell_type": "code", - "execution_count": 62, + "execution_count": 22, "id": "60c8dfd9", "metadata": {}, "outputs": [ @@ -844,20 +896,20 @@ "name": "stdout", "output_type": "stream", "text": [ - "Оценка λ: 0.3586\n", + "Оценка λ: 0.4490\n", + "\n", + "Таблица частот после объединения:\n", + "Интервал \tНаблюдаемо\tОжидаемо\tВклад в χ²\n", + "[0.00, 1.40)\t\t 30 \t23.33\t\t1.9049\n", + "[1.40, 4.20)\t\t 10 \t19.08\t\t4.3222\n", + "[4.20, 16.80)\t\t 10 \t7.56\t\t0.7883\n", "\n", "Критерий χ² для сложной гипотезы:\n", - " Статистика χ²: 10.9186\n", - " Критическое значение (α=0.1): 4.6052\n", - " p-значение: 0.0043\n", - " Степени свободы: 2\n", - "Гипотеза отвергается на уровне 0.1\n", - "\n", - "Таблица частот:\n", - "[0.00, 1.40): O=30, E=19.74\n", - "[1.40, 2.80): O=4, E=11.95\n", - "[2.80, 4.20): O=6, E=7.23\n", - "[4.20, 16.80): O=10, E=10.97\n" + " Статистика χ²: 7.0154\n", + " Критическое значение (α=0.1): 2.7055\n", + " p-значение: 0.0081\n", + " Степени свободы: 1\n", + "Гипотеза отвергается на уровне 0.1\n" ] } ], @@ -865,19 +917,8 @@ "import numpy as np\n", "from scipy.stats import chi2, expon\n", "\n", - "# Данные\n", - "data = np.array([\n", - " 0.18, 0.10, 3.34, 0.67, 0.85, 1.17, 0.24, 0.15, 1.31, 0.00, 0.49, 2.37, 14.94, 2.44, 3.13, 0.06, 2.98, 9.25, 6.84, 3.96,\n", - " 0.07, 6.72, 11.83, 0.50, 0.11, 6.50, 0.29, 0.17, 0.03, 0.06, 1.02, 0.49, 15.68, 3.03, 0.24, 11.40, 0.53, 0.59, 4.55, 3.57,\n", - " 8.33, 0.12, 2.58, 2.77, 0.12, 1.11, 0.31, 0.36, 1.31, 0.57\n", - "])\n", - "\n", - "alpha2 = 0.10\n", - "h = 1.40\n", - "n = len(data)\n", - "\n", "# Оценка параметра λ методом максимального правдоподобия\n", - "lambda_ml = 1 / np.mean(data)\n", + "lambda_ml = 0.4490\n", "print(f\"Оценка λ: {lambda_ml:.4f}\")\n", "\n", "# Построение интервалов с шагом h\n", @@ -900,9 +941,9 @@ "bin_edges_list = bin_edges.tolist()\n", "\n", "i = 0\n", - "while i < len(expected_list):\n", - " if expected_list[i] < 5:\n", - " if i == 0 or i == len(expected_list) - 1:\n", + "while i < len(observed_list):\n", + " if observed_list[i] <= 5:\n", + " if i == 0 or i == len(observed_list) - 1:\n", " # Объединяем с соседним интервалом\n", " if i == 0:\n", " expected_list[i+1] += expected_list[i]\n", @@ -941,6 +982,15 @@ "chi2_crit = chi2.ppf(1 - alpha2, df)\n", "p_value = 1 - chi2.cdf(chi2_stat, df)\n", "\n", + "print(\"\\nТаблица частот после объединения:\")\n", + "print(\"Интервал \\tНаблюдаемо\\tОжидаемо\\tВклад в χ²\")\n", + "for i in range(len(observed_list)):\n", + " lower = bin_edges_list[i]\n", + " upper = bin_edges_list[i+1]\n", + " contrib = (observed_list[i]-expected_list[i])**2/expected_list[i]\n", + " print(f\"[{lower:.2f}, {upper:.2f})\\t\\t{observed_list[i]:^10}\\t{expected_list[i]:.2f}\\t\\t{contrib:.4f}\")\n", + "\n", + "\n", "# Вывод результатов\n", "print(f\"\\nКритерий χ² для сложной гипотезы:\")\n", "print(f\" Статистика χ²: {chi2_stat:.4f}\")\n", @@ -951,12 +1001,113 @@ "if chi2_stat > chi2_crit:\n", " print(f\"Гипотеза отвергается на уровне {alpha2}\")\n", "else:\n", - " print(f\"Нет оснований отвергнуть гипотезу на уровне {alpha2}\")\n", + " print(f\"Нет оснований отвергнуть гипотезу на уровне {alpha2}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "28fd5d5e", + "metadata": {}, + "outputs": [], + "source": [ + "# import numpy as np\n", + "# from scipy.stats import chi2, expon\n", "\n", - "# Вывод таблицы с интервалами\n", - "print(\"\\nТаблица частот:\")\n", - "for i in range(len(bin_edges_list) - 1):\n", - " print(f\"[{bin_edges_list[i]:.2f}, {bin_edges_list[i+1]:.2f}): O={observed_list[i]}, E={expected_list[i]:.2f}\")" + "\n", + "# # Функция для вычисления χ² и объединения интервалов\n", + "# def calculate_chi2(lambda_, data, h):\n", + "# n = len(data)\n", + "# max_val = np.max(data)\n", + "# bins = np.arange(0, max_val + h, h)\n", + "# if bins[-1] < max_val:\n", + "# bins = np.append(bins, bins[-1] + h)\n", + "# observed, _ = np.histogram(data, bins=bins)\n", + " \n", + "# # Теоретические вероятности и ожидаемые частоты\n", + "# probs = []\n", + "# for i in range(len(bins) - 1):\n", + "# cdf_lower = expon.cdf(bins[i], scale=1/lambda_)\n", + "# cdf_upper = expon.cdf(bins[i+1], scale=1/lambda_)\n", + "# prob = cdf_upper - cdf_lower\n", + "# probs.append(prob)\n", + "# expected = np.array(probs) * n\n", + " \n", + "# # Объединение интервалов с expected < 5\n", + "# observed_list = observed.tolist()\n", + "# expected_list = expected.tolist()\n", + "# bin_edges_list = bins.tolist()\n", + " \n", + "# i = 0\n", + "# while i < len(expected_list):\n", + "# if expected_list[i] < 5:\n", + "# if i == len(expected_list) - 1: # Последний интервал\n", + "# expected_list[i-1] += expected_list[i]\n", + "# observed_list[i-1] += observed_list[i]\n", + "# del expected_list[i]\n", + "# del observed_list[i]\n", + "# del bin_edges_list[i+1]\n", + "# i -= 1 # Проверить новый интервал\n", + "# else: # Объединить со следующим интервалом\n", + "# expected_list[i] += expected_list[i+1]\n", + "# observed_list[i] += observed_list[i+1]\n", + "# del expected_list[i+1]\n", + "# del observed_list[i+1]\n", + "# del bin_edges_list[i+1]\n", + "# else:\n", + "# i += 1\n", + " \n", + "# # Вычисление статистики χ²\n", + "# chi2_val = 0\n", + "# for obs, exp in zip(observed_list, expected_list):\n", + "# if exp > 0: # Избегаем деления на ноль\n", + "# chi2_val += (obs - exp) ** 2 / exp\n", + "# return chi2_val, bin_edges_list, observed_list, expected_list\n", + "\n", + "# # Поиск λ, минимизирующей χ² в диапазоне [0.2, 0.5]\n", + "# lambda_grid = np.arange(0.2, 0.51, 0.001)\n", + "# min_chi2 = float('inf')\n", + "# lambda_min = None\n", + "# results = {}\n", + "\n", + "# for lambda_val in lambda_grid:\n", + "# chi2_val, bin_edges, observed, expected = calculate_chi2(lambda_val, data, h)\n", + "# results[lambda_val] = (chi2_val, bin_edges, observed, expected)\n", + "# if chi2_val < min_chi2:\n", + "# min_chi2 = chi2_val\n", + "# lambda_min = lambda_val\n", + "\n", + "# # Финализация результатов для оптимального λ\n", + "# chi2_stat, bin_edges_final, observed_final, expected_final = results[lambda_min]\n", + "# k = len(observed_final)\n", + "# df = k - 2 # Степени свободы (k интервалов - 1 - 1 параметр)\n", + "# chi2_crit = chi2.ppf(1 - alpha2, df)\n", + "# p_value = 1 - chi2.cdf(chi2_stat, df)\n", + "\n", + "# # Вывод результатов\n", + "# print(f\"Оптимальное λ, минимизирующее χ²: {lambda_min:.10f}\")\n", + "# print(f\"Минимальная статистика χ²: {chi2_stat:.4f}\")\n", + "# print(f\"Степени свободы: {df}\")\n", + "# print(f\"Критическое значение χ² (α={alpha2}): {chi2_crit:.4f}\")\n", + "# print(f\"p-значение: {p_value:.4f}\")\n", + "\n", + "# # Проверка гипотезы\n", + "# if chi2_stat > chi2_crit:\n", + "# print(f\"Гипотеза отвергается на уровне {alpha2}\")\n", + "# else:\n", + "# print(f\"Нет оснований отвергнуть гипотезу на уровне {alpha2}\")\n", + "\n", + "# # Наибольший уровень значимости, при котором гипотеза не отвергается\n", + "# print(f\"Наибольший уровень значимости: {p_value:.4f}\")\n", + "\n", + "# # Вывод таблицы частот\n", + "# print(\"\\nИнтервалы и частоты после объединения:\")\n", + "# print(\"Интервал \\tНаблюдаемо\\tОжидаемо\\tВклад в χ²\")\n", + "# for i in range(len(observed_final)):\n", + "# lower = bin_edges_final[i]\n", + "# upper = bin_edges_final[i+1]\n", + "# contrib = (observed_final[i] - expected_final[i])**2 / expected_final[i] if expected_final[i] > 0 else 0\n", + "# print(f\"[{lower:.2f}, {upper:.2f})\\t{observed_final[i]:10}\\t{expected_final[i]:10.2f}\\t{contrib:.4f}\")" ] }, { diff --git a/idz4/ИДЗ 4_1 Артём.ipynb b/idz4/ИДЗ 4_1 Артём.ipynb index 9a99f67..7d2cbc4 100644 --- a/idz4/ИДЗ 4_1 Артём.ipynb +++ b/idz4/ИДЗ 4_1 Артём.ipynb @@ -541,7 +541,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "id": "72be1710", "metadata": {}, "outputs": [ @@ -588,7 +588,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 23, "id": "bd170677", "metadata": {}, "outputs": [ @@ -596,50 +596,179 @@ "name": "stdout", "output_type": "stream", "text": [ - "Хи-квадрат статистика: 2.7737\n", - "Критическое значение: 13.3882\n", - "p-value: 0.7348\n", - "Не отвергаем H0: распределение нормальное\n" + "Оценка μ (среднее остатков): -0.0000\n", + "Оценка σ (стандартное отклонение отстатков): 2.3980\n", + "[ 5 8 15 8 14] [ 5.312141 8.34056556 11.88152698 11.5684878 12.28662505]\n", + "\n", + "Результаты после объединения интервалов:\n", + "Интервалы\t\tНаблюдаемо\tОжидаемо\tВклад в χ²\n", + "[-5.9095, -2.9020]\t 5 \t5.31\t\t0.0183\n", + "[-2.9020, -1.3982]\t 8 \t8.34\t\t0.0139\n", + "[-1.3982, 0.1055]\t 15 \t11.88\t\t0.8185\n", + "[0.1055, 1.6093]\t 8 \t11.57\t\t1.1008\n", + "[1.6093, 6.1205]\t 14 \t12.29\t\t0.2389\n", + "\n", + "Количество интервалов после объединения: 5\n", + "Степени свободы: 2 (5 интервалов - 1 - 2 оцененных параметра)\n", + "Статистика χ²: 2.1904\n", + "Критическое значение (α=0.02): 7.8240\n", + "p-value: 0.3345\n", + "Не отвергаем H₀: распределение нормальное\n" ] } ], "source": [ - "# Разбиение на интервалы (используем 6 интервалов для примера)\n", + "import numpy as np\n", + "from scipy.stats import chi2, norm\n", + "# Вычисление параметров нормального распределения\n", "mu = np.mean(residuals) # Среднее остатков\n", "std = np.std(residuals, ddof=1) # Стандартное отклонение (несмещенное)\n", + "print(f\"Оценка μ (среднее остатков): {mu:.4f}\")\n", + "print(f\"Оценка σ (стандартное отклонение отстатков): {std:.4f}\")\n", + "\n", + "# Разбиение на интервалы\n", "observed_freq, bins = np.histogram(residuals, bins=8)\n", + "n = len(residuals)\n", "n_bins = len(observed_freq)\n", "\n", - "# Ожидаемые частоты для нормального распределения\n", - "expected_freq = []\n", + "# Вычисление ожидаемых частот\n", + "expected_freq = np.zeros(n_bins)\n", "for i in range(n_bins):\n", - " bin_start = bins[i]\n", - " bin_end = bins[i+1]\n", - " cdf_start = norm.cdf(bin_start, mu, std)\n", - " cdf_end = norm.cdf(bin_end, mu, std)\n", - " expected_freq.append(len(residuals) * (cdf_end - cdf_start))\n", + " cdf_start = norm.cdf(bins[i], mu, std)\n", + " cdf_end = norm.cdf(bins[i+1], mu, std)\n", + " expected_freq[i] = n * (cdf_end - cdf_start)\n", + "\n", + "# Объединение интервалов с expected_freq < 5\n", + "observed_list = observed_freq.tolist()\n", + "expected_list = expected_freq.tolist()\n", + "bins_list = bins.tolist()\n", + "\n", + "# i = 0\n", + "# while i < len(expected_list):\n", + "# if expected_list[i] < 5:\n", + "# # Если последний интервал, объединить с предыдущим\n", + "# if i == len(expected_list) - 1:\n", + "# if i > 0: # Если есть предыдущий интервал\n", + "# expected_list[i-1] += expected_list[i]\n", + "# observed_list[i-1] += observed_list[i]\n", + "# del expected_list[i]\n", + "# del observed_list[i]\n", + "# del bins_list[i+1]\n", + "# i -= 1 # Вернуться к новому интервалу\n", + "# else:\n", + "# # Объединить с следующим интервалом\n", + "# expected_list[i] += expected_list[i+1]\n", + "# observed_list[i] += observed_list[i+1]\n", + "# del expected_list[i+1]\n", + "# del observed_list[i+1]\n", + "# del bins_list[i+1]\n", + "# else:\n", + "# i += 1\n", + "i = 0\n", + "while i < len(expected_list):\n", + " if expected_list[i] < 5:\n", + " if i == len(expected_list) - 1: # последний интервал\n", + " if i > 0:\n", + " expected_list[i-1] += expected_list[i]\n", + " observed_list[i-1] += observed_list[i]\n", + " del expected_list[i]\n", + " del observed_list[i]\n", + " del bins_list[i] # ← исправлено\n", + " i -= 1\n", + " else: # объединяем с правым\n", + " expected_list[i] += expected_list[i+1]\n", + " observed_list[i] += observed_list[i+1]\n", + " del expected_list[i+1]\n", + " del observed_list[i+1]\n", + " del bins_list[i+1] # граница между i и i+1\n", + " else:\n", + " i += 1\n", + "\n", + "\n", + "# Проверка минимального количества интервалов\n", + "if len(expected_list) < 3:\n", + " # Если интервалов меньше 3, принудительно объединяем все интервалы\n", + " expected_list = [sum(expected_list)]\n", + " observed_list = [sum(observed_list)]\n", + " bins_list = [bins[0], bins[-1]]\n", + " print(\"Внимание: после объединения остался только 1 интервал!\")\n", + "\n", + "# Преобразование в массивы NumPy\n", + "observed_merged = np.array(observed_list)\n", + "expected_merged = np.array(expected_list)\n", + "\n", + "# Расчет статистики хи-квадрат\n", + "chi2_stat = np.sum((observed_merged - expected_merged)**2 / expected_merged)\n", + "\n", + "# Степени свободы\n", + "k = len(observed_merged) # Количество интервалов после объединения\n", + "dof = k - 1 - 2 # k-1 интервалов минус 2 оцененных параметра (μ и σ)\n", + "if dof < 1:\n", + " dof = 1 # Минимум 1 степень свободы\n", "\n", - "# Критерий хи-квадрат\n", - "chi2_stat = sum((observed_freq - expected_freq)**2 / expected_freq)\n", - "dof = n_bins - 1 - 2 # 2 параметра (mu, std) оценены по данным\n", - "alpha = 0.02\n", "critical_value = chi2.ppf(1 - alpha, dof)\n", "p_value = 1 - chi2.cdf(chi2_stat, dof)\n", + "print(observed_merged, expected_merged)\n", + "# Вывод результатов\n", + "print(\"\\nРезультаты после объединения интервалов:\")\n", + "print(\"Интервалы\\t\\tНаблюдаемо\\tОжидаемо\\tВклад в χ²\")\n", "\n", - "print(f\"Хи-квадрат статистика: {chi2_stat:.4f}\")\n", - "print(f\"Критическое значение: {critical_value:.4f}\")\n", + "for i in range(len(observed_merged)):\n", + " lower = bins_list[i]\n", + " upper = bins_list[i+1]\n", + " contrib = (observed_merged[i] - expected_merged[i])**2 / expected_merged[i]\n", + " print(f\"[{lower:.4f}, {upper:.4f}]\\t{observed_merged[i]:^12}\\t{expected_merged[i]:.2f}\\t\\t{contrib:.4f}\")\n", + "\n", + "print(f\"\\nКоличество интервалов после объединения: {k}\")\n", + "print(f\"Степени свободы: {dof} (5 интервалов - 1 - 2 оцененных параметра)\")\n", + "print(f\"Статистика χ²: {chi2_stat:.4f}\")\n", + "print(f\"Критическое значение (α={alpha}): {critical_value:.4f}\")\n", "print(f\"p-value: {p_value:.4f}\")\n", "\n", - "# Визуальная оценка\n", + "# Проверка гипотезы\n", "if chi2_stat > critical_value:\n", - " print(\"Отвергаем H0: распределение не нормальное\")\n", + " print(\"Отвергаем H₀: распределение не нормальное\")\n", "else:\n", - " print(\"Не отвергаем H0: распределение нормальное\")" + " print(\"Не отвергаем H₀: распределение нормальное\")\n", + "\n", + "\n", + "# # Разбиение на интервалы (используем 6 интервалов для примера)\n", + "# mu = np.mean(residuals) # Среднее остатков\n", + "# std = np.std(residuals, ddof=1) # Стандартное отклонение (несмещенное)\n", + "# observed_freq, bins = np.histogram(residuals, bins=8)\n", + "# n_bins = len(observed_freq)\n", + "\n", + "# # Ожидаемые частоты для нормального распределения\n", + "# expected_freq = []\n", + "# for i in range(n_bins):\n", + "# bin_start = bins[i]\n", + "# bin_end = bins[i+1]\n", + "# cdf_start = norm.cdf(bin_start, mu, std)\n", + "# cdf_end = norm.cdf(bin_end, mu, std)\n", + "# expected_freq.append(len(residuals) * (cdf_end - cdf_start))\n", + "\n", + "# # Критерий хи-квадрат\n", + "# chi2_stat = sum((observed_freq - expected_freq)**2 / expected_freq)\n", + "# dof = n_bins - 1 - 2 # 2 параметра (mu, std) оценены по данным\n", + "# alpha = 0.02\n", + "# critical_value = chi2.ppf(1 - alpha, dof)\n", + "# p_value = 1 - chi2.cdf(chi2_stat, dof)\n", + "# print(observed_freq, expected_freq);\n", + "\n", + "# print(f\"Хи-квадрат статистика: {chi2_stat:.4f}\")\n", + "# print(f\"Критическое значение: {critical_value:.4f}\")\n", + "# print(f\"p-value: {p_value:.4f}\")\n", + "\n", + "# # Визуальная оценка\n", + "# if chi2_stat > critical_value:\n", + "# print(\"Отвергаем H0: распределение не нормальное\")\n", + "# else:\n", + "# print(\"Не отвергаем H0: распределение нормальное\")" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 11, "id": "f498c322", "metadata": {}, "outputs": [], @@ -731,7 +860,7 @@ }, { "cell_type": "code", - "execution_count": 69, + "execution_count": 12, "id": "ca6842f7", "metadata": {}, "outputs": [ @@ -776,7 +905,7 @@ }, { "cell_type": "code", - "execution_count": 86, + "execution_count": 13, "id": "68365ffd", "metadata": {}, "outputs": [ @@ -836,7 +965,7 @@ }, { "cell_type": "code", - "execution_count": 90, + "execution_count": 14, "id": "76e0b82d", "metadata": {}, "outputs": [], @@ -868,7 +997,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 15, "id": "f791c572", "metadata": {}, "outputs": [ @@ -933,7 +1062,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 16, "id": "9b48da35", "metadata": {}, "outputs": [ @@ -965,7 +1094,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 17, "id": "b34812e2", "metadata": {}, "outputs": [], @@ -1030,13 +1159,13 @@ "- $H_1$: Зависимость нелинейна ($\\beta_3 \\neq 0$).\n", "\n", "#### Гипотеза независимости\n", - "- $H_0$: $Y$ не зависит от $X$ линейна ($\\beta_2 = \\beta_3 = 0$).\n", - "- $H_1$: $Y$ зависит от $X$ линейна (хотя бы один из $\\beta_2, \\beta_3 \\neq 0$)." + "- $H_0$: $Y$ не зависит от $X$ линейно ($\\beta_2 = \\beta_3 = 0$).\n", + "- $H_1$: $Y$ зависит от $X$ линейно (хотя бы один из $\\beta_2, \\beta_3 \\neq 0$)." ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 18, "id": "1fde6d40", "metadata": {}, "outputs": [], @@ -1070,7 +1199,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 19, "id": "405456a9", "metadata": {}, "outputs": [ @@ -1080,55 +1209,67 @@ "text": [ "\n", "Проверка гипотез:\n", - "Проверка гипотезы линейности (H₀: β₃ = 0):\n", - "t-статистика: 0.6775\n", - "p-значение: 0.5014\n", - "Нет оснований отвергать гипотезу о линейности\n", "\n", - "Проверка гипотезы независимости (H₀: β₂ = 0):\n", - "t-статистика: -0.8509\n", - "p-значение: 0.3991\n", - "Нет оснований отвергать гипотезу о независимости\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "C:\\Users\\margaery\\AppData\\Local\\Temp\\ipykernel_15852\\3529389018.py:4: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", - " print(f\"t-статистика: {model_poly.tvalues[2]:.4f}\")\n", - "C:\\Users\\margaery\\AppData\\Local\\Temp\\ipykernel_15852\\3529389018.py:5: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", - " print(f\"p-значение: {model_poly.pvalues[2]:.4f}\")\n", - "C:\\Users\\margaery\\AppData\\Local\\Temp\\ipykernel_15852\\3529389018.py:6: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", - " if model_poly.pvalues[2] < alpha:\n", - "C:\\Users\\margaery\\AppData\\Local\\Temp\\ipykernel_15852\\3529389018.py:13: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", - " print(f\"t-статистика: {model_poly.tvalues[1]:.4f}\")\n", - "C:\\Users\\margaery\\AppData\\Local\\Temp\\ipykernel_15852\\3529389018.py:14: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", - " print(f\"p-значение: {model_poly.pvalues[1]:.4f}\")\n", - "C:\\Users\\margaery\\AppData\\Local\\Temp\\ipykernel_15852\\3529389018.py:15: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", - " if model_poly.pvalues[1] < alpha:\n" + "Проверка гипотезы линейности (H₀: β₃ = 0):\n", + "F-статистика: 0.4591\n", + "p-значение: 0.5014\n", + "Нет оснований отвергать гипотезу о линейности (alpha=0.02)\n", + "\n", + "Проверка гипотезы независимости (H₀: β₂ = β₃ = 0):\n", + "F-статистика: 0.5770\n", + "p-значение: 0.5655\n", + "Нет оснований отвергать гипотезу о независимости (alpha=0.02)\n" ] } ], "source": [ - "print(\"\\nПроверка гипотез:\")\n", - "# Тест на линейность (значимость β₃)\n", - "print(\"Проверка гипотезы линейности (H₀: β₃ = 0):\")\n", - "print(f\"t-статистика: {model_poly.tvalues[2]:.4f}\")\n", - "print(f\"p-значение: {model_poly.pvalues[2]:.4f}\")\n", - "if model_poly.pvalues[2] < alpha:\n", - " print(f\"Гипотеза о линейности отвергается\")\n", - "else:\n", - " print(f\"Нет оснований отвергать гипотезу о линейности\")\n", + "model_const = sm.OLS(df['Y'], sm.add_constant(np.ones(len(df)))).fit()\n", "\n", - "# Тест на независимость (значимость β₂)\n", - "print(\"\\nПроверка гипотезы независимости (H₀: β₂ = 0):\")\n", - "print(f\"t-статистика: {model_poly.tvalues[1]:.4f}\")\n", - "print(f\"p-значение: {model_poly.pvalues[1]:.4f}\")\n", - "if model_poly.pvalues[1] < alpha:\n", - " print(f\"Гипотеза о независимости отвергается\")\n", + "# print(\"\\nПроверка гипотез:\")\n", + "# # Тест на линейность (значимость β₃)\n", + "# print(\"Проверка гипотезы линейности (H₀: β₃ = 0):\")\n", + "# print(f\"t-статистика: {model_poly.tvalues[2]:.4f}\")\n", + "# print(f\"p-значение: {model_poly.pvalues[2]:.4f}\")\n", + "# if model_poly.pvalues[2] < alpha:\n", + "# print(f\"Гипотеза о линейности отвергается\")\n", + "# else:\n", + "# print(f\"Нет оснований отвергать гипотезу о линейности\")\n", + "\n", + "# # Тест на независимость (значимость β₂)\n", + "# print(\"\\nПроверка гипотезы независимости (H₀: β₂ = 0):\")\n", + "# print(f\"t-статистика: {model_poly.tvalues[1]:.4f}\")\n", + "# print(f\"p-значение: {model_poly.pvalues[1]:.4f}\")\n", + "# if model_poly.pvalues[1] < alpha:\n", + "# print(f\"Гипотеза о независимости отвергается\")\n", + "# else:\n", + "# print(f\"Нет оснований отвергать гипотезу о независимости\")\n", + "\n", + "from statsmodels.stats.anova import anova_lm\n", + "\n", + "print(\"\\nПроверка гипотез:\")\n", + "\n", + "# Сравнение моделей для проверки гипотез\n", + "# 1. Сравнение линейной и квадратичной моделей (проверка линейности)\n", + "anova_results_lin = anova_lm(model_lin, model_poly)\n", + "print(\"\\nПроверка гипотезы линейности (H₀: β₃ = 0):\")\n", + "print(f\"F-статистика: {anova_results_lin.F[1]:.4f}\")\n", + "print(f\"p-значение: {anova_results_lin['Pr(>F)'][1]:.4f}\")\n", + "\n", + "if anova_results_lin['Pr(>F)'][1] < alpha:\n", + " print(\"Гипотеза о линейности отвергается (alpha=0.02)\")\n", "else:\n", - " print(f\"Нет оснований отвергать гипотезу о независимости\")\n" + " print(\"Нет оснований отвергать гипотезу о линейности (alpha=0.02)\")\n", + "\n", + "# 2. Сравнение константной и квадратичной моделей (проверка независимости)\n", + "anova_results_const = anova_lm(model_const, model_poly)\n", + "print(\"\\nПроверка гипотезы независимости (H₀: β₂ = β₃ = 0):\")\n", + "print(f\"F-статистика: {anova_results_const.F[1]:.4f}\")\n", + "print(f\"p-значение: {anova_results_const['Pr(>F)'][1]:.4f}\")\n", + "\n", + "if anova_results_const['Pr(>F)'][1] < alpha:\n", + " print(\"Гипотеза о независимости отвергается (alpha=0.02)\")\n", + "else:\n", + " print(\"Нет оснований отвергать гипотезу о независимости (alpha=0.02)\")" ] }, { @@ -1153,7 +1294,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 21, "id": "c00ff024", "metadata": {}, "outputs": [ @@ -1166,20 +1307,21 @@ "--------------------------------------\n", "Модель AIC BIC\n", "Линейная 232.83 236.66\n", - "Квадратичная 234.35 240.08\n" + "Квадратичная 234.35 240.08\n", + "Константная 231.56 233.47\n" ] } ], "source": [ "# f) AIC и BIC\n", "# Добавляем константную модель для сравнения\n", - "model_const = sm.OLS(df['Y'], sm.add_constant(np.ones(len(df)))).fit()\n", "\n", "print(\"\\nСравнение моделей по AIC и BIC:\")\n", "print(\"--------------------------------------\")\n", "print(\"Модель AIC BIC\")\n", "print(f\"Линейная {model_lin.aic:.2f} {model_lin.bic:.2f}\")\n", - "print(f\"Квадратичная {model_poly.aic:.2f} {model_poly.bic:.2f}\")\n" + "print(f\"Квадратичная {model_poly.aic:.2f} {model_poly.bic:.2f}\")\n", + "print(f\"Константная {model_const.aic:.2f} {model_const.bic:.2f}\")\n" ] }, { @@ -1187,57 +1329,8 @@ "id": "d66aad56", "metadata": {}, "source": [ - "**AIC/BIC** линейной модели меньше, она лучше описывает данные." - ] - }, - { - "cell_type": "markdown", - "id": "a6887b63", - "metadata": {}, - "source": [ - "## Пункт g)\n", - "### Характер зависимости $Y$ от $X$\n", - "- **Линейная модель:**\n", - " $$\n", - " Y = 15.59 - 0.25X,\\ R^2 = 0.014.\n", - " $$\n", - " - Крайне низкий $R^2$ (1.4%) указывает на отсутствие линейной зависимости.\n", - " - Коэффициент $\\beta_2 = -0.25$ статистически незначим (доверительный интервал [−4.29, 2.05] включает ноль).\n", - "\n", - "- **Квадратичная модель:**\n", - " $$\n", - " Y = 16.87 - 1.12X + 0.13X^2,\\ R^2 = 0.024.\n", - " $$\n", - " - $R^2 = 2.4\\%$ показывает, что модель объясняет лишь незначительную часть вариации.\n", - " - Коэффициенты:\n", - " - $\\beta_2 = -1.12$ (линейный член): интервал [−4.29, 2.05] включает ноль.\n", - " - $\\beta_3 = 0.13$ (квадратичный член): интервал [−0.33, 0.59] включает ноль.\n", - "\n", - "### Проверка гипотез\n", - "Остатки близки к нормальному распределению. Критерий $\\chi^2$ не выявил значимых отклонений от нормальности на уровне $\\alpha=0.02$.\n", - "\n", - "*Предположение о нормальности ошибок выполняется.*\n", - "\n", - "### AIC/BIC\n", - "| Модель | AIC | BIC |\n", - "|----------------|--------|--------|\n", - "| Линейная | 232.83 | 236.66 |\n", - "| Квадратичная | 234.35 | 240.08 |\n", - "\n", - "- **Линейная модель** имеет более низкие AIC/BIC, чем квадратичная.\n", - "\n", - "### Аномалии в результатах\n", - "\n", - "**Парадокс низкого $R^2$:**\n", - " - Обе модели объясняют менее 3% вариации, что ставит под сомнение их практическую применимость.\n", - "\n", - "### Итоговый вывод\n", - "- **Отсутствие значимой связи:** Ни линейная, ни квадратичная модели не демонстрируют статистически значимой зависимости $Y$ от $X$ на уровне $\\alpha=0.02$.\n", - "- **Артефакты анализа:** Низкий $R^2$, незначимые коэффициенты и противоречивые результаты тестов указывают на то, что переменная $X$ не является релевантным предиктором для $Y$ в данном наборе данных.\n", - "- **Рекомендации:** \n", - " - Проверить данные на наличие выбросов или ошибок.\n", - " - Рассмотреть другие предикторы или преобразования.\n", - " - Увеличить объем данных для повышения надежности тестов." + "**AIC/BIC** константной модели меньше, она лучше описывает данные. \n", + "Это согласуется с F-тестами: константная модель лучше описывает данные" ] } ], diff --git a/idz4/ИДЗ 4_2 Артём.ipynb b/idz4/ИДЗ 4_2 Артём.ipynb index 0cc9f55..d59eb7f 100644 --- a/idz4/ИДЗ 4_2 Артём.ipynb +++ b/idz4/ИДЗ 4_2 Артём.ipynb @@ -10,7 +10,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "id": "a34b5583", "metadata": {}, "outputs": [ @@ -84,7 +84,7 @@ "4 11.59 1 1" ] }, - "execution_count": 3, + "execution_count": 2, "metadata": {}, "output_type": "execute_result" } @@ -135,7 +135,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "id": "31f5b8b6", "metadata": {}, "outputs": [ @@ -185,7 +185,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "id": "7594c82a", "metadata": {}, "outputs": [ @@ -213,7 +213,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "id": "db397206", "metadata": {}, "outputs": [ @@ -241,7 +241,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "id": "ca70b1e2", "metadata": {}, "outputs": [ @@ -303,7 +303,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "id": "382c3054", "metadata": {}, "outputs": [], @@ -319,7 +319,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "id": "c77e2e2e", "metadata": {}, "outputs": [ @@ -360,7 +360,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 9, "id": "78ffd74b", "metadata": {}, "outputs": [ @@ -386,7 +386,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "id": "87f134c2", "metadata": {}, "outputs": [ @@ -431,7 +431,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 12, "id": "a3830347", "metadata": {}, "outputs": [ @@ -439,22 +439,223 @@ "name": "stdout", "output_type": "stream", "text": [ - "Таблица ANOVA:\n", + "\n", + "Таблица ANOVA для полной модели:\n", " df sum_sq mean_sq F PR(>F)\n", "C(A) 1.0 478.108752 478.108752 631.694471 4.061068e-26\n", "C(B) 3.0 153.241356 51.080452 67.489330 1.051893e-15\n", "C(A):C(B) 3.0 178.558140 59.519380 78.639144 8.022881e-17\n", + "Residual 40.0 30.274683 0.756867 NaN NaN\n", + "\n", + "Отсортированная таблица дисперсионного анализа:\n", + " df sum_sq mean_sq F PR(>F)\n", + "C(A) 1.0 478.108752 478.108752 631.694471 4.061068e-26\n", + "C(A):C(B) 3.0 178.558140 59.519380 78.639144 8.022881e-17\n", + "C(B) 3.0 153.241356 51.080452 67.489330 1.051893e-15\n", "Residual 40.0 30.274683 0.756867 NaN NaN\n" ] } ], "source": [ "from statsmodels.stats.anova import anova_lm\n", + "import pandas as pd\n", + "import statsmodels.api as sm\n", + "from statsmodels.formula.api import ols\n", "\n", - "# ANOVA с взаимодействием\n", - "anova_table = anova_lm(model_full)\n", - "print(\"Таблица ANOVA:\")\n", - "print(anova_table)" + "# Предположим, df - ваш DataFrame с колонками 'Y', 'A', 'B'\n", + "# где A и B - категориальные факторы\n", + "\n", + "# 1. Построение моделей\n", + "model_full = ols('Y ~ C(A) + C(B) + C(A):C(B)', data=df).fit()\n", + "model_additive = ols('Y ~ C(A) + C(B)', data=df).fit()\n", + "model_onlyA = ols('Y ~ C(A)', data=df).fit()\n", + "model_onlyB = ols('Y ~ C(B)', data=df).fit()\n", + "model_const = ols('Y ~ 1', data=df).fit()" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "4d18ccf9", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\margaery\\AppData\\Local\\Temp\\ipykernel_18780\\711296804.py:36: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation.\n", + " return pd.concat([df, pd.DataFrame([row])], ignore_index=True)\n", + "C:\\Users\\margaery\\AppData\\Local\\Temp\\ipykernel_18780\\711296804.py:57: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation.\n", + " AOV = pd.concat([AOV, pd.DataFrame([error_row])], ignore_index=True)\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
HRSSDfMRSSFx_alphaPr(>F)
0H_(12)208.8328233.069.61094178.6391443.6674218.022881e-17
1H_(1)362.0741796.060.34569773.0642372.8765875.459874e-20
2H_(2)686.9415754.0171.735394216.9029763.2953721.537422e-26
3H_(0)840.1829317.0120.026133152.8685562.7448378.275092e-27
4Err30.27468340.00.756867NaNNaNNone
\n", + "
" + ], + "text/plain": [ + " H RSS Df MRSS F x_alpha Pr(>F)\n", + "0 H_(12) 208.832823 3.0 69.610941 78.639144 3.667421 8.022881e-17\n", + "1 H_(1) 362.074179 6.0 60.345697 73.064237 2.876587 5.459874e-20\n", + "2 H_(2) 686.941575 4.0 171.735394 216.902976 3.295372 1.537422e-26\n", + "3 H_(0) 840.182931 7.0 120.026133 152.868556 2.744837 8.275092e-27\n", + "4 Err 30.274683 40.0 0.756867 NaN NaN None" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\n", + "from scipy.stats import f\n", + "\n", + "# Функция для сравнения моделей и извлечения нужных данных\n", + "def compare_models(reduced_model, full_model):\n", + " anova_result = sm.stats.anova_lm(reduced_model, full_model)\n", + " return {\n", + " 'df_diff': anova_result['df_diff'].iloc[1], # Разница степеней свободы\n", + " 'ss_diff': anova_result['ss_diff'].iloc[1], # Разница в RSS\n", + " 'rss_reduced': reduced_model.ssr, # RSS редуцированной модели\n", + " 'rss_full': full_model.ssr, # RSS полной модели\n", + " 'F': anova_result['F'].iloc[1], # F-статистика\n", + " 'p_value': anova_result['Pr(>F)'].iloc[1] # P-значение\n", + " }\n", + "\n", + "# Сравнение моделей\n", + "anova_additive = compare_models(model_additive, model_full)\n", + "anova_A = compare_models(model_onlyA, model_full)\n", + "anova_B = compare_models(model_onlyB, model_full)\n", + "anova_null = compare_models(model_const, model_full)\n", + "\n", + "# Создание DataFrame для результатов\n", + "columns = ['H', 'RSS', 'Df', 'Sum of Sq', 'F', 'x_alpha', 'Pr(>F)']\n", + "AOV = pd.DataFrame(columns=columns)\n", + "\n", + "# Функция для добавления строк в AOV\n", + "def add_aov_row(df, name, anova_result, rdf, alpha):\n", + " row = {\n", + " 'H': name,\n", + " 'RSS': anova_result['rss_reduced'], # Теперь берём RSS редуцированной модели\n", + " 'Df': anova_result['df_diff'],\n", + " 'Sum of Sq': anova_result['ss_diff'],\n", + " 'F': anova_result['F'],\n", + " 'x_alpha': f.ppf(1 - alpha, anova_result['df_diff'], rdf),\n", + " 'Pr(>F)': f\"{anova_result['p_value']:.7g}\", \n", + " }\n", + " return pd.concat([df, pd.DataFrame([row])], ignore_index=True)\n", + "\n", + "# Получение residual degrees of freedom из полной модели\n", + "rdf = model_full.df_resid\n", + "\n", + "# Добавление строк в AOV\n", + "AOV = add_aov_row(AOV, 'H_(12)', anova_additive, rdf, alpha)\n", + "AOV = add_aov_row(AOV, 'H_(1)', anova_A, rdf, alpha)\n", + "AOV = add_aov_row(AOV, 'H_(2)', anova_B, rdf, alpha)\n", + "AOV = add_aov_row(AOV, 'H_(0)', anova_null, rdf, alpha)\n", + "\n", + "# Добавление строки ошибок (RSS полной модели)\n", + "error_row = {\n", + " 'H': 'Err',\n", + " 'RSS': model_full.ssr,\n", + " 'Df': rdf,\n", + " 'Sum of Sq': None,\n", + " 'F': None,\n", + " 'x_alpha': None,\n", + " 'Pr(>F)': None\n", + "}\n", + "AOV = pd.concat([AOV, pd.DataFrame([error_row])], ignore_index=True)\n", + "\n", + "# Вычисление MRSS (Mean Residual Sum of Squares)\n", + "AOV['MRSS'] = AOV['RSS'] / AOV['Df']\n", + "\n", + "# Финальная обработка AOV (перестановка столбцов)\n", + "AOV1 = AOV[['H', 'RSS', 'Df', 'MRSS', 'F', 'x_alpha', 'Pr(>F)']]\n", + "\n", + "# Вывод результатов\n", + "AOV1\n" ] }, { @@ -480,7 +681,7 @@ " $$\n", "\n", "- **Вывод:**\n", - " На уровне значимости $\\alpha=0.02$ все факторы (A, B) и их взаимодействие **значимо** ($p < 0.02$). Это означает, что влияние фактора A на Y зависит от уровня фактора B, и наоборот." + " На уровне значимости $\\alpha=0.02$ все факторы (A, B) и их взаимодействие **значимо**. Это означает, что влияние фактора A на Y зависит от уровня фактора B, и наоборот." ] }, { @@ -493,20 +694,32 @@ "- AIC оценивает баланс между качеством подгонки модели и её сложностью, накладывая штраф за избыточное количество параметров.\n", "- BIC работает аналогично AIC, но применяет более строгий штраф за сложность, особенно при больших объемах данных.\n", "\n", - "Сравниваем две модели:\n", + "Сравниваем модели:\n", "1. **Полная модель** (с взаимодействием): \n", " $$\n", - " Y \\sim A + b + A : B.\n", + " Y \\sim A + B + AB.\n", " $$\n", "2. **Аддитивная модель** (без взаимодействия):\n", " $$\n", " Y \\sim A + B.\n", + " $$\n", + "3. **Только А**:\n", + " $$\n", + " Y \\sim A\n", + " $$\n", + "4. **Только В**:\n", + " $$\n", + " Y \\sim B\n", + " $$\n", + "5. **Константная**:\n", + " $$\n", + " Y \\sim 1\n", " $$" ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 12, "id": "2db6d2ce", "metadata": {}, "outputs": [ @@ -514,23 +727,39 @@ "name": "stdout", "output_type": "stream", "text": [ - "Модель\t\tAIC\tBIC\n", - "Полная \t130.10\t145.07\n", - "Аддитивная \t216.79\t226.15\n" + "\n", + "Сравнение моделей по AIC и BIC:\n", + " Модель AIC BIC\n", + "0 Полная 130.095418 145.065027\n", + "1 Аддитивная 216.794085 226.150090\n", + "2 Только A 237.209208 240.951610\n", + "3 Только B 271.948414 279.433218\n", + "4 Константная 275.614194 277.485395\n" ] } ], "source": [ - "# Список моделей\n", + "\n", + "# 3. Таблица AIC/BIC для всех моделей (как во втором изображении)\n", "models = {\n", - " \"Полная\": model_full,\n", - " \"Аддитивная\": model_additive\n", + " 'Полная': model_full,\n", + " 'Аддитивная': model_additive,\n", + " 'Только A': model_onlyA,\n", + " 'Только B': model_onlyB,\n", + " 'Константная': model_const\n", "}\n", "\n", - "# Вывод AIC и BIC\n", - "print(\"Модель\\t\\tAIC\\tBIC\")\n", + "results = []\n", "for name, model in models.items():\n", - " print(f\"{name} \\t{model.aic:.2f}\\t{model.bic:.2f}\")" + " results.append({\n", + " 'Модель': name,\n", + " 'AIC': model.aic,\n", + " 'BIC': model.bic\n", + " })\n", + "\n", + "aic_bic_table = pd.DataFrame(results)\n", + "print(\"\\nСравнение моделей по AIC и BIC:\")\n", + "print(aic_bic_table)" ] }, { @@ -541,66 +770,11 @@ "### Вывод о сравнении моделей\n", "\n", "- **Результаты AIC и BIC:**\n", - " - **AIC:** Полная модель имеет AIC = 130.10, в то время как аддитивная модель имеет AIC = 216.79. Это указывает на значительное преимущество полной модели.\n", - " - **BIC:** Полная модель также имеет BIC = 145.07, а аддитивная модель — BIC = 226.15. Разница подтверждает выбор полной модели.\n", + " - **AIC:** Минимален у полной модели (не меньше чем на 86.7 меньше, чем у остальных).\n", + " - **BIC:** Минимален у полной модели (не меньше чем на 81.08 меньше, чем у остальных).\n", "\n", "- **Заключение:**\n", - " - Полная модель **предпочтительнее**, так как она лучше соответствует данным, что подтверждается меньшими значениями AIC и BIC.\n", - " - Аддитивная модель не учитывает взаимодействие факторов." - ] - }, - { - "cell_type": "markdown", - "id": "6dff2300", - "metadata": {}, - "source": [ - "## Пункт f)\n", - "### 1. Основные эффекты факторов A и B\n", - "- **Фактор A:** \n", - " Оказал **сильное статистически значимое влияние** на $Y$ ($F=631.69, p<0.001$). \n", - "\n", - "\n", - "- **Фактор B:** \n", - " Также **значимо влияет** на $Y$ ($F=67.49, p<0.001$). \n", - "\n", - "### 2. Взаимодействие факторов $A \\times B$\n", - "- **Статистическая значимость:** \n", - " Взаимодействие **значимо** ($F=78.64, p<0.001$).\n", - " \n", - "- **Визуальное подтверждение:** \n", - " График зависимости $Y$ от $A$ при фиксированных $B$ показывает пересечение линий (особенно для $B=4$), что указывает на **неаддитивность эффектов**.\n", - "\n", - "\n", - "### 3. Выбор оптимальной модели\n", - "- **AIC/BIC:** \n", - " | Модель | AIC | BIC |\n", - " |-----------------|--------|--------|\n", - " | Полная (с взаимодействием) | 130.10 | 145.07 |\n", - " | Аддитивная | 216.79 | 226.15 |\n", - "\n", - " - Разница $\\Delta AIC = 86.69$ и $\\Delta BIC = 81.08$ **явно указывает на преимущество полной модели**. \n", - " - Аддитивная модель не учитывает взаимодействие, что приводит к потере информации.\n", - "\n", - "\n", - "### 4. Нормальность остатков\n", - "- **Тест Шапиро-Уилка:** \n", - " $$p\\text{-value} = 0.949 \\implies \\text{гипотеза о нормальности остатков не отвергается}.$$\n", - "- **Графическая проверка:** \n", - " - Гистограмма остатков близка к нормальной форме. \n", - " - Q-Q график показывает совпадение точек с линией $y = x$.\n", - "\n", - "\n", - "- **Рекомендации:** \n", - " Для прогнозирования $Y$ **необходимо учитывать взаимодействие** $A \\times B$, так как его игнорирование приведет к систематической ошибке.\n", - "\n", - "\n", - "## Итоговый вывод\n", - "1. **Полная модель с взаимодействием** предпочтительна по критериям AIC/BIC и объясняет данные лучше аддитивной. \n", - "2. **Нормальность остатков** подтверждена тестами и графиками.\n", - "\n", - "**Рекомендации:** \n", - "- Проверить данные на наличие выбросов для уровня $B=4$. \n", - "- Использовать полную модель для прогнозирования и анализа эффектов." + " - Полная модель **предпочтительнее**, так как она лучше соответствует данным, что подтверждается меньшими значениями AIC и BIC." ] }, {