cmkl_assignments/fall-2024/math/mat-205/00020/MAT-205_00020.ipynb

286 lines
51 KiB
Plaintext
Raw Normal View History

2024-11-27 20:52:31 +07:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"import scipy.stats as stats\n",
"from scipy.stats import norm, expon\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"DataWhr2024 = pd.read_csv(\"DataWhr2024.csv\")\n",
"UnM49 = pd.read_csv(\"UnM49.csv\", sep=';')\n",
"\n",
"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)\n",
"\n",
"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\"\n",
"\n",
"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\"\n",
"\n",
"_ = 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)\n",
"\n",
"# 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": "markdown",
"metadata": {},
"source": [
"**Question 1**"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(a) P(20 < X ≤ 40): 0.625475\n",
"(b) P(X > 60): 0.000179\n",
"(c) P(X ≤ 30): 0.097141\n",
"(d) P(X = 35): 0.057498\n"
]
}
],
"source": [
"# Define parameters for the normal distribution\n",
"# Assuming p (breast cancer rate) is known; set an example p = 0.0005 (0.05%)\n",
"p = 38/100000\n",
"n = 100000\n",
"mu = n * p # Mean\n",
"sigma = (n * p * (1 - p))**0.5 # Standard deviation\n",
"\n",
"# Define the normal distribution\n",
"normal_dist = norm(loc=mu, scale=sigma)\n",
"\n",
"# (a) P(20 < X ≤ 40)\n",
"prob_a = normal_dist.cdf(40) - normal_dist.cdf(20)\n",
"\n",
"# (b) P(X > 60)\n",
"prob_b = 1 - normal_dist.cdf(60)\n",
"\n",
"# (c) P(X ≤ 30)\n",
"prob_c = normal_dist.cdf(30)\n",
"\n",
"# (d) P(X = 35)\n",
"# For continuous distributions, the probability at an exact point is 0, \n",
"# so use the probability density function (PDF) to approximate it.\n",
"prob_d = normal_dist.pdf(35)\n",
"\n",
"# Print results with 6 decimal places\n",
"print(f\"(a) P(20 < X ≤ 40): {prob_a:.6f}\")\n",
"print(f\"(b) P(X > 60): {prob_b:.6f}\")\n",
"print(f\"(c) P(X ≤ 30): {prob_c:.6f}\")\n",
"print(f\"(d) P(X = 35): {prob_d:.6f}\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 2**"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAHHCAYAAAC88FzIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB7I0lEQVR4nO3dd1xV9R/H8dcFZDgAB4Ib91bcqalZJJaZZuWoHOQot6ENK7fmNjVXuU1LW/prkkZpZeZM0zJz4wInIA5QOL8/rl69ggYInAu8n4/HeXDvuWe8z72MD9/zPd9jMQzDQERERCSHcjI7gIiIiIiZVAyJiIhIjqZiSERERHI0FUMiIiKSo6kYEhERkRxNxZCIiIjkaCqGREREJEdTMSQiIiI5moohERERydFUDEmG8/f3p1u3bmbHkHtYsmQJFouFI0eOmB3FVBn1vXrkyBEsFgtLlixJ923fKbnP0t/fnyeeeCLD9w2wfv16LBYL69evz5T9ZUX6eXM8KoYkVW7+EG/bti3Z1x966CGqVat23/v59ttvGTly5H1vR9LXyJEjsVgs+Pr6cvny5SSvZ+YfXTNZLBbb5OLiQoECBahTpw4DBw7k77//Trf9zJkzJ1MKqLRwxGwPPfSQ7XNxcnLC09OTihUr0rlzZ9atW2d2vHtyxPczJ1ExJBlu3759zJ8/P1XrfPvtt4waNSrDMsn9OX36NHPnzjU7hqkeffRRPvzwQxYvXszo0aOpXbs2S5cupWbNmkybNs1u2VKlSnHlyhU6d+6cqn2k5Q9k586duXLlCqVKlUrVeql1t2xNmzblypUrNG3aNEP3fzfFixfnww8/ZNmyZUyePJknn3yS3377jRYtWtChQweuXbtmSq7bJfcZqRgyl4vZAST7c3NzMztCql26dIk8efKYHcNhBQQEMHnyZPr06YOHh0eG7CMxMZH4+Hjc3d0zZPv3q0KFCrzwwgt28yZMmEDr1q0ZPHgwlSpV4vHHH4cbLUkZfRw3v2ednZ1xdnbO0H3di5OTk6mfmZeXV7Kfy4ABA5gzZw7+/v5MnDjRtHyA6Z+RJKWWIclwd/bDuHbtGqNGjaJ8+fK4u7tTsGBBHnzwQVszdrdu3Zg9ezbccTripkuXLjF48GBKlCiBm5sbFStWZMqUKRiGYbffK1euMGDAAAoVKkS+fPl48sknOXHiBBaLxe4U3M1TP3///TfPPfcc+fPn58EHHwTgzz//pFu3bpQpUwZ3d3f8/Px48cUXOXfunN2+bm7j33//5YUXXsDLywsfHx+GDRuGYRgcO3aMNm3a4OnpiZ+fH1OnTk3Re7d48WIefvhhChcujJubG1WqVEm2Rebm6alff/2V+vXr4+7uTpkyZVi2bFmSZf/66y8efvhhPDw8KF68OGPHjiUxMTFFeW4aPnw4kZGRKWodSunnZbFY6NevHytWrKBq1aq4ubkRGhpqOzX766+/MmDAAHx8fPD29uall14iPj6eqKgounTpQv78+cmfPz+vvfZakm1PmTKFRo0aUbBgQTw8PKhTpw6fffZZqo45JQoWLMjKlStxcXFh3LhxtvnJ9RmKiIggODiY4sWL4+bmRpEiRWjTpo2tH4m/vz9//fUXGzZssP0MPPTQQ3Db6eoNGzbQp08fChcuTPHixe1eS64/ytq1awkICMDd3Z0qVarwxRdf2L1+8/v4Tndu817Z7tZn6NNPP6VOnTp4eHhQqFAhXnjhBU6cOGG3TLdu3cibNy8nTpygbdu25M2bFx8fH4YMGUJCQkIaPhErZ2dnZs6cSZUqVZg1axbR0dF2ry9fvtyWrUCBAnTs2JFjx47ZLXOzC8Dff/9N8+bNyZ07N8WKFWPSpElJ9vfee+9RtWpVcufOTf78+albty4fffRRqt/PQ4cOYbFYePfdd5Ps47fffsNisfDxxx+n+X2RW9QyJGkSHR3N2bNnk8xPSRP0yJEjGT9+PD169KB+/frExMSwbds2duzYwaOPPspLL73EyZMnWbduHR9++KHduoZh8OSTT/LTTz/RvXt3AgIC+P7773n11Vc5ceKE3S+Nbt268cknn9C5c2ceeOABNmzYQKtWre6a69lnn6V8+fK88847tj+m69at49ChQwQHB+Pn58dff/3FBx98wF9//cXvv/+e5A9Hhw4dqFy5MhMmTOCbb75h7NixFChQgPfff5+HH36YiRMnsmLFCoYMGUK9evX+81TC3LlzqVq1Kk8++SQuLi589dVX9OnTh8TERPr27Wu37IEDB3jmmWfo3r07Xbt2ZdGiRXTr1o06depQtWpVuPEHuHnz5ly/fp033niDPHny8MEHH6S6dadJkyY8/PDDTJo0id69e991/dR8XgA//vgjn3zyCf369aNQoUL4+/uzc+dOAPr374+fnx+jRo3i999/54MPPsDb25vffvuNkiVL8s477/Dtt98yefJkqlWrRpcuXWzbnTFjBk8++STPP/888fHxrFy5kmeffZavv/76nt8TaVGyZEmaNWvGTz/9RExMDJ6ensku9/TTT/PXX3/Rv39//P39OX36NOvWrSM8PBx/f3+mT59O//79yZs3L2+99RYAvr6+dtvo06cPPj4+DB8+nEuXLt0z1/79++nQoQMvv/wyXbt2ZfHixTz77LOEhoby6KOPpuoYU5LtdkuWLCE4OJh69eoxfvx4IiMjmTFjBhs3buSPP/7A29vbtmxCQgJBQUE0aNCAKVOm8MMPPzB16lTKli1L7969U5Xzds7OznTq1Ilhw4bx66+/2j73cePGMWzYMNq3b0+PHj04c+YM7733Hk2bNk2S7cKFC7Rs2ZJ27drRvn17PvvsM15//XWqV6/OY489BsD8+fMZMGAAzzzzDAMHDuTq1av8+eefbN68meeeey5V72eZMmVo3LgxK1as4JVXXrFbZ8WKFeTLl482bdqk+T2R2xgiqbB48WIDuOdUtWpVu3VKlSpldO3a1fa8Zs2aRqtWre65n759+xrJfXuuWbPGAIyxY8fazX/mmWcMi8ViHDhwwDAMw9i+fbsBGIMGDbJbrlu3bgZgjBgxwjZvxIgRBmB06tQpyf4uX76cZN7HH39sAMbPP/+cZBu9evWyzbt+/bpRvHhxw2KxGBMmTLDNv3DhguHh4WH3ntxNcvsPCgoyypQpYzevVKlSSTKdPn3acHNzMwYPHmybN2jQIAMwNm/ebLecl5eXARiHDx++Z56bx3nmzBljw4YNBmBMmzbNLsftn21KPy/DWn0aTk5Oxl9//WW37M3vuaCgICMxMdE2v2HDhobFYjFefvll27yb73mzZs3u+T7Gx8cb1apVMx5++OEk72NKPhfA6Nu3711fHzhwoAEYu3btMgzDMA4fPmwAxuLFiw3jxvcAYEyePPme+6latWqSYzFue08efPBB4/r168m+dvtnefP74/PPP7fNi46ONooUKWLUqlXLNu/m53u3/d2+zbtl++mnnwzA+OmnnwzjxntduHBho1q1asaVK1dsy3399dcGYAwfPtw2r2vXrgZgjB492m6btWrVMurUqXPP98owDKNZs2ZJfv/cbvXq1QZgzJgxwzAMwzhy5Ijh7OxsjBs3zm653bt3Gy4uLnbzmzVrZgDGsmXLbPPi4uIMPz8/4+mnn7bNa9OmzT0zGKl8P99//30DMPbu3WubFx8fbxQqVChF36uSMjpNJmkye/Zs1q1bl2SqUaPGf67r7e3NX3/9xf79+1O932+//RZnZ2cGDBhgN3/w4MEYhsF3330HQGhoKNz4z/l2/fv3v+u2X3755STzbm/xuHr1KmfPnuWBBx4AYMeOHUmW79Gjh+2xs7MzdevWxTAMunfvbpvv7e1NxYoVOXTo0H8e7+37v9ka16xZMw4dOpSkqb9KlSo0adLE9tzHxyfJfr799lseeOAB6tevb7fc888//59Z7tS0aVOaN2/OpEmTuHLlSrLLpPTzuqlZs2ZUqVIl2W11797driWuQYMGSd7bm+/5ne/t7e/jhQsXiI6OpkmTJsl+hukhb968AFy8eDHZ1z08PHB1dWX9+vVcuHAhzfvp2bNnivueFC1
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"(a) Probability Life Ladder > 5 (Normal): 0.707072\n",
"(a) Probability Life Ladder > 5 (Empirical): 0.702899\n",
"(b) Probability 4 ≤ Life Ladder ≤ 7 (Normal): 0.809476\n",
"(b) Probability 4 ≤ Life Ladder ≤ 7 (Empirical): 0.833333\n"
]
}
],
"source": [
"Dat1 = Dat2023['Life Ladder']\n",
"\n",
"meanLL = Dat1.mean()\n",
"stdLL = Dat1.std()\n",
"\n",
"normal_dist = norm(loc=meanLL, scale=stdLL)\n",
"\n",
"# Calculate mean and standard deviation of Life Ladder\n",
"meanLL = Dat1.mean()\n",
"stdLL = Dat1.std()\n",
"\n",
"# Define the normal distribution with calculated parameters\n",
"normal_dist = norm(loc=meanLL, scale=stdLL)\n",
"\n",
"# (a) Probability that Life Ladder > 5\n",
"prob_normal_a = 1 - normal_dist.cdf(5) # Using the normal distribution\n",
"prob_empirical_a = (Dat1 > 5).mean() # Empirical probability from dataset\n",
"\n",
"# (b) Probability that Life Ladder is between 4 and 7\n",
"prob_normal_b = normal_dist.cdf(7) - normal_dist.cdf(4) # Using the normal distribution\n",
"prob_empirical_b = ((Dat1 >= 4) & (Dat1 <= 7)).mean() # Empirical probability from dataset\n",
"\n",
"# (c) Histogram and density plot\n",
"plt.hist(Dat1, bins=20, density=True, alpha=0.6, color='skyblue', label='Histogram')\n",
"x_values = np.linspace(min(Dat1), max(Dat1), 1000)\n",
"plt.plot(x_values, normal_dist.pdf(x_values), color='red', label='Normal Density')\n",
"plt.xlabel('Life Ladder')\n",
"plt.ylabel('Density')\n",
"plt.title('Histogram and Normal Distribution Density')\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# Print results\n",
"print(f\"(a) Probability Life Ladder > 5 (Normal): {prob_normal_a:.6f}\")\n",
"print(f\"(a) Probability Life Ladder > 5 (Empirical): {prob_empirical_a:.6f}\")\n",
"print(f\"(b) Probability 4 ≤ Life Ladder ≤ 7 (Normal): {prob_normal_b:.6f}\")\n",
"print(f\"(b) Probability 4 ≤ Life Ladder ≤ 7 (Empirical): {prob_empirical_b:.6f}\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question 3**"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Group 1 (stain=1) - σ: 127.08, λ: 0.01\n",
"Group 2 (stain=2) - σ: 83.72, λ: 0.01\n",
"Probability that Group 1 survives more than 10 years: 0.3889\n",
"Probability that Group 2 survives more than 10 years: 0.2385\n",
"Probability that Group 1 dies within 3 years: 0.2467\n",
"Probability that Group 2 dies within 3 years: 0.3495\n"
]
}
],
"source": [
"# Load the dataset (replace 'PrognosisWomenBreastCancer.csv' with the actual file path)\n",
"DatBreastCancer = pd.read_csv(\"PrognosisWomenBreastCancer.csv\")\n",
"\n",
"# Check the structure of the dataset to find relevant columns\n",
"# Assuming 'time_of_death' is the column representing time of death in months and 'stain' is the group column\n",
"#print(DatBreastCancer.head())\n",
"\n",
"# Separate the data into two groups based on 'stain' (1 or 2)\n",
"group_1 = DatBreastCancer[DatBreastCancer['stain'] == 1]['time']\n",
"group_2 = DatBreastCancer[DatBreastCancer['stain'] == 2]['time']\n",
"\n",
"# Calculate the mean (σ) for each group, which corresponds to the scale parameter of the exponential distribution\n",
"sigma_1 = group_1.mean()\n",
"sigma_2 = group_2.mean()\n",
"\n",
"# Calculate the rate parameter λ (inverse of σ) for each group\n",
"lambda_1 = 1 / sigma_1\n",
"lambda_2 = 1 / sigma_2\n",
"\n",
"# (b) Calculate the probability that each survives more than 10 years (120 months)\n",
"# P(X > 120) = e^(-120 / σ)\n",
"survival_1 = np.exp(-120 / sigma_1)\n",
"survival_2 = np.exp(-120 / sigma_2)\n",
"\n",
"# (c) Calculate the probability that each dies within the next 3 years (36 months)\n",
"# P(X <= 36) = 1 - e^(-36 / σ)\n",
"death_1 = 1 - np.exp(-36 / sigma_1)\n",
"death_2 = 1 - np.exp(-36 / sigma_2)\n",
"\n",
"# Print the results\n",
"print(f\"Group 1 (stain=1) - σ: {sigma_1:.2f}, λ: {lambda_1:.2f}\")\n",
"print(f\"Group 2 (stain=2) - σ: {sigma_2:.2f}, λ: {lambda_2:.2f}\")\n",
"print(f\"Probability that Group 1 survives more than 10 years: {survival_1:.4f}\")\n",
"print(f\"Probability that Group 2 survives more than 10 years: {survival_2:.4f}\")\n",
"print(f\"Probability that Group 1 dies within 3 years: {death_1:.4f}\")\n",
"print(f\"Probability that Group 2 dies within 3 years: {death_2:.4f}\")\n"
]
}
],
"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
}