univer-6th-semester/matstat
Code Actions Releases 2

Некоторые исправления (спасибо моей Ритке :)

Browse Source
This commit is contained in:
Артём Тищенко
2025-06-07 11:23:04 +03:00
parent 2064ae46b0
commit 85fc09ad29
3 changed files with 718 additions and 300 deletions
Download Patch File Download Diff File Expand all files Collapse all files

315
idz3/ИДЗ 3_2 Артём.ipynb
View File

@@ -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}\")"
]
},
{