cmkl/fall-2024/math/mat-206/00010/MAT-206-00010.ipynb

867 lines
237 KiB
Plaintext
Raw Normal View History

2024-11-29 23:54:08 +07:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"import seaborn as sns\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from scipy.stats import norm, t, chi2\n",
"from scipy.stats import iqr, median_abs_deviation"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"# Significane level\n",
"ALPHA = 0.05"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Reading and preprocessing data"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"DataWhr2024 = pd.read_csv(\"DataWhr2024.csv\")\n",
"UnM49 = pd.read_csv(\"UnM49.csv\", sep=';')"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"DataWhr2024.loc[DataWhr2024[\"Country name\"].str.startswith(\"Hong\"), \"Country name\"] = \"Hong Kong\"\n",
"DataWhr2024.loc[DataWhr2024[\"Country name\"].str.startswith(\"Somaliland\"), \"Country name\"] = \"Somaliland\"\n",
"DataWhr2024.loc[DataWhr2024[\"Country name\"].str.startswith(\"Taiwan\"), \"Country name\"] = \"Taiwan\""
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"UnM49 = UnM49[['Country or Area', 'Sub-region Name', 'Region Name']]\n",
"UnM49 = UnM49.rename({'Country or Area':'Country name', 'Sub-region Name':'Subregion', 'Region Name':'Continent'}, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"UnM49.loc[97, \"Country name\"] = \"Bolivia\"\n",
"UnM49.loc[33, \"Country name\"] = \"Congo (Brazzaville)\"\n",
"UnM49.loc[34, \"Country name\"] = \"Congo (Kinshasa)\"\n",
"UnM49.loc[124, \"Country name\"] = \"Hong Kong\"\n",
"UnM49.loc[125, \"Country name\"] = \"Macao\"\n",
"UnM49.loc[126, \"Country name\"] = \"North Korea\"\n",
"UnM49.loc[145, \"Country name\"] = \"Iran\"\n",
"UnM49.loc[46, \"Country name\"] = \"Ivory Coast\"\n",
"UnM49.loc[133, \"Country name\"] = \"Laos\"\n",
"UnM49.loc[129, \"Country name\"] = \"South Korea\"\n",
"UnM49.loc[173, \"Country name\"] = \"Moldova\"\n",
"UnM49.loc[217, \"Country name\"] = \"Netherlands\"\n",
"UnM49.loc[175, \"Country name\"] = \"Russia\"\n",
"UnM49.loc[164, \"Country name\"] = \"Syria\"\n",
"UnM49.loc[26, \"Country name\"] = \"Tanzania\"\n",
"UnM49.loc[116, \"Country name\"] = \"United States\"\n",
"UnM49.loc[193, \"Country name\"] = \"United Kingdom\"\n",
"UnM49.loc[111, \"Country name\"] = \"Venezuela\"\n",
"UnM49.loc[140, \"Country name\"] = \"Vietnam\""
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"_ = pd.DataFrame(\n",
" {\n",
" \"Country name\": [\"Kosovo\", \"Somaliland\", \"Taiwan\"],\n",
" \"Subregion\": [\"Southern Europe\", \"Sub-Saharan Africa\", \"Eastern Asia\"],\n",
" \"Continent\": [\"Europe\", \"Africa\", \"Asia\"],\n",
" }\n",
")\n",
"\n",
"UnM49 = pd.concat([UnM49, _], axis=0)\n",
"UnM49 = UnM49.reset_index(drop=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Merging the datasets"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"# Data\n",
"Dat = pd.merge(DataWhr2024, UnM49)\n",
"\n",
"# Data of 2023\n",
"Dat2023 = Dat[Dat['year'] == 2023]\n",
"Dat2023 = Dat2023.reset_index(drop=True)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"# Extract the 'Life Ladder' data\n",
"Data1 = Dat2023['Life Ladder']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 1**"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Mean (X̄n) Median (X̃n) Std. Dev (σ̂1) MAD (σ̂2) IQR (σ̂3)\n",
"0 5.6208 5.863 1.1395 1.1764 1.3399\n"
]
}
],
"source": [
"# Estimators for μ\n",
"mean_estimator = np.mean(Data1)\n",
"median_estimator = np.median(Data1)\n",
"\n",
"# Estimators for σ\n",
"# Sample standard deviation (1 degree of freedom)\n",
"std_dev_estimator = np.std(Data1, ddof=1)\n",
"\n",
"# Median Absolute Deviation (MAD), scaled\n",
"mad = np.median(np.abs(Data1 - median_estimator))\n",
"mad_estimator = mad * 1.4826\n",
"\n",
"# Interquartile Range (IQR), scaled\n",
"q1 = np.percentile(Data1, 25)\n",
"q3 = np.percentile(Data1, 75)\n",
"iqr = q3 - q1\n",
"iqr_estimator = iqr * 0.7413\n",
"\n",
"# Create a results table\n",
"results = pd.DataFrame({\n",
" 'Mean (X̄n)': [round(mean_estimator, 4)],\n",
" 'Median (X̃n)': [round(median_estimator, 4)],\n",
" 'Std. Dev (σ̂1)': [round(std_dev_estimator, 4)],\n",
" 'MAD (σ̂2)': [round(mad_estimator, 4)],\n",
" 'IQR (σ̂3)': [round(iqr_estimator, 4)]\n",
"})\n",
"\n",
"# Print the results\n",
"print(results)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 2**"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" MLE for μ MLE for σ\n",
"0 5.6208 1.1353\n"
]
}
],
"source": [
"# MLE for μ (mean)\n",
"mle_mu = np.mean(Data1)\n",
"\n",
"# MLE for σ (standard deviation)\n",
"mle_sigma = np.sqrt(np.mean((Data1 - mle_mu) ** 2))\n",
"\n",
"# Report the results rounded to 4 decimal places\n",
"results = pd.DataFrame({\n",
" 'MLE for μ': [round(mle_mu, 4)],\n",
" 'MLE for σ': [round(mle_sigma, 4)]\n",
"})\n",
"\n",
"# Print the results\n",
"print(results)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 3**"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>CI Type</th>\n",
" <th>Lower Bound</th>\n",
" <th>Upper Bound</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>Mean (σ̂1)</td>\n",
" <td>5.4307</td>\n",
" <td>5.8110</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>Mean (σ̂2)</td>\n",
" <td>5.4246</td>\n",
" <td>5.8171</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>Mean (σ̂3)</td>\n",
" <td>5.3973</td>\n",
" <td>5.8444</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>Median (σ̂1)</td>\n",
" <td>5.6247</td>\n",
" <td>6.1013</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>Median (σ̂2)</td>\n",
" <td>5.6170</td>\n",
" <td>6.1090</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>Median (σ̂3)</td>\n",
" <td>5.5828</td>\n",
" <td>6.1432</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" CI Type Lower Bound Upper Bound\n",
"0 Mean (σ̂1) 5.4307 5.8110\n",
"1 Mean (σ̂2) 5.4246 5.8171\n",
"2 Mean (σ̂3) 5.3973 5.8444\n",
"3 Median (σ̂1) 5.6247 6.1013\n",
"4 Median (σ̂2) 5.6170 6.1090\n",
"5 Median (σ̂3) 5.5828 6.1432"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = len(Data1) # Sample size\n",
"Z = 1.96 # Critical value for 95% confidence level\n",
"\n",
"# Mean and median estimators\n",
"mean_estimator = np.mean(Data1)\n",
"median_estimator = np.median(Data1)\n",
"\n",
"# Standard deviation estimators\n",
"std_dev_estimator = np.std(Data1, ddof=1) # σ̂1\n",
"mad = np.median(np.abs(Data1 - median_estimator))\n",
"mad_estimator = mad * 1.4826 # σ̂2\n",
"q1 = np.percentile(Data1, 25)\n",
"q3 = np.percentile(Data1, 75)\n",
"iqr = q3 - q1\n",
"iqr_estimator = iqr * 0.7413 # σ̂3\n",
"\n",
"# CI for the mean\n",
"def CI_mu_mean(hat_mu, hat_sigma, n, alpha=0.05):\n",
" z_critical = norm.ppf(1 - alpha / 2) # Two-tailed Z-critical value\n",
" margin_error = z_critical * (hat_sigma / np.sqrt(n))\n",
" return hat_mu - margin_error, hat_mu + margin_error\n",
"\n",
"# CI for the median\n",
"def CI_mu_median(hat_mu, hat_sigma, n, alpha=0.05):\n",
" z_critical = norm.ppf(1 - alpha / 2)\n",
" density_at_median = 1 / (np.sqrt(2 * np.pi) * hat_sigma) # f( ̃μ)\n",
" std_error = 1 / (4 * n * density_at_median**2)**0.5 # Standard error for the median\n",
" margin_error = z_critical * std_error\n",
" return hat_mu - margin_error, hat_mu + margin_error\n",
"\n",
"ci_mean_sigma1 = CI_mu_mean(mean_estimator, std_dev_estimator, n)\n",
"ci_mean_sigma2 = CI_mu_mean(mean_estimator, mad_estimator, n)\n",
"ci_mean_sigma3 = CI_mu_mean(mean_estimator, iqr_estimator, n)\n",
"ci_median_sigma1 = CI_mu_median(median_estimator, std_dev_estimator, n)\n",
"ci_median_sigma2 = CI_mu_median(median_estimator, mad_estimator, n)\n",
"ci_median_sigma3 = CI_mu_median(median_estimator, iqr_estimator, n)\n",
"\n",
"# Prepare the results in a table, rounded to 4 decimal places\n",
"result_table = pd.DataFrame({\n",
" 'CI Type': ['Mean (σ̂1)', 'Mean (σ̂2)', 'Mean (σ̂3)', \n",
" 'Median (σ̂1)', 'Median (σ̂2)', 'Median (σ̂3)'],\n",
" 'Lower Bound': [round(ci_mean_sigma1[0], 4), round(ci_mean_sigma2[0], 4), round(ci_mean_sigma3[0], 4),\n",
" round(ci_median_sigma1[0], 4), round(ci_median_sigma2[0], 4), round(ci_median_sigma3[0], 4)],\n",
" 'Upper Bound': [round(ci_mean_sigma1[1], 4), round(ci_mean_sigma2[1], 4), round(ci_mean_sigma3[1], 4),\n",
" round(ci_median_sigma1[1], 4), round(ci_median_sigma2[1], 4), round(ci_median_sigma3[1], 4)]\n",
"})\n",
"\n",
"result_table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 4**"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"95% Confidence Interval for σ: (np.float64(1.0191), np.float64(1.2924))\n"
]
}
],
"source": [
"# Degrees of freedom\n",
"df = n - 1\n",
"\n",
"# Critical values for chi-squared distribution\n",
"alpha = 0.05\n",
"chi2_lower = chi2.ppf(alpha / 2, df)\n",
"chi2_upper = chi2.ppf(1 - alpha / 2, df)\n",
"\n",
"# Confidence interval for σ²\n",
"lower_bound_variance = (df * (std_dev_estimator ** 2)) / chi2_upper\n",
"upper_bound_variance = (df * (std_dev_estimator ** 2)) / chi2_lower\n",
"\n",
"# Confidence interval for σ (square root of variance bounds)\n",
"lower_bound_sigma = np.sqrt(lower_bound_variance)\n",
"upper_bound_sigma = np.sqrt(upper_bound_variance)\n",
"\n",
"# Print results rounded to 4 decimals\n",
"ci_sigma = (round(lower_bound_sigma, 4), round(upper_bound_sigma, 4))\n",
"print(f\"95% Confidence Interval for σ: {ci_sigma}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 5**"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>CI Type</th>\n",
" <th>Lower Bound</th>\n",
" <th>Upper Bound</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>Mean (σ̂1)</td>\n",
" <td>3.3794</td>\n",
" <td>7.8623</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>Mean (σ̂2)</td>\n",
" <td>3.3067</td>\n",
" <td>7.9350</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>Mean (σ̂3)</td>\n",
" <td>2.9852</td>\n",
" <td>8.2565</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>Median (σ̂1)</td>\n",
" <td>3.6170</td>\n",
" <td>8.1090</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>Median (σ̂2)</td>\n",
" <td>3.5441</td>\n",
" <td>8.1819</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>Median (σ̂3)</td>\n",
" <td>3.2219</td>\n",
" <td>8.5041</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>IQR</td>\n",
" <td>2.9573</td>\n",
" <td>8.2097</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" CI Type Lower Bound Upper Bound\n",
"0 Mean (σ̂1) 3.3794 7.8623\n",
"1 Mean (σ̂2) 3.3067 7.9350\n",
"2 Mean (σ̂3) 2.9852 8.2565\n",
"3 Median (σ̂1) 3.6170 8.1090\n",
"4 Median (σ̂2) 3.5441 8.1819\n",
"5 Median (σ̂3) 3.2219 8.5041\n",
"6 IQR 2.9573 8.2097"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = len(Data1) # Sample size\n",
"Z = 1.96 # Critical value for 95% confidence level\n",
"\n",
"# Mean and median estimators\n",
"mean_estimator = np.mean(Data1)\n",
"median_estimator = np.median(Data1)\n",
"\n",
"# Standard deviation estimators\n",
"std_dev_estimator = np.std(Data1, ddof=1) # σ̂1\n",
"mad = np.median(np.abs(Data1 - median_estimator))\n",
"mad_estimator = mad * 1.4826 # σ̂2\n",
"q1 = np.percentile(Data1, 25)\n",
"q3 = np.percentile(Data1, 75)\n",
"iqr = q3 - q1\n",
"iqr_estimator = iqr * 0.7413 # σ̂3\n",
"\n",
"\n",
"def PI_mu_mean(hat_mu, hat_sigma, n, alpha=0.05):\n",
" z_critical = norm.ppf(1 - alpha / 2)\n",
" margin_error = z_critical * np.sqrt(hat_sigma**2 + (hat_sigma**2 / n))\n",
" return hat_mu - margin_error, hat_mu + margin_error\n",
"\n",
"def PI_mu_median(hat_mu, hat_sigma, n, alpha=0.05):\n",
" z_critical = norm.ppf(1 - alpha / 2)\n",
" density_at_median = 1 / (np.sqrt(2 * np.pi) * hat_sigma) # f(μ̃)\n",
" std_error = np.sqrt(1 / (4 * n * density_at_median**2) + hat_sigma**2)\n",
" margin_error = z_critical * std_error\n",
" return hat_mu - margin_error, hat_mu + margin_error\n",
"\n",
"def PI_iqr_nonparametric(q1, q3, iqr, alpha=0.05):\n",
" # Critical Z-value\n",
" z_critical = norm.ppf(1 - alpha / 2)\n",
" # Inverse CDF of the normal distribution at 0.75\n",
" phi_inv_3_4 = norm.ppf(0.75)\n",
" # Delta calculation\n",
" delta = 0.5 * ((z_critical / phi_inv_3_4) - 1)\n",
" # Prediction interval\n",
" lower_bound = q1 - delta * iqr\n",
" upper_bound = q3 + delta * iqr\n",
" return lower_bound, upper_bound\n",
"\n",
"pi_mean_sigma1 = PI_mu_mean(mean_estimator, std_dev_estimator, n)\n",
"pi_mean_sigma2 = PI_mu_mean(mean_estimator, mad_estimator, n)\n",
"pi_mean_sigma3 = PI_mu_mean(mean_estimator, iqr_estimator, n)\n",
"pi_median_sigma1 = PI_mu_median(median_estimator, std_dev_estimator, n)\n",
"pi_median_sigma2 = PI_mu_median(median_estimator, mad_estimator, n)\n",
"pi_median_sigma3 = PI_mu_median(median_estimator, iqr_estimator, n)\n",
"\n",
"pi_iqr = PI_iqr_nonparametric(q1, q3, iqr)\n",
"\n",
"# Prepare the results in a table, rounded to 4 decimal places\n",
"result_table = pd.DataFrame({\n",
" 'CI Type': ['Mean (σ̂1)', 'Mean (σ̂2)', 'Mean (σ̂3)', \n",
" 'Median (σ̂1)', 'Median (σ̂2)', 'Median (σ̂3)',\n",
" 'IQR'],\n",
" 'Lower Bound': [round(pi_mean_sigma1[0], 4), round(pi_mean_sigma2[0], 4), round(pi_mean_sigma3[0], 4),\n",
" round(pi_median_sigma1[0], 4), round(pi_median_sigma2[0], 4), round(pi_median_sigma3[0], 4),\n",
" round(pi_iqr[0], 4)],\n",
" 'Upper Bound': [round(pi_mean_sigma1[1], 4), round(pi_mean_sigma2[1], 4), round(pi_mean_sigma3[1], 4),\n",
" round(pi_median_sigma1[1], 4), round(pi_median_sigma2[1], 4), round(pi_median_sigma3[1], 4),\n",
" round(pi_iqr[1], 4)]\n",
"})\n",
"\n",
"result_table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 6**"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAKyCAYAAAAEvm1SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gU5dfG8e+mNxJaQieEAKF3QaT3jqJ0CyACgiiKiq+gFEVR9Idgo1gAxQqKioI0AQERkN5r6C0ESWgJSXbeP5ZdsySB9EnC/bmuvXYy9czszm5y8jznsRiGYSAiIiIiIiIiIpKNXMwOQERERERERERE7j5KSomIiIiIiIiISLZTUkpERERERERERLKdklIiIiIiIiIiIpLtlJQSEREREREREZFsp6SUiIiIiIiIiIhkOyWlREREREREREQk2ykpJSIiIiIiIiIi2U5JKRERERERERERyXZKSomIpFGZMmWwWCzMnj072eURERHUrVsXi8VClSpVOHXqFACrVq3CYrE4Pdzd3SlYsCAVKlSgW7duTJkyhfPnz6d47KNHjybZR0qPo0ePZtk1sLt06RJPPfUUwcHBeHh4YLFYaNas2R23mz17NhaLhTJlyqT6WPbzSs5ff/1FmzZtKFiwIC4uLrd9fTIq8WuQHdf4TnGk5Rra3e5a5mQbNmxgyJAhVKlShfz58+Ph4UFQUBBNmzZlwoQJHD9+3OwQOXXqFI8++ijFixfHzc0Ni8VCv379AOjXr1+63pv2+8W+H7GZO3cu9evXx9fXF39/fypVqsSgQYNYv359uveZ0mesr68vlSpVYtiwYYSHh2fqeWREs2bNsFgsrFq1yml+et9r6ZEb3p+Jv38lbezXLjXf7SIi6eFmdgAiInnJiRMnaN26Nfv376devXosWrSIQoUKJVmvb9++ABiGQXR0NCdOnOCXX37hhx9+YOTIkbz00kuMGTMGd3f3FI/10EMP4efnl+Ly2y3LLIMGDWLevHmUKVOGBx98EC8vLypWrJjlx03s9OnTdOzYkaioKBo1akSZMmVwcXGhXLly2RqHZJ1r167xxBNP8M033wBQtGhRGjVqREBAABcuXGDjxo38+eefvP7663z33Xc88MADpsRpGAYPPvggGzdupHLlyjRv3hx3d3caNWpkSjx52auvvsqECROwWCw0bdqUokWLsmfPHj755BNu3LhBgwYNMnyMxJ+xp06dYsOGDXz00UfMmTOHRYsW0bhx40w4k5zt6NGjhISEEBwcbGoSPifq168fc+bMYdasWTk6ISciktMpKSUikkn2799P69atOXHiBK1atWLBggUpJoaS++/1pUuX+PDDD3nttdeYMGECBw8e5JtvvknxP7vvvvtuulrJZJa4uDgWLFiAl5cX27dvx9/fP0uPt3fv3mTnL126lEuXLtGnTx+++uqrLI1Bsl9cXBxt27Zl7dq1FCtWjOnTp9OlSxendeLj41mwYAGjRo0y9Q/nY8eOsXHjRkqXLs327dtxc3P+NWvixIn83//9H8WKFTMtxrzg5MmTTJw4ETc3N5YsWUKLFi0cy3bv3s2OHTsy5Ti3fsaeOXOGDh06sG3bNvr27cuBAweSvMY5RXa+17p27cq9995LQEBAlh9LRETyHnXfExHJBFu2bKFx48acOHGChx56iN9++y3NLZXy58/PK6+8wo8//ojFYuG7775j7ty5WRZzRp05c4b4+HiKFCmS5QkpgIoVKybbCsveZat8+fJZHoNkv9dff521a9eSP39+1q1blyQhBeDm5kb37t3ZunUrTZs2NSVOEr0XQ0JCkk1WFCtWjIoVK+qP9wz6+++/SUhIoGbNmk4JKYAqVarQu3fvLDlusWLFeO+99wAIDw/nn3/+yZLjZIbsfK8FBARQsWJFJVtFRCRdlJQSEcmgP//8k+bNmxMREcHAgQP5/vvv8fDwSPf+OnXqRLdu3QCYNGlSJkZ6e/v27aN///4EBwfj6elJwYIFadmyJd9//32SdS0WC8HBwXCzdUji2iu31jbJLLfWA7HXMRk7diwA48ePd6xzawuy69ev87///Y97772X/Pnz4+XlRVhYGCNHjiQyMjJL4k1sz549jB07loYNG1KiRAk8PDwoVKgQrVq1Svb6Jvbrr7/StGlT8uXLR0BAAI0bN+bnn3++4zHXr19P+/btyZ8/P35+ftStW5fPP//8jtul9Volridz8eJFnn32WUJDQ/H09MxwDZLLly8zdepUAMaMGUNISMht1/fz86NWrVpJ5i9ZsoROnToRFBSEh4cHxYsXp2fPnikmFRLX6dm2bRsPPvgghQsXxtPTk8qVK/O///0PwzAc69vrENkTYqtXr062vtvt6vzEx8czZcoUqlWrhpeXF4GBgTz00EPs3LnzjtfpwIEDDB48mNDQULy8vAgICKBJkyYpJrXTen63+uOPP+jevTslS5bE09OTwMBA7rnnHsaOHZvseySt8d2JPeF38uRJbty4ka59pFedOnUc0/bXNXHNnWvXrjFmzBgqVaqEj49Pks+izZs38/DDD1O6dGnH52zbtm1ZtGhRisc8ceIEjz/+OMWKFcPLy4vy5cszevRorl+/nuI2d6optXnzZvr27UtISAheXl4ULFiQGjVq8OKLL3Ls2DHHPuz33K2f88l9FqfUhW3jxo306NGD4sWLO+rAde7cmWXLlt0x9vDwcB599FGKFi2Kp6cnoaGhvPLKK8TGxqZ47mllrxN59OhRVq5cSZs2bShQoADe3t7Url2bL774wml9+/0+Z84cAPr37+90XcaNG+e0fmZ/pi5ZsgSLxUKlSpVSPKf4+HiKFi2KxWJh+/btjvkbN25k5MiR1KtXj6JFi+Lh4UGRIkXo3Lkzy5cvT/O127x5Mz179qRkyZJ4eHjg7+9P2bJleeihh1L1PSUiArb6ByIikgbBwcEGYMyaNctYuHCh4eXlZQDGSy+9dNvtVq5caQBGaj56f/75Z8e6Z86cccwPDw93zA8PD8+U8zEMw/j1118d5xEWFmb06tXLaNGiheHq6moAxuOPP+60ft++fY2HHnrIAAxfX1+jb9++jsfevXvveLxZs2YZgBEcHJzqGG+9dmvWrDH69u1r1KhRwwCMGjVqOGJ4/vnnHeudOnXKqFatmgEYBQsWNFq1amV07drV8TqWKVPGOHr0aKrjSM9rMGDAAAMwKlasaLRt29bo2bOn0aBBA8PFxcUAjOeeey7Z7SZPnuw4Vr169YzevXsbdevWNQBjxIgRKV7D77//3vHaVa1a1ejdu7fRqFEjw2KxOLZL7n2Ynmtlfy07duxohISEGAUKFDC6dOlidO/e3Xj44Ycd6yW+b1LLfh9YLBbjwoULqd4usVdeecWxj4YNGxq9e/c2atasaQCGq6ur8dlnnyXZpmnTpgZg/N///Z/h4eFhVKpUyejVq5fRtGlTx3UdPny4Y/2IiAijb9++Rtu2bQ3AKFKkiNM9ERERYRg375vkrkFCQoLxwAMPGIDh4eFhtGnTxujZs6dRpkwZw8vLyxg6dKgBGH379k0S6/fff++4dytWrGh07drVaNGiheHr62sARv/+/TN8fok9/fTTjvdPzZo1jV69ehnt27c3ypYtawDGypUrMxzfnVy8eNGx/dNPP53m7W/nTvf3yZMnHct//PFHw0j02V6/fn3jnnvuMXx9fY327dsbPXv2NFq1auXYdsqUKY57vmbNmka3bt2MRo0aGR4eHgZgjB8/Psnx9u7dawQFBRmAUaxYMaN79+5Ghw4dDG9vb6NBgwZGgwYNkr3uKb3XDMMwJk2a5IijQoUKRo8ePYzOnTsblSpVctrmk08+SfFzPvF70f4ZkNz7c+bMmY5j1apVy+jdu7dx3333Oa7huHHjkmxjj3348OGGv7+/ERwcbPTo0cNo1aqV4e3tbQDGAw88kIpX8z+3+/61fza9+uqrhsViMerUqWP06tXLuPfeex3bvPfee4717fd7aGioARgNGzZ0ui4LFixwrJsVn6kJCQlGyZIlDcBYv359suf7yy+/GIBRu3Ztp/ktW7Y0XFxcjGrVqhk
"text/plain": [
"<Figure size 1200x700 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Calculate basic statistics\n",
"mean_xn = np.mean(Data1)\n",
"std1 = np.std(Data1, ddof=1) # Sample SD\n",
"z_critical = norm.ppf(0.975) # 95% confidence level\n",
"\n",
"# Confidence interval for the mean (μ)\n",
"ci_lower = mean_xn - z_critical * (std1 / np.sqrt(len(Data1)))\n",
"ci_upper = mean_xn + z_critical * (std1 / np.sqrt(len(Data1)))\n",
"\n",
"# Prediction interval for Life Ladder\n",
"pi_lower = mean_xn - z_critical * std1\n",
"pi_upper = mean_xn + z_critical * std1\n",
"\n",
"# Prepare the KDE and prediction distribution\n",
"x_vals = np.linspace(min(Data1), max(Data1), 1000)\n",
"pdf = norm.pdf(x_vals, loc=mean_xn, scale=std1) # Prediction distribution (normal)\n",
"\n",
"# Plot the KDE\n",
"plt.figure(figsize=(12, 7))\n",
"sns.kdeplot(Data1, color=\"blue\", label=\"KDE of Data1\", alpha=0.5)\n",
"\n",
"# Overlay prediction distribution\n",
"plt.plot(x_vals, pdf, color=\"orange\", linestyle=\"--\", label=\"Prediction Distribution\")\n",
"\n",
"# Plot intervals with shaded regions\n",
"plt.hlines(y=0.01, xmin=ci_lower, xmax=ci_upper, color=\"green\", linewidth=2, label=\"Confidence Interval (μ)\")\n",
"plt.hlines(y=0.005, xmin=pi_lower, xmax=pi_upper, color=\"red\", linewidth=2, label=\"Prediction Interval\")\n",
"plt.fill_betweenx([0.01], ci_lower, ci_upper, color=\"green\", alpha=0.2)\n",
"plt.fill_betweenx([0.005], pi_lower, pi_upper, color=\"red\", alpha=0.2)\n",
"\n",
"# Mark mean as a prominent dot for both intervals\n",
"plt.scatter(mean_xn, 0.01, color=\"green\", zorder=5, s=50)\n",
"plt.scatter(mean_xn, 0.005, color=\"red\", zorder=5, s=50)\n",
"\n",
"\n",
"# Add labels, legend, and grid\n",
"plt.title(\"KDE of Life Ladder: Confidence & Prediction Intervals\", fontsize=16)\n",
"plt.ylabel(\"Density\", fontsize=14)\n",
"plt.legend(fontsize=12)\n",
"plt.grid(alpha=0.3)\n",
"\n",
"# Show the plot\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 7**"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"ALPHA = 0.05\n",
"\n",
"N = len(Data1)\n",
"MU = Data1.mean()\n",
"SIGMA = Data1.std()\n",
"\n",
"np.random.seed(111)\n",
"YSample = norm.rvs(MU, SIGMA, size=N)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
"def RelativeLogLikelihood(mu, sigma, data):\n",
" n = len(data)\n",
" sigma_hat = np.std(data)\n",
"\n",
" return n * (np.log(sigma_hat / sigma) + 0.5 * (1 - (np.mean(data**2) - 2 * mu * np.mean(data) + mu**2) / sigma**2))\n",
"\n",
"def RelativeLikelihood(mu, sigma, data):\n",
" return np.exp(RelativeLogLikelihood(mu, sigma, data))"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"def RelativeLogLikelihood(mu, sigma, data):\n",
" n = len(data)\n",
" sigma_hat = np.std(data)\n",
"\n",
" # Vectorized computation of the log-likelihood\n",
" mean_data = np.mean(data)\n",
" mean_squared_data = np.mean(data**2)\n",
" \n",
" return n * (np.log(sigma_hat / sigma) + 0.5 * (1 - (mean_squared_data - 2 * mu * mean_data + mu**2) / sigma**2))\n",
"\n",
"# Define RelativeLikelihood based on the log-likelihood\n",
"def RelativeLikelihood(mu, sigma, data):\n",
" # Vectorized calculation\n",
" return np.exp(RelativeLogLikelihood(mu, sigma, data))\n"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"# Calculate the MLE of mu and sigma\n",
"Mu_MLE = np.mean(YSample)\n",
"Sigma_MLE = np.std(YSample)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"# Find confidence intervals for mu and sigma\n",
"S = np.std(YSample, ddof=1)\n",
"\n",
"LowMu, UppMu = t.ppf(ALPHA/2, N-1, Mu_MLE, S/np.sqrt(N)), t.ppf(1 - ALPHA/2, N-1, Mu_MLE, S/np.sqrt(N))\n",
"LowSigma, UppSigma = S * np.sqrt(N-1) / np.sqrt(chi2.ppf(1-ALPHA/2, N-1)), S * np.sqrt(N-1) / np.sqrt(chi2.ppf(ALPHA/2, N-1))\n",
"\n",
"# Find prediction interval for Y\n",
"LowY, UppY = t.ppf(ALPHA/2, N-1, Mu_MLE, np.sqrt(S**2/N + S**2)), t.ppf(1 - ALPHA/2, N-1, Mu_MLE, np.sqrt(S**2/N + S**2))"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"Probabilities = np.array([0.1, 0.5, 0.75, 0.89, 0.95])\n",
"Probabilities = Probabilities[::-1] # We need to write the probabilities in a dicreasing order\n",
"\n",
"Levels = np.exp(-0.5 * chi2.ppf(Probabilities, 2))"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"mu_vect = np.linspace(5.3, 6, 100) # Extending the range to cover 5.62\n",
"sigma_vect = np.linspace(0.9, 1.3, 100) # Adjusting to cover 1.139\n",
"\n",
"mu_grid, sigma_grid = np.meshgrid(mu_vect, sigma_vect)"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAHZCAYAAAClwGDeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAADpvUlEQVR4nOydd3gU1feH39lNNpteSAgkhARCb6FDkCpNFKQooqAUAQvoF0RRUKRZsAsCKlhA5QdSBRtFkCq99xJIQktCes8m2Z3fH7MJLYEEZnez4b7PM0+S2Tv3ntnMzn7m3HPPkWRZlhEIBAKBQCAow2hsbYBAIBAIBALB3RCCRSAQCAQCQZlHCBaBQCAQCARlHiFYBAKBQCAQlHmEYBEIBAKBQFDmEYJFIBAIBAJBmUcIFoFAIBAIBGUeIVgEAoFAIBCUeYRgEQgEAoFAUOYRgkUgEAgEAkGZRwgWgUAgEAgEZR4hWASC+2Dfvn20adMGV1dXJEni8OHDLFy4EEmSiIqKuuOxU6dORZIkq9laXijp+2sriromBOJ6F9w/DrY2QCC4X86fP88nn3zCP//8w9WrV9HpdDRs2JCnnnqKF154AWdnZ4uMm5eXR//+/dHr9Xz55Ze4uLgQHBwsvqAeYIq7JgQCwf0jBIvArvnrr7/o378/Tk5ODB48mAYNGpCbm8uOHTsYP348J06cYP78+RYZ+/z580RHR/Pdd98xYsSIwv3PPfccTz/9NE5OThYZ90GnLL+/xV0TAoHg/hGCRWC3REZG8vTTTxMcHMy///5L5cqVC18bPXo0ERER/PXXXxYb/9q1awB4eXndtF+r1aLVai02rj2QmZmJq6urRfouy+9vcdfE/WDJ91IgsCdEDIvAbvnkk0/IyMjghx9+uEmsFFCjRg3GjBlT+PehQ4fo0aMHHh4euLm50blzZ3bv3n3TMQXz7BEREQwdOhQvLy88PT0ZNmwYWVlZhe2GDh1Khw4dAOjfvz+SJNGxY0coJsZix44dtGjRAr1eT2hoKPPmzSv2vK5cucLzzz+Pv78/Tk5O1K9fnx9//PGe7Lyxz+HDhxMQEICTkxPVqlXj5ZdfJjc3t1TjFkWBLSdPnmTgwIF4e3vTtm3bUve7ZcsWmjdvftN7VFTcQ1Hvb0n+t6V539LT0xk7diwhISE4OTlRsWJFunbtysGDB4t9H+50TZTUxru9l0VRvXp1nn322dv2d+rUqdCe4oiOjmbUqFHUrl0bZ2dnKlSoQP/+/W+LDyrt9Vaa67045s+fT9OmTXFxcUGSpJu26tWrl7o/gf0jPCwCu+WPP/6gevXqtGnT5q5tT5w4Qbt27fDw8ODNN9/E0dGRefPm0bFjR7Zu3UqrVq1uav/UU09RrVo1ZsyYwcGDB/n++++pWLEiH3/8MQAvvvgigYGBfPjhh/zvf/+jRYsW+Pv7Fzn2sWPH6NatG35+fkydOpX8/HymTJlSZPu4uDhat26NJEm88sor+Pn5sXbtWoYPH05aWhpjx44tlZ0AV69epWXLlqSkpPDCCy9Qp04drly5wooVK8jKykKn05V63KLo378/NWvW5MMPP0SW5VKdz6FDh3jkkUeoXLky06ZNw2g0Mn36dPz8/FT/35bkfXvppZdYsWIFr7zyCvXq1SMxMZEdO3Zw6tQpmjZtWqQdd7omSmtjUe9lUWRkZBAVFcXLL79822tHjx5l4MCBd3zv9u3bx86dO3n66aepUqUKUVFRfPPNN3Ts2JGTJ0/i4uJSqveNUl7vxfHaa68xc+ZMunXrxrBhw7h8+TJffvkleXl59OzZk2bNmpW4L0E5QhYI7JDU1FQZkHv37l2i9n369JF1Op18/vz5wn1Xr16V3d3d5fbt2xfumzJligzIzz///E3H9+3bV65QocJN+zZv3iwD8vLly2/av2DBAhmQIyMjC8fW6/VydHR0YZuTJ0/KWq1WvvUjOHz4cLly5cpyQkLCTfuffvpp2dPTU87Kyiq1nYMHD5Y1Go28b9++294Xk8lUqnGLosCWZ5555rbXStpvr169ZBcXF/nKlSuFbc6dOyc7ODjc9h4V9f6W5H9bmvfN09NTHj16dLHnXBzFXROlvf6Kei+LYteuXTIgr1+//qb9ly5dkgF5/vz5dzy+qP9rQZ8///zzbXaV5HorzfVeFNu2bZMB+eWXX75p/7Rp02RA3rt37137EJRPxJSQwC5JS0sDwN3d/a5tjUYjGzZsoE+fPje5kitXrszAgQPZsWNHYX8FvPTSSzf93a5dOxITE29rV5Kx169fT58+fahatWrh/rp169K9e/eb2sqyzMqVK+nVqxeyLJOQkFC4de/endTU1NumJO5mp8lkYvXq1fTq1YvmzZvfZp8kSfc0blHcaktJ+zUajWzcuJE+ffoQEBBQeHyNGjXo0aPHXd/f0v5vS/K+eXl5sWfPHq5evXrX874balx/xXH8+HEAwsLCbtp/5MgRABo1anTH429cQZeXl0diYiI1atTAy8uryP/53d630lzvxfHll1/i4+PDp59+etP+gumts2fPlqgfQflDCBaBXeLh4QHmWIO7ER8fT1ZWFrVr177ttbp162Iymbh06dJN+2+82QJ4e3sDkJycXCo74+Pjyc7OpmbNmre9dqs98fHxpKSkMH/+fPz8/G7ahg0bBjcEdZbUzvj4eNLS0mjQoMEdbSztuEVRrVq1e+r32rVrZGdnU6NGjdv6LGrfrWOU9n9LCd63Tz75hOPHjxMUFETLli2ZOnUqFy5cuOt7oJaNt76XxXHs2DH8/f1vm245evQoGo3mjv93gOzsbCZPnkxQUBBOTk74+vri5+dHSkoKqampt7UvyfVW0uu9KPLz8/nnn3/o0aPHbYHGBfFWBZ99wYOHiGER2CUeHh4EBAQUPmGqTXGrUO4UT3C/mEwmAJ599lmGDBlSZJtbn5jVsPNexi2KW/PdlLRfo9FYYlvV4m7v21NPPUW7du347bff2LBhA59++ikff/wxq1atuqvXRw1Kmjvo+PHjt3lXAA4fPkz16tXvurro1VdfZcGCBYwdO5bw8HA8PT2RJImnn3668P93I5b+XERFRZGRkVGk0Dpw4ACYRZ7gwUQIFoHd0rNnT+bPn8+uXbsIDw8vtp2fnx8uLi6cOXPmttdOnz6NRqMhKCjIIjb6+fnh7OzMuXPnbnvtVnv8/Pxwd3fHaDTSpUsX1cb38PC4o7CzxLil6ddoNKLX64mIiLjttaL23TqGpf63lStXZtSoUYwaNYpr167RtGlTPvjgg1ILFkvaeOzYMQYMGHDTPpPJxL///kv79u3vevyKFSsYMmQIn3/+eeG+nJwcUlJS7sme0lzvRVHgMdXpdDftl2WZ5cuXU79+/bt63QTlFzElJLBb3nzzTVxdXRkxYgRxcXG3vX7+/HlmzZqFVqulW7durFmz5qblmnFxcSxevJi2bdtazM2s1Wrp3r07q1ev5uLFi4X7T506xfr1629r+8QTT7By5coiBUZ8fHypx9doNPTp04c//viD/fv33/a6LMsWGZdSnI9Wq6VLly6sXr36ppiRiIgI1q5de9cx1P7fGo3G26ZDKlasSEBAAAaDoVR9WcpGzNNp8fHxxMTE3LT/q6++IiEhgYYNG5bItlu9I7Nnz75nr1dprveiKJhy2rhx4037Z86cycGDB5k4ceI92SUoHwgPi8BuCQ0NZfHixQwYMIC6develOl2586dLF++nKFDhwLw/vvv888//9C2bVtGjRqFg4MD8+bNw2Aw8Mknn1jUzmnTprFu3TratWvHqFGjyM/PZ/bs2dSvX5+jR4/e1Pajjz5i8+bNtGrVipEjR1KvXj2SkpI4ePAgGzduJCkpqdTjf/jhh2zYsIEOHTrwwgsvULduXWJiYli+fDk7duzAy8vLIuOW5nymTp3Khg0beOihh3j55ZcxGo3MmTOHBg0a3LXUgdr/2/T0dKpUqcKTTz5JWFgYbm5ubNy4kX379t3kiSgNlrj+jh07BsCGDRsYNWoUderUYffu3YXC4MCBA+zZs6fIZd0F9OzZk19++QVPT0/q1avHrl272Lh
"text/plain": [
"<Figure size 600x500 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(6, 5))\n",
"\n",
"cnt = ax.contour(mu_grid, sigma_grid, RelativeLikelihood(mu_grid, sigma_grid, YSample), Levels)\n",
"ax.scatter(x=Mu_MLE, y=Sigma_MLE, color='k')\n",
"\n",
"ax.plot([LowMu, LowMu, UppMu, UppMu, LowMu], [LowSigma, UppSigma, UppSigma, LowSigma, LowSigma], color='r', ls='--')\n",
"\n",
"ax.set_title(r\"Confidence regions for $\\mu$ and $\\sigma$\")\n",
"ax.set_xlabel(r\"$\\mu$\")\n",
"ax.set_ylabel(r\"$\\sigma$\")\n",
"\n",
"_, labels = cnt.legend_elements()\n",
"ax.legend(_, Probabilities, loc=\"upper right\", frameon=False)\n",
"\n",
"ax.spines[['right', 'top']].set_visible(False)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB86klEQVR4nO3dd1yVdf/H8ddh772RpSLiAic5Sk1ylSMbVpajst9ty27KzEpt2O1Iy4ZpWc670oaVWZneJO5tigMRFBSVIciQPc71+wM4eRQUEbiA83k+HudR55zvua73hQofvtd3aBRFURBCCCGEMCBGagcQQgghhGhoUgAJIYQQwuBIASSEEEIIgyMFkBBCCCEMjhRAQgghhDA4UgAJIYQQwuBIASSEEEIIgyMFkBBCCCEMjhRAVVAUhZycHGSNSCGEEKJ5kgKoCleuXMHe3p4rV66oHUUIIYQQ9UAKICGEEEIYHCmAhBBCCGFwpAASQgghhMGRAkgIIYQQBkcKICGEEEIYHCmAhBBCCGFwpAASQgghhMGRAkgIIYQQBkcKICGEEEIYHCmAhBBCCGFwpAASQgghhMGRAkgIIYQQBkcKICGEEEIYHCmAhBBC3Jbx48czcuRIvdd++OEHLCwsWLBgAePHj0ej0aDRaDA1NcXd3Z177rmHZcuWodVq9T7n7++va3v1Y86cOQ18VaK5kwJICCFEnfryyy8ZM2YMixcv5uWXXwZg8ODBJCcnk5iYyB9//EH//v2ZPHky9913H6WlpXqff+edd0hOTtZ7vPDCCypdjWiuTNQOIERzcurgRfZvOs3l5Fz82rvSZ0RbnDxs2LNnD/b29gQHB6sdUTQhiqJQXFhag5Z1z8zCBI1Gc8ufmzdvHjNnzmTNmjXcf//9utfNzc3x8PAAwNvbmy5dunDHHXcwYMAAVqxYwdNPP61ra2trq2srRH2RAkiIOqAoCr99eYhfPz+oe+3YriQ2rNzFOYv/8eeWXwF48MEHmTFjBh07dlQxrWgqigtLefHO5aqc++PtEzC3NL2lz0ydOpXPPvuMDRs2MGDAgJu2v/vuuwkJCWHdunV6BZAQDUFugQlRB7asPa4rfu4YGsjoV3qSY32K1fvf5M8tv2JkZIRGo+GHH36gX79+FBQUqB1ZiDr1xx9/MG/ePH755ZcaFT+V2rZtS2Jiot5rU6dOxcbGRu+xffv2ekgtDJn0AAlxmy6eyeTHj/cCMOrFMAaNDSm/BbBtAQCOFl70bzmeJ98KZ8V3nxESEoKlpaXKqUVTYGZhwsfbJ6h27lvRqVMn0tPTmTlzJj169MDGxqZGn1MU5bpbbVOmTGH8+PF6r3l7e99SHiFuRgogIW7T9x/sprS4jA69fBj4RCcAXn/9dSwtLcnPz8erqBd/R55j+4qLrPhmNVZ25rrPVvXNX4hKGo3mlm9DqcXb25sffviB/v37M3jwYP744w9sbW1v+rmYmBgCAgL0XnNxcaF169b1mFYIuQUmxG1JOJbGiT3nMTLW8OjU3rpixtzcnKlTp/LOO+8wfsbduPnak5max3cLduvabNiwgdDQUC5cuKDyVQhRN/z8/Ni6dSspKSkMHjyYK1eu3LD9X3/9xdGjR3nggQcaLKMQlaQAEuI2/LHsbwDChgTi4m1HTEyMbkpvZaFjYW3Gk+/0R6OBPb/HEXvgIoqi8N577xEdHa2bJixEc+Dj40NUVBRpaWkMGjSInJwcAIqKikhJSeHChQscOnSI//znP4wYMYL77ruPsWPH6h3jypUrpKSk6D0qjyNEXZECSIhaykzLI3r7WQAGjw8lNzeXu+++m44dO3L69Gm9tgEd3LjrwXYAfDtvJ4pWYfHixRgZGbF27Vo2b96syjUIUR9atGhBVFQU6enpuiJo48aNeHp64u/vz+DBg9myZQsff/wxv/zyC8bGxnqfnzFjBp6ennqPV199VbXrEc2TjAESopb2/RGHokDrzh54+Dswc+ZMUlJSsLKyokWLFte1H/lsdw5sOk3ymUz2/B5Hr2GhPP/883z88cc899xzHD16FHNz8yrPJURjtmLFiute8/b25tSpU7d8rGtnhAlRX6QHSIhaUBSFPb/FAdDz3jZcuHCB999/H4C5c+dWWchY2ZozeHwoAL9+fpCS4jLeeecdPDw8iIuLY8mSJQ18FUIIYbikABKiFi7EX+bimUxMzIzpGt6Sd999l4KCAnr37n3DAZ39HmqPg5s1l1Ny2fbjCezt7XnzzTcBWL5cnQXvhBDCEEkBJEQtRG8rH/vT/o4WaEy0fPvtt1Cxh9GNprWbWZhw38QuAPz+1d8U5hXz6KOPYmZmxpEjR4iOjm6gKxBCCMPWKAqgRYsW4e/vj4WFBWFhYezbt6/atuvWraNbt244ODhgbW1NaGgoq1ev1mtz9c7DlY/Bgwc3wJUIQxG94xwAHe/05Y8//iAnJ4cWLVrQr1+/m36217Ag3H3tyc0qZMt3x3FycmL+/Pls3ryZ9u3bN0B6IYQQqhdAa9euJSIigpkzZ3Lo0CFCQkIYNGgQaWlpVbZ3cnLijTfeYPfu3URHRzNhwgQmTJjAn3/+qdeucufhykflb+hC3K6cywUkHiv/+9mhty9RUVEAPPLIIxgZ3fyflLGJEUOfLu8F+mvNMYoLS3nhhRcIDw+/bjaMEEKI+qF6AfTBBx8wceJEJkyYQLt27ViyZAlWVlYsW7asyvb9+vXj/vvvJzg4mFatWjF58mQ6derEjh079NpV7jxc+XB0dGygKxLN3fFdSSgK+AQ54+hmzcKFC/n77795/vnna3yM7gNb4eRhQ05GAbs33PpMGSGEELdH1QKouLiYgwcPEh4e/k8gIyPCw8PZvXv3TT+vKAqRkZHExsZy11136b0XFRWFm5sbQUFBTJo0iYyMjHq5BmF4YvaeB6BDLx+oWPAwNDQUPz+/Gh/D2MSIex4v3zZj83+j0ZZpOXv2LK+88gqTJ0+up+RCCCEqqboOUHp6OmVlZbi7u+u97u7uzsmTJ6v9XHZ2Nt7e3hQVFWFsbMxnn33GPffco3t/8ODBjBo1ioCAAE6fPs3rr7/OkCFD2L17d5W3GIqKiigqKtI9lxVHRXUUReHUoWQAgrp7U1JSgqlp7fZq6j0iiA1LD3LpfA6H/koAx0wWLFiAhYUF7777LnZ2dnWcXgghRCXVb4HVhq2tLYcPH2b//v289957RERE6MZhUDEWY/jw4XTs2JGRI0eyYcMG9u/fr9fmarNnz8be3l738PHxacCrEU1J+oUrZKbmYWxihI07uLm5MXbsWIqLi2/5WOaWptz9SAcANq44TJcuXQgODqawsJAffvihHtIL0fSNHz+ekSNH6p7369ePl1566baOWRfHqI2oqCg0Gg1ZWVlQsaCkg4NDnZ8nMTERjUbD4cOHqzxvfZ6rMVO1AHJxccHY2JjU1FS911NTU/Hw8Kj2c0ZGRrRu3ZrQ0FBefvllHnzwQWbPnl1t+5YtW+Li4kJ8fHyV70+bNo3s7GzdIykp6TauSjRnpw5ehIqtLX759SeysrKIiYnBzMysVsfr91B7TM2NSYrN4Ex0mm5PpFWrVtVpbiHq09Uzb83MzGjdujXvvPOObl+8+rRu3TrefffdGrWt7gf/rRyjPo0ePbrGq2ffSrHk4+NDcnIyHTp0uM2E+q4tRuvzXPVB1QLIzMyMrl27EhkZqXtNq9USGRlJz549a3wcrVardwvrWufPnycjIwNPT88q3zc3N8fOzk7vIURVTh0sv/3VpqunbmbhY489Vuvj2ThYEDY0EIAta4/x+OOPo9Fo2Lp1q2wJIJqUypm3cXFxvPzyy7z11lu61dGvVZse0+o4OTlha2ur+jHqgqWlJW5ubnV6zOLiYoyNjfHw8MDEpP5HvTTkuW6X6rfAIiIiWLp0KStXriQmJoZJkyaRl5fHhAkTABg7dizTpk3TtZ89ezabN2/mzJkzxMTEsGDBAlavXs3jjz8OQG5uLlOmTGHPnj0kJiYSGRnJiBEjaN26NYMGDVLtOkXzEH8kBQBL9xJ27dqFRqNh9OjRt3XMfg+Vr/1z6K8EbMw
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"y_vect = np.linspace(4, 16, 100)\n",
"\n",
"ax = sns.kdeplot(YSample, color=\"rebeccapurple\", label=\"KDE\")\n",
"ax.plot(y_vect, t.pdf(y_vect, N-1, Mu_MLE, np.sqrt(S**2/N + S**2)), color=\"k\", ls=\"--\", label=\"Prediction distribution\")\n",
"\n",
"ax.hlines(0.025, LowY, UppY, color='red', label=\"Prediction interval\")\n",
"ax.hlines(0.05, LowMu, UppMu, color='darkorange', label=\"Confidence interval\")\n",
"ax.scatter(Mu_MLE, 0.025, color='red')\n",
"ax.scatter(Mu_MLE, 0.05, color='darkorange')\n",
"\n",
"ax.legend(frameon=False)\n",
"ax.spines[['right', 'top']].set_visible(False)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}