2024-10-16 10:21:36 +07:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
2024-11-19 23:03:48 +07:00
"outputs": [
{
"ename": "ModuleNotFoundError",
"evalue": "No module named 'numpy'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[1], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;66;03m#Numerical Python\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;66;03m# To plot\u001b[39;00m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyplot\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mplt\u001b[39;00m\n",
"\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'numpy'"
]
}
],
2024-10-16 10:21:36 +07:00
"source": [
"#Numerical Python\n",
"import numpy as np\n",
"\n",
"# To plot\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"\n",
"# Gamma function\n",
"from scipy.special import gamma\n",
"\n",
"# To calculate statistics\n",
"from statistics import mode\n",
"from scipy.stats import norm\n",
"from scipy.stats import hmean, trim_mean, iqr, median_abs_deviation, skew, kurtosis\n",
"from scipy.stats.mstats import gmean, winsorize\n",
"from statsmodels.distributions import ECDF"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"ALPHA = 0.2"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"colors = plt.rcParams['axes.prop_cycle'].by_key()['color']"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"Dat = np.array([0.0, 0.8, 1.0, 1.2, 1.3, 1.3, 1.4, 1.8, 2.4, 4.6])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Empirical cumulative distribution function (ECDF)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"ecdf = ECDF(Dat)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-0.05, 4.6499999999999995)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1EAAAINCAYAAADfvvWSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABHhElEQVR4nO3dfVxUdf7//+eAMigKYlx5QYJGmdeCStKapii0ZrldYG6lkrVWympollbatWWuq1um1me9aPuWqK1t24XGQmqprRdoF5aXYXYFyKqDgoLOnN8f/ZwcuZCDwAzwuN9uc7s573mfc17ncGY8z3mfc8ZiGIYhAAAAAECleLm7AAAAAACoSwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAA9cayZctksVicD19fX7Vu3VoJCQn629/+phMnTlRpvps3b9aTTz6p48ePV2/BAIA6iRAFAKh3nn76af3jH//QwoULlZKSIkmaNGmSunbtqi+//NL0/DZv3qynnnqKEAUAkCQ1cncBAABUtxtuuEG9evVyPp82bZoyMzN144036qabbtK3336rJk2auLFCAEBdxkgUAKBBGDhwoJ544gl9//33evPNNyVJX375pcaMGaP27dvL19dXYWFhuueee/S///3POd2TTz6phx9+WJIUGRnpPFXw0KFDkqSlS5dq4MCBCgkJkdVqVadOnbRw4cJaXz8AQO1hJAoA0GDcfffdmj59uj7++GPdd999Sk9P13fffafk5GSFhYVp9+7deu2117R79259/vnnslgsuuWWW7Rv3z69/fbb+utf/6qgoCBJUnBwsCRp4cKF6ty5s2666SY1atRI//73v/Xggw/K4XBo/Pjx7lxdAEANsRiGYbi7CAAAqsOyZcuUnJysbdu2uZzOd74WLVqoffv2ysrK0qlTp0qd1rdixQqNHDlSGzduVL9+/SRJc+bM0cMPP6zs7GxFRES49C9rHomJidq/f78OHjxYfSsHAPAYnM4HAGhQmjVr5rxL3/nh5/Tp08rPz9c111wjScrKyqrU/M6fh81mU35+vvr376/vvvtONputGisHAHgKQhQAoEE5efKkmjdvLkk6evSoJk6cqNDQUDVp0kTBwcGKjIyUpEoHoE2bNik+Pl5+fn5q0aKFgoODNX36dFPzAADULVwTBQBoMH788UfZbDZdccUVkqSkpCRt3rxZDz/8sHr06KFmzZrJ4XAoMTFRDofjovM7ePCgBg0apI4dO2ru3LkKDw+Xj4+PPvzwQ/31r3+t1DwAAHUPIQoA0GD84x//kCQlJCTo2LFjysjI0FNPPaUZM2Y4++zfv7/UdBaLpcz5/fvf/1ZxcbHee+89XX755c72Tz75pJorBwB4Ek7nAwA0CJmZmXrmmWcUGRmpO++8U97e3pKkC++vNG/evFLT+vn5SVKpH9stax42m01Lly6txsoBAJ6GkSgAQL3z0Ucfac+ePTp79qxyc3OVmZmp9PR0tWvXTu+99558fX3l6+ur6667TrNnz9aZM2fUpk0bffzxx8rOzi41v5iYGEnSY489pjvuuEONGzfWsGHDNGTIEPn4+GjYsGEaN26cTp48qddff10hISH65Zdfanu1AQC1hBAFAKh3zp2e5+Pjo5YtW6pr166aN2+ekpOTnTeVkKS33npLKSkpWrBggQzD0JAhQ/TRRx+pdevWLvPr3bu3nnnmGS1atEhr166Vw+FQdna2rrrqKq1evVqPP/64pkyZorCwMD3wwAMKDg7WPffcU6vrDACoPfxOFAAAAACYwDVRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACY0MjdBXg0w5CKin79d9OmksXi3nrgXp6wP9SHGhwOKT//138HBUleVfguxxO2gydgO3gO/hYA4Nmq+XOakaiKFBVJzZr9+ji30dFwecL+UB9qyM+XQkN/fZwLU7VdQ33BdvAc/C0AwLNV8+c0IQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACY0cncB7mYYhk6cOFH2i4WFv/27oECy22unKHgmT9gf6kMN57/fTpyQfH1rv4b6gu3gOfhbAIBnM/k53bx5c1kslnJftxiGYVRXbXVRQUGBAgIC3F0GAAAAAA9hs9nk7+9f7usNPkRVOBKlX0NWeHi4fvjhhwo3JFBb2Cfhadgn4WnYJ+FJ2B/rpouNRDX40/ksFkuldmh/f392fHgU9kl4GvZJeBr2SXgS9sf6hRtLAAAAAIAJhCgAAAAAMIEQdRFWq1UzZ86U1Wp1dymAJPZJeB72SXga9kl4EvbH+qnB31gCAAAAAMxgJAoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhSAem39+vWyWCxav359pfoPGDBAAwYMqLF6IiIiNGbMmBqbv6cwu90ry2Kx6Mknn6zWedaEQ4cOyWKxaNmyZTW+rGXLlslisejQoUPOtoiICN144401vmyp5v7W9dWTTz4pi8Xi7jIAXCJCFIBqc+5grrzH559/7u4S65TTp0/rr3/9q2JjYxUQECBfX19deeWVmjBhgvbt2+fu8mrMhx9+6HFB6fz9uFGjRmrZsqViYmI0ceJEffPNN9W2nFdffbVWgldVeGJtAwYMKPfzZs+ePW6rq6ioSE8++STBEqjHLIZhGO4uAkD9sGzZMiUnJ+vpp59WZGRkqdcTExMVFBRUqzU5HA6VlJTIx8dHXl4X/96opKREkuTj41Mj9URERGjAgAEXPRjNz89XYmKiduzYoRtvvFHx8fFq1qyZ9u7dqxUrVignJ8dZqydav369rr/+en3yySemR/YmTJigBQsWqKz/nk6fPq1GjRqpUaNG1VRp5VgsFg0ePFijRo2SYRiy2Wz64osvtGrVKhUWFurFF19Uamqqs79hGCouLlbjxo3l7e1d6eV06dJFQUFBpg6+7Xa7zpw5I6vV6hzhiIiIUJcuXfT+++9Xej5Vrc3se6w6DRgwQAcPHtSsWbNKvXbTTTfJ39+/Vus5Jz8/X8HBwZo5c2apLwTOnj2rs2fPytfX1y21Aagetfu/EIAG4YYbblCvXr3cXYYkycvLq1IHK0VFRWratGmNhSezxowZo507d2r
"text/plain": [
"<Figure size 1000x600 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Divide the figure in two subfigures\n",
"f, ax = plt.subplots(2, 1, figsize=(10,6), sharex=True, height_ratios=[1, 3]) # The subfigures share the same x axis\n",
"\n",
"# In the first subfigure we show our observations\n",
"ax[0].set_yticks([]) # Do not show the ticks in the vertical axis\n",
"ax[0].set_ylabel(\"\") # Do not show any label in the vertical axis\n",
"ax[0].spines[['left', 'right', 'top']].set_visible(False) # Do not plot the left, right and top margins\n",
"ax[0].text(2.2, 0.75, \"Data\", fontsize=12)\n",
"\n",
"for i in range(len(Dat)):\n",
" ax[0].axvline(x=Dat[i], ymax=0.2, c='r') # Plot the observations as a rugplot\n",
"\n",
"# In the second subfigure we show the ECDF\n",
"ax[1].spines[['right', 'top']].set_visible(False) # Do not plot the right and top margins\n",
"\n",
"# Plot the ECDF\n",
"ax[1].scatter(ecdf.x, ecdf.y, color='rebeccapurple') \n",
"ax[1].step(ecdf.x, ecdf.y, where='post', color='rebeccapurple')\n",
"\n",
"# Add title and labels to the plot\n",
"ax[1].set_title(\"Empirical Cumulative Distribution Function\", fontsize=12)\n",
"ax[1].set_xlabel(r'$x$')\n",
"ax[1].set_ylabel(r'$F_n(x)$')\n",
"\n",
"# Modify the limits of the vertical and horizontal axes\n",
"ax[1].set_ylim(0, 1.05)\n",
"ax[1].set_xlim(Dat.min()-0.05, Dat.max()+0.05)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Quantiles, quartiles, percentiles"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.792 , 1.3689085, 2.422 ])"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Quantiles of specific probabilities\n",
"p = [0.11, 1-np.exp(-1), 0.89]\n",
"np.quantile(Dat, p)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1.05, 1.3 , 1.7 ])"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Quartiles\n",
"p = [0.25, 0.5, 0.75]\n",
"np.quantile(Dat, p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Central tendency statistics"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"Mode = mode(Dat)\n",
"Median = np.median(Dat)\n",
"Mean = Dat.mean()\n",
"GeomMean = gmean(Dat)\n",
"HarmMean = hmean(Dat)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"Weights = norm.pdf(Dat, Dat.mean(), Dat.std()) # Weights for the weighted mean"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGxCAYAAACKvAkXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwPUlEQVR4nO3df1RVdb7/8deBhnP8iTIEB4wCf4wOo0iK8CUzf5HgtFx6b9NS1zQaq2szJN38nltNNAqazmDmeKmrA42z7IdleqcpZ9k4mHMSyxuJQU6Z1soGr5ocQB1BcYTinO8ffj15BjBAOPvAfj7W2ivO53z257x3p+TlZ3/23haPx+MRAACAiQQZXQAAAIC/EYAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIBa4fF4VF9fL+4RCQBA70QAasX58+cVGhqq8+fPG10KAADoBgERgDZs2KDY2FjZbDalpKSorKyszb6vv/66kpKSNGjQIPXr10+JiYnavHmzT5/77rtPFovFZ8vIyOjuwwAAAD3EDUYXsG3bNjkcDhUVFSklJUUFBQVKT0/XZ599poiIiBb9w8LC9Itf/EKjRo1SSEiI3nzzTWVmZioiIkLp6enefhkZGXr++ee9r61Wq1+OBwAABD6L0Q9DTUlJ0YQJE7R+/XpJktvtVkxMjB566CE9/vjj7Rpj3Lhxuuuuu7Ry5UpJl2eAzp07p+3bt3eqpvr6eoWGhqqurk4DBw7s1BgAACBwGXoKrKmpSeXl5UpLS/O2BQUFKS0tTaWlpd+6v8fjkdPp1GeffaY77rjD572SkhJFRERo5MiRysrK0pkzZ9ocp7GxUfX19T4bAADovQw9BXb69Gk1NzcrMjLSpz0yMlKffvppm/vV1dVpyJAhamxsVHBwsH7zm9/ozjvv9L6fkZGhf/3Xf1VcXJy++OILPfHEE5o5c6ZKS0sVHBzcYrz8/HytWLGi6w4MAAAENMPXAHXGgAEDdPDgQV24cEFOp1MOh0NDhw7VlClTJEnz5s3z9h0zZowSEhI0bNgwlZSUaPr06S3Gy8nJkcPh8L6ur69XTExMtx8HAAAwhqEBKDw8XMHBwaqurvZpr66ult1ub3O/oKAgDR8+XJKUmJioI0eOKD8/3xuA/tnQoUMVHh6uo0ePthqArFYri6QBADARQ9cAhYSEaPz48XI6nd42t9stp9Op1NTUdo/jdrvV2NjY5vsnT57UmTNnFBUVdV31AgCA3sHwU2AOh0MLFy5UUlKSkpOTVVBQoIaGBmVmZkqSFixYoCFDhig/P1/S5fU6SUlJGjZsmBobG7Vz505t3rxZhYWFkqQLFy5oxYoVuvvuu2W32/XFF1/oscce0/Dhw30ukwcAAOZleACaO3euamtrlZubK5fLpcTERBUXF3sXRh8/flxBQd9MVDU0NOjBBx/UyZMn1adPH40aNUovv/yy5s6dK0kKDg7WRx99pBdffFHnzp1TdHS0ZsyYoZUrV3KaCwAASAqA+wAFIu4DBABA72b4DBCAbzS7PSqrPKua85cUMcCm5LgwBQdZjC4LAHodAhAQIIoPVWnFjsOqqrvkbYsKtSlvVrwyRrOAHwC6UkA8DBUwu+JDVcp6ucIn/EiSq+6Ssl6uUPGhKoMq+0az26PSL87ojwe/VOkXZ9Ts5uw5gJ6LGSDAYM1uj1bsOKzW4oRHkkXSih2HdWe83bDTYcxOAehtmAECDFZWebbFzM/VPJKq6i6prPKs/4q6Sk+YnQKAjiIAAQarOd92+OlMv670bbNT0uXZKU6HAehpOAUG+EHT125tLj2m/z17UbeE9dVPUmMVcsPlv39EDLC1a4z29utKHZmdSh32Xf8VBgDXiQAEdLP8nYe18d1KXT1J8sudR7RoUpxyfhiv5LgwRYXa5Kq71OpMi0WSPfTyJfH+FsizUwBwPTgFBnSj/J2H9dw7vuFHktwe6bl3KpW/87CCgyzKmxUv6XLYudqV13mz4g1ZAB3Is1MAcD0IQEA3afrarY3vVl6zz8Z3K9X0tVsZo6NUeO84RQz0fVyLPdSmwnvHGXal1ZXZqbail0WXrwYzYnYKAK4HAQjoJptLj7WY+flnbs/lfpKUMTpKf3FM9r73QuYE7fv5NEMvMw/k2SkAuB4EIKCb/O/Zix3ud3WQCJTHYATq7BQAXA8WQQPd5Jawvl3az0gZo6M0cXi4xix/S9Ll2alJI24MiIAGAJ3BDBBwHa71eIifpMbq2/JBkOVyv54gEGenAKCzmAECOunbHg8RckOQFk2K03PvtL0QetGkOO/9gAAA/sOfvEAntPfxEDk/jNdP74hrMRMUZJF+esfl+wABAPyPAAR0UEcfD5Hzw3jt+/lUb5//mzZCn66cSfgBAAMRgIAO6szDS68+zfXj/3Mzp70AwGD8KQx0EI+HAICejwAEdBCPhwCAno8ABHQQj4cAgJ6PAAR0EI+HAICejwAEdAKPhwCAno0bIQKdxOMhAKDnYgYIuA48HgIAeiYCEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB2eBQbTaXZ7VFZ5VjXnLyligI1HWACACRGAYCrFh6q0YsdhVdVd8rZFhdqUNyueJ7gDgIlwCgymUXyoSlkvV/iEH0ly1V1S1ssVKj5UZVBlAAB/IwDBFJrdHq3YcVieVt670rZix2E1u1vrAQDobQIiAG3YsEGxsbGy2WxKSUlRWVlZm31ff/11JSUladCgQerXr58SExO1efNmnz4ej0e5ubmKiopSnz59lJaWps8//7y7DwMBrKzybIuZn6t5JFXVXVJZ5Vn/FQUAMIzhAWjbtm1yOBzKy8tTRUWFxo4dq/T0dNXU1LTaPywsTL/4xS9UWlqqjz76SJmZmcrMzNSuXbu8fdasWaNnn31WRUVF2r9/v/r166f09HRdutT2L0D0bjXn2/fdt7cfAKBnMzwArVu3TosWLVJmZqbi4+NVVFSkvn37atOmTa32nzJliv7lX/5F3//+9zVs2DA9/PDDSkhI0L59+yRdnv0pKCjQ0qVLNXv2bCUkJOill17SqVOntH37dj8eGQJJxABbl/YDAPRshgagpqYmlZeXKy0tzdsWFBSktLQ0lZaWfuv+Ho9HTqdTn332me644w5JUmVlpVwul8+YoaGhSklJaXPMxsZG1dfX+2zoXZLjwhQValNbF7tbdPlqsOS4MH+WBQAwiKEB6PTp02publZkZKRPe2RkpFwuV5v71dXVqX///goJCdFdd92l//qv/9Kdd94pSd79OjJmfn6+QkNDvVtMTMz1HBYCUHCQRXmz4iWpRQi68jpvVjz3AwIAkzD8FFhnDBgwQAcPHtSBAwf0y1/+Ug6HQyUlJZ0eLycnR3V1dd7txIkTXVcsAkbG6CgV3jtOEQOtPu32UJsK7x3HfYAAwEQMvRFieHi4goODVV1d7dNeXV0tu93e5n5BQUEaPny4JCkxMVFHjhxRfn6+pkyZ4t2vurpaUVHf/EKrrq5WYmJiq+NZrVZZrdZW30PvkjE6ShOHh2vM8rckSS9kTtCkETcy8wMAJmPoDFBISIjGjx8vp9PpbXO73XI6nUpNTW33OG63W42NjZKkuLg42e12nzHr6+u1f//+Do2J3uvqsMNjMADAnAx/FIbD4dDChQuVlJSk5ORkFRQUqKGhQZmZmZKkBQsWaMiQIcrPz5d0eb1OUlKShg0bpsbGRu3cuVObN29WYWGhJMlisWjJkiVatWqVRowYobi4OC1btkzR0dGaM2eOUYcJAAACiOEBaO7cuaqtrVVubq5cLpcSExNVXFzsXcR8/PhxBQV9M1HV0NCgBx98UCdPnlSfPn00atQovfzyy5o7d663z2OPPaaGhgY98MADOnfunG6//XYVFxfLZuMSZwAAEAABSJKys7OVnZ3d6nv/vLh51apVWrVq1TXHs1gsevLJJ/Xkk092VYkAAKAX6ZFXgQEAAFwPAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADA
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plot the weight for each observation\n",
"ax = plt.subplot(111)\n",
"\n",
"ax.stem(Dat, Weights)\n",
"ax.set_xlabel('x')\n",
"ax.set_ylabel('Weight')\n",
"\n",
"ax.spines[['right', 'top']].set_visible(False)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"WeightedMean = np.average(Dat, weights=Weights)\n",
"TrimmedMean = trim_mean(Dat, 0.1)\n",
"InterquartileMean = trim_mean(Dat, 0.25)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"masked_array(data=[0.8, 0.8, 1. , 1.2, 1.3, 1.3, 1.4, 1.8, 2.4, 2.4],\n",
" mask=False,\n",
" fill_value=1e+20)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"DatWinsorized = winsorize(Dat, [0.1,0.1]) # Winsorized the data\n",
"DatWinsorized"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"WinsorizedMean = DatWinsorized.mean()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"CentralTendency = [Mode, Median, Mean, GeomMean, HarmMean, WeightedMean, TrimmedMean, InterquartileMean, WinsorizedMean]"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(1.75, 0.75, 'Central tendency statistics')"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAxoAAAFfCAYAAADJSrqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAi4UlEQVR4nO3deXRW9Z3A4W8Ak0ASQGUX2a0KiKmIVq2giAsF6lZRaxWhtpyKCqPFVh13e6jLVCwqamdGrEp1RHGpg8iIyCh2hqqMihtqcBlRBCQQUBFy5w8PGWNAAv7wJeF5zsmR3Pfm3t99740nn9wleVmWZQEAAJBQg1wPAAAAqH+EBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQGwnZo0aVLk5eVVfRQWFka7du3iyCOPjD/+8Y+xcuXKLVrunDlz4rLLLovly5enHTAAdYrQANjOXXHFFXHnnXfGxIkT4+yzz46IiDFjxsRee+0VL7744mYvb86cOXH55ZcLDYDtXKNcDwCA3Bo4cGDsu+++VZ9fcMEFMXPmzBg8eHD8+Mc/jldffTUaN26cwxECUBc5owFADf3794+LL7443nnnnbjrrrsiIuLFF1+M008/Pbp06RKFhYXRpk2bGDFiRCxdurTq6y677LIYO3ZsRER07ty56rKshQsXRkTE7bffHv37949WrVpFQUFBdO/ePSZOnPidbx8AW58zGgBs0KmnnhoXXnhhPP744/GLX/wiZsyYEW+//XYMHz482rRpE/Pnz4/bbrst5s+fH3/7298iLy8vjjvuuHjjjTfiL3/5S1x//fXRokWLiIho2bJlRERMnDgxevToET/+8Y+jUaNG8cgjj8SZZ54ZlZWVMWrUqFxuLgCJ5WVZluV6EAB89yZNmhTDhw+PuXPnVrt06quaN28eXbp0ieeffz4+/fTTGpdQ3XPPPXHyySfH7Nmz4+CDD46IiOuuuy7Gjh0bZWVl0alTp2rzb2gZRx11VCxYsCDeeuutdBsHQM65dAqAjSouLq56+tRXA+Gzzz6LJUuWxA9+8IOIiHj++edrtbyvLqO8vDyWLFkS/fr1i7fffjvKy8sTjhyAXBMaAGxURUVFlJSURETEsmXLYvTo0dG6deto3LhxtGzZMjp37hwRUetIeOaZZ2LAgAFRVFQUzZs3j5YtW8aFF164WcsAoG5wjwYAG/T+++9HeXl5dOvWLSIihg4dGnPmzImxY8dGaWlpFBcXR2VlZRx11FFRWVm5yeW99dZbcdhhh8Uee+wRf/jDH2LXXXeN/Pz8+Pd///e4/vrra7UMAOoOoQHABt15550REXHkkUfGJ598Ek888URcfvnlcckll1TNs2DBghpfl5eXt8HlPfLII/H555/Hww8/HB06dKia/uSTTyYeOQDbApdOAVDDzJkz48orr4zOnTvHKaecEg0bNoyIiK8/P2T8+PE1vraoqCgiosYf7NvQMsrLy+P2229POHIAthXOaABs56ZNmxavvfZarF27Nj766KOYOXNmzJgxIzp27BgPP/xwFBYWRmFhYfTt2zeuueaa+OKLL2KXXXaJxx9/PMrKymosr3fv3hERcdFFF8VJJ50UO+ywQwwZMiSOOOKIyM/PjyFDhsTIkSOjoqIi/vSnP0WrVq1i0aJF3/VmA7CVCQ2A7dz6S6Hy8/Njp512ir322ivGjx8fw4cPr7oRPCJi8uTJcfbZZ8dNN90UWZbFEUccEdOmTYt27dpVW16fPn3iyiuvjFtuuSUee+yxqKysjLKysth9991jypQp8Y//+I/x61//Otq0aRO/+tWvomXLljFixIjvdJsB2Pr8HQ0AACA592gAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJNcr1ADZLlkWsXv3lv5s0icjLy+142Hpyta9zeYx9m3VXVkYsWfLlv1u0iGiwmb9D2B6/t7bHbc4l7zfAlquj/w+tW2c0Vq+OKC7+8mP9m039lKt9nctj7Nuse8mSiNatv/xYHxzf1brrqu1xm3PJ+w2w5ero/0PrVmgAAAB1gtAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByjWozU5ZlsXLlyq09lk1bter//71iRcS6dbkbC1tXrvZ1Lo+xb7Pur35/rlwZUVj43a27rtoetzmXvN8AW24b/X9oSUlJ5OXlbfT1vCzLsk0tZMWKFdGsWbOkAwMAAOqu8vLyaNq06UZfr1VobDNnNOLL6Nl1113jvffe+8YNgy3lGOO74Dhja3OMsbU5xtjUGY1aXTqVl5e3zR1ATZs23ebGRP3iGOO74Dhja3OMsbU5xtgYN4MDAADJCQ0AACC5OhcaBQUFcemll0ZBQUGuh0I95Rjju+A4Y2tzjLG1OcbYlFrdDA4AALA56twZDQAAYNsnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAKjjTj/99OjUqVOuh1ErdWms26qFCxdGXl5eTJo0KcnyZs2aFXl5eTFr1qwkywNYT2gA24233norRo4cGV26dInCwsJo2rRpHHTQQXHDDTfEp59+utXW+8EHH8Rll10W8+bN22rrqI2bb7452Q+nfLNXXnklLrvssli4cOEWL2Py5Mkxfvz4ZGOy/4HvWqNcDwDgu/Doo4/GCSecEAUFBXHaaadFz549Y82aNfH000/H2LFjY/78+XHbbbdtlXV/8MEHcfnll0enTp2itLR0q6yjNm6++eZo0aJFnH766Tkbw/bilVdeicsvvzwOOeSQLT6DM3ny5Hj55ZdjzJgx1aZ37NgxPv3009hhhx02a3kb2/99+/aNTz/9NPLz87donAAbIzSAeq+srCxOOumk6NixY8ycOTPatm1b9dqoUaPizTffjEcffTSHI6xu9erV0aRJk1wPg21UXl5eFBYWJltegwYNki4PYD2XTgH13jX
"text/plain": [
"<Figure size 1000x400 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Divide the figure in two subfigures\n",
"f, ax = plt.subplots(2, 1, figsize=(10,4), sharex=True) # The subfigures share the same x axis\n",
"\n",
"# In the first subfigure we show our observations\n",
"ax[0].set_yticks([]) # Do not show the ticks in the vertical axis\n",
"ax[0].set_ylabel(\"\") # Do not show any label in the vertical axis\n",
"ax[0].spines[['left', 'right', 'top']].set_visible(False) # Do not plot the left, right and top margins\n",
"ax[0].text(2.2, 0.75, \"Data\", fontsize=12)\n",
"\n",
"for i in range(len(Dat)):\n",
" ax[0].axvline(x=Dat[i], ymax=0.2, c='r') # Plot the observations as a rugplot\n",
"\n",
"# In the second subfigure we show the central tendency statistics\n",
"ax[1].spines[['left', 'right', 'top']].set_visible(False) # Do not plot the left, right and top margins\n",
"\n",
"for i in range(len(CentralTendency)):\n",
" ax[1].axvline(x=CentralTendency[i], ymax=0.2)\n",
"\n",
"ax[1].set_yticks([])\n",
"ax[1].set_ylabel(\"\")\n",
"ax[1].text(1.75, 0.75, \"Central tendency statistics\", fontsize=12)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# DispersionStatistics"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"S2Biased = np.var(Dat)\n",
"S2Unbiased = np.var(Dat, ddof=1)\n",
"\n",
"S1 = np.sqrt(S2Biased)\n",
"S2 = np.sqrt(S2Unbiased)\n",
"R = Dat.max() - Dat.min()\n",
"IQR = iqr(Dat)\n",
"MAD = median_abs_deviation(Dat)\n",
"AAD = abs(Dat-Dat.mean()).mean()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def an(n):\n",
" return np.sqrt((n-1)/2) * gamma((n-1)/2) / gamma(n/2)\n",
"\n",
"def c4(n):\n",
" return 1/an(n)\n",
"\n",
"# d2 for n in 2,...,25\n",
"d2 = {2:1.128, 3:1.693, 4:2.059, 5:2.326, 6:2.534, 7:2.704, 8:2.847, 9:2.970, 10:3.078, 11:3.173, 12:3.258, 13:3.336, 14:3.407, 15:3.472, 16:3.532, 17:3.588, 18:3.640, 19:3.689, 20:3.735, 21:3.778, 22:3.819, 23:3.858, 24:3.895, 25:3.931}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate unbiased estimators for the standard deviation"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"N = len(Dat)\n",
"\n",
"sigma_1 = S2 / c4(N)\n",
"sigma_2 = R / d2[N]\n",
"sigma_3 = IQR / (2 * norm.ppf(0.75))\n",
"sigma_4 = MAD / norm.ppf(0.75)\n",
"sigma_5 = AAD * np.sqrt(np.pi/2)\n",
"\n",
"sigma_list = [sigma_1, sigma_2, sigma_3, sigma_4, sigma_5]\n",
"labels = ['Deviation', 'Range', 'IQR', 'MAD', 'AAD']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we assume that the data follows a normal distribution, a prediction interval of prob. $(1-\\alpha)\\times100\\%$ would have limits corresponding with the quantiles of the normal distribution."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-0.05, 0.5)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAyEAAAH5CAYAAACWHDGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA850lEQVR4nO3deXyNZ/7/8ffJIZtEUmRBqag9ljaKalXVGlUjU2Mrra01tU2NtlQXSztUaYuO0poJulh/xtJ2aq+lQlGqqFgbqoTYt4gl5/794ZszTpMQJNcdyev5eJzHI7nOfV/35z7nFud97uu6b4dlWZYAAAAAwBAvuwsAAAAAkL8QQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQDYYurUqXI4HO6Hr6+vSpQooWbNmumjjz7SuXPnbqvftWvXaujQoTp9+nT2FgwAyDaEEACArd5++2198cUXmjhxovr27StJ6tevn6pVq6atW7fecn9r167VsGHDCCEAkIsVsLsAAED+1rx5cz300EPu3wcNGqTvvvtOTz31lP70pz8pPj5efn5+NlYIAMhunAkBAOQ6DRs21FtvvaUDBw7oyy+/lCRt3bpVXbp0UdmyZeXr66vw8HB169ZNJ06ccK83dOhQvfrqq5KkiIgI91Cv/fv3S5KmTJmihg0bKjQ0VD4+PqpSpYomTpxofP8AIL/jTAgAIFd69tln9frrr2vJkiV64YUXtHTpUv3666/q2rWrwsPD9csvv2jSpEn65Zdf9MMPP8jhcOjpp5/W7t27NWPGDI0ZM0bFihWTJIWEhEiSJk6cqMjISP3pT39SgQIF9PXXX6tXr15yuVzq3bu3nbsLAPmKw7Isy+4iAAD5z9SpU9W1a1dt3LjRYzjW9YKDg1W2bFlt3rxZFy9eTDcsa+bMmerQoYNWr16txx57TJL0/vvv69VXX1VCQoLKlCnjsXxGfURHR2vPnj3at29f9u0cAOCGGI4FAMi1AgIC3FfJuj48pKSk6Pjx43r44YclSZs3b85Sf9f3cebMGR0/flyPP/64fv31V505cyYbKwcA3AghBACQa50/f16BgYGSpJMnT+qll15SWFiY/Pz8FBISooiICEnKcoCIi4tT48aNVahQIQUHByskJESvv/76LfUBALhzzAkBAORKv//+u86cOaNy5cpJktq2bau1a9fq1Vdf1QMPPKCAgAC5XC5FR0fL5XLdtL99+/apUaNGqlSpkj788EOVKlVK3t7e+vbbbzVmzJgs9QEAyB6EEABArvTFF19Ikpo1a6ZTp05p+fLlGjZsmAYPHuxeZs+ePenWczgcGfb39ddf69KlS/rqq69UunRpd/uKFSuyuXIAwM0wHAsAkOt89913eueddxQREaGOHTvK6XRKkv54LZWxY8emW7dQoUKSlO5mhRn1cebMGU2ZMiUbKwcAZAVnQgAAtlq4cKF27typq1ev6ujRo/ruu++0dOlS3Xffffrqq6/k6+srX19f1a9fX6NGjdKVK1dUsmRJLVmyRAkJCen6q1mzpiTpjTfeUPv27VWwYEG1bNlSTZs2lbe3t1q2bKm//vWvOn/+vP71r38pNDRUiYmJpncbAPI1QggAwFZpw6u8vb1VpEgRVatWTWPHjlXXrl3dk9Ilafr06erbt68+/vhjWZalpk2bauHChSpRooRHf7Vq1dI777yjTz75RIsWLZLL5VJCQoIqVqyoOXPm6M0339Qrr7yi8PBw9ezZUyEhIerWrZvRfQaA/I77hAAAAAAwijkhAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCqgN0FII+yLCk5+drP/v6Sw2FvPUjP1HuU27fjcknHj1/7uVgxySsL383kteM7r+2PSbx2AOx0F/8N4kwIckZyshQQcO2R9o8DuYup9yi3b+f4cSks7NojLYzk1LZyq7y2Pybx2gGw0138N4gQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCqwJ12YFmWzp07lx21IC+5cOF/P589K6Wm2lcLMmbqPcrt27n+79e5c5Kvb85tK7fKa/tjEq8dADvl4r9BgYGBcjgcmT7vsCzLupMNnD17VkFBQXfSBQAAAIA85MyZMypcuHCmz99xCLn+TMjZs2dVqlQpHTx48IYbRf7A8YA0HAu4HscD0nAsIA3HQt5zszMhdzwcy+FwpDtYChcuzAEEN44HpOFYwPU4HpCGYwFpOBbyDyamAwAAADCKEAIAAADAqGwNIT4+PhoyZIh8fHyys1vcpTgekIZjAdfjeEAajgWk4VjIf+54YjoAAAAA3AqGYwEAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCqgN0FAAAAIH9ITU3VlStX7C4Dd6BgwYJyOp133A8hBAAAADnKsiwdOXJEp0+ftrsUZIPg4GCFh4fL4XDcdh+EEAAAAOSotAASGhoqf3//O/rwCvtYlqXk5GQlJSVJkooXL37bfRFCAAAAkGNSU1PdAaRo0aJ2l4M75OfnJ0lKSkpSaGjobQ/NYmI6AAAAckzaHBB/f3+bK0F2SXsv72R+DyEEAAAAOY4hWHlHdryXhBAAAAAARhFCAAAAAButXLlSDofjjq8ell39mEAIAQAAADLQpUsXORwOORwOFSxYUGFhYWrSpIkmT54sl8uVbdt55JFHlJiYqKCgoCyv06BBA/Xr1++O+7ELIQQAAADIRHR0tBITE7V//34tXLhQTzzxhF566SU99dRTunr1arZsw9vb+47vu5Gd/ZhACAGAfGro0KEe/1FNnTpVDodD+/fvz5b+s7u/u2XbpuWnfQXs4OPjo/DwcJUsWVJRUVF6/fXXtWDBAi1cuFBTp06VJJ0+fVrPP/+8QkJCVLhwYTVs2FA///yzJGn37t1yOBzauXO
"text/plain": [
"<Figure size 1000x600 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"f, ax = plt.subplots(2, 1, figsize=(10,6), sharex=True, height_ratios=[1, 3])\n",
"\n",
"for i in range(len(Dat)):\n",
" ax[0].axvline(x=Dat[i], ymax=0.2, c='r')\n",
"\n",
"ax[0].set_yticks([])\n",
"ax[0].set_ylabel(\"\")\n",
"ax[0].spines[['left', 'right', 'top']].set_visible(False)\n",
"ax[0].text(2, 0.75, \"Data\", fontsize=12)\n",
"\n",
"for i, sigma in enumerate(sigma_list):\n",
" ax[1].hlines(\n",
" 0.1*(4-i), \n",
" norm.ppf(ALPHA/2, Median, np.sqrt(1+np.pi/(2*N))*sigma), \n",
" norm.ppf(1-ALPHA/2, Median, np.sqrt(1+np.pi/(2*N))*sigma),\n",
" color=colors[i],\n",
" label=labels[i])\n",
" ax[1].scatter(Median, 0.1*(4-i))\n",
"\n",
"ax[1].legend()\n",
"ax[1].set_yticks([])\n",
"ax[1].set_ylabel(\"\")\n",
"ax[1].spines[['left', 'right', 'top']].set_visible(False)\n",
"ax[1].text(1.5, 0.45, r\"Prediction intervals, $\\alpha={}$\".format(ALPHA), fontsize=12)\n",
"ax[1].set_ylim(-0.05, 0.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Skewness"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Empirical central moments and k-statistics"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"m2 = S2Biased\n",
"m3 = ((Dat-Mean)**3).mean()\n",
"\n",
"k2 = S2Unbiased\n",
"k3 = N**2 / ((N-1)*(N-2))*m3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Coefficient of skewness based on central moments"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"g1_byhand = m3 / m2**(3/2)\n",
"g1 = skew(Dat)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(np.float64(1.472), np.float64(1.472))"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g1_byhand.round(4), g1.round(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Coefficient of skewness based on $k$-statistics"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"G1_byhand = k3 / k2**(3/2)\n",
"G1 = skew(Dat, bias=False)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(np.float64(1.7455), np.float64(1.7455))"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"G1_byhand.round(4), G1.round(4)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<Axes: ylabel='Density'>"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABW1klEQVR4nO3deXhTZdoG8DtJm3RP97SlhS7sa6EttexIZXMBUQecBeiM6Oc2OtVPxRlhHHUKigyfwoDLKOAG4zqOSxUrlcWytbKvLXShbdKNJm3aJm1yvj/SBCot0vVkuX/XlUs4OTl9ItDcPe/zvq9EEAQBRERERC5EKnYBRERERH2NAYiIiIhcDgMQERERuRwGICIiInI5DEBERETkchiAiIiIyOUwABEREZHLYQAiIiIil+MmdgH2yGw2o6ysDL6+vpBIJGKXQ0RERNdBEATU1dUhIiICUum17/EwALWjrKwMUVFRYpdBREREXVBSUoLIyMhrnsMA1A5fX18Alv+Bfn5+IldDRERE10On0yEqKsr2OX4tDEDtsA57+fn5MQARERE5mOtpX2ETNBEREbkcBiAiIiJyOQxARERE5HIYgIiIiMjlMAARERGRy2EAIiIiIpfDAEREREQuhwGIiIiIXA4DEBEREbkcBiAiIiJyOQxARERE5HIYgIiIiMjlMAARERGRy+Fu8ETUoTPqOuw+V4n9F2qgN7Sgf6AXBql8cee4SCi93MUuj4ioyxiAiOgqNXojXvjyFD7Ou9jm+I8F1QCAdd+dxbLJsbhncgy85Pw2QkSOh9+5iKiN3ecq8ccPfsKlhmZIJMCUQSFIiQtCiI8CxTUNyDyuxhlNHdbuOIsdJzV4a2kSQnwVYpdNRNQpEkEQBLGLsDc6nQ5KpRJarRZ+fn5il0PUZ344W4llWw/B2GLG0DBf/H3BKIzrH9DmHLNZwH+PluFv/z2Jar0RUYGe2JI2HrEhPiJVTURk0ZnPbzZBExEAYNcV4WfmcBU+f2jSVeEHAKRSCebF98PH909A/0AvlNQ04tdv7EdlnUGEqomIuoYBiIhQUtOAB9/Ps4Wf9b8eB7nbtb89RAd74+P7JyA2xBtqXRMefC8PzSZzH1VMRNQ9DEBELq7ZZMYft/2EuqYWjO3vf13hxyrEV4HXf5cIH4UbDhTW4IUvT/VytUREPYMBiMjFrd1xFj8V18LXww2vLBp73eHHamCoD9b+agwAYPOPhfjhbGVvlElE1KMYgIhc2PFSLTb9UAAAWH3HaEQFenXpOjNHhGHphGgAwIr/HEdTs6mnSiQi6hUMQEQuShAE/PXzExAE4NYxEZg7Krxb13ts5mCo/BQoqm7Ahp35PVQlEVHvYAAiclGfHynDoaJL8HSX4em5Q7t9PV8Pd/z11hEAgE0/FCC/or7b1yQi6i0MQEQuqMHYglVfnwYAPDAtDuFKzx657uyRYZg2JATNJgH/2HG2R65JRNQbGICIXNDWnCKUa5sQGeCJZVNie+y6EokET80ZCokE+PJYOU6W6Xrs2kREPYkBiMjFNBpNeHP3eQDAo6mD4eEu69HrDw3zw82t/UTrvuNdICKyT3YRgDZs2IDo6Gh4eHggOTkZBw4c6PDcTz75BImJifD394e3tzfi4+PxzjvvtDln6dKlkEgkbR6zZ8/u7bdB5BA+OFCMqnojIgM8MS8+ole+xqOpgyCVAN+e1ODYRW2vfA0iou4QPQBt374d6enpWLlyJfLy8jBmzBjMmjULFRUV7Z4fGBiIP//5z8jJycHRo0eRlpaGtLQ0fPPNN23Omz17NsrLy22PDz74oC/eDpFda2o24bVdlmnvD0wbCHdZ73wLGBjqi3nx/QAAr35/rle+BhFRd4gegNauXYtly5YhLS0Nw4cPx6ZNm+Dl5YW33nqr3fOnTZuG22+/HcOGDUNcXBweeeQRjB49Gnv27GlznkKhQFhYmO0REHD1nkZEruaj3IvQ6AwIV3rgjoR+vfq1Hpw+EACw45QGxdUNvfq1iIg6S9QAZDQakZubi9TUVNsxqVSK1NRU5OTk/OLrBUFAVlYWzpw5gylTprR5Ljs7G6GhoRgyZAjuv/9+VFdXd3gdg8EAnU7X5kHkbARBwNt7LwAAlk2OhcKtZ3t/fm5gqA+mDg6BIFhWiCYisieiBqCqqiqYTCaoVKo2x1UqFdRqdYev02q18PHxgVwux80334xXX30VN910k+352bNnY+vWrcjKysLq1avxww8/YM6cOTCZ2l+dNiMjA0ql0vaIiorqmTdIZEd+LKhGQaUe3nIZ7kqM7JOv+ftJMQCAfx8qQV1Tc598TSKi6+EmdgFd4evri8OHD6O+vh5ZWVlIT09HbGwspk2bBgBYtGiR7dxRo0Zh9OjRiIuLQ3Z2NmbMmHHV9ZYvX4709HTb73U6HUMQOZ2tOYUAgDsSIuHr4d4nX3PKoGDEhXijoFKPDw9dtAUiIiKxiXoHKDg4GDKZDBqNps1xjUaDsLCwDl8nlUoxcOBAxMfH47HHHsOdd96JjIyMDs+PjY1FcHAw8vPbX55foVDAz8+vzYPImZTWNmLHScu/s9/dMKDPvq5EIkHaREvo2fxjIcxmoc++NhHRtYgagORyORISEpCVlWU7ZjabkZWVhZSUlOu+jtlshsFg6PD5ixcvorq6GuHh3dvriMhRvbevCGYBmBAXhEEq3z792gvG9YOvwg3FNQ3Yd77jXjwior4k+iyw9PR0vPHGG9iyZQtOnTqF+++/H3q9HmlpaQCAxYsXY/ny5bbzMzIysGPHDpw/fx6nTp3Cyy+/jHfeeQe//e1vAQD19fX43//9X+zbtw+FhYXIysrCvHnzMHDgQMyaNUuU90gkpmaTGf8+dBEAsDil7+7+WHnJ3XBr63pD/z5U0udfn4ioPaL3AC1cuBCVlZVYsWIF1Go14uPjkZmZaWuMLi4uhlR6Oafp9Xo88MADuHjxIjw9PTF06FC8++67WLhwIQBAJpPh6NGj2LJlC2praxEREYGZM2fiueeeg0KhEOU9Eolp19lKVNUbEOwjx4xhql9+QS/4VWIU3t9fjK+Pq/FsYzOUnn3Tg0RE1BGJIAgclP8ZnU4HpVIJrVbLfiByePe/m4uvj6vxh0kxeOaW4aLUIAgCZq/bjTOaOjw3f2Sf9iERkevozOe36ENgRNR7LumNyDplWVX9zoS+mfreHolEYpt6/++DHAYjIvExABE5sf8eLYPRZMbwcD8MCxf3buaCcZFwl0lwrFSL02ouNkpE4mIAInJiH+damp/FvPtjFegtx/QhoQCAzw+XiVwNEbk6BiAiJ1VQWY8jF7Vwk0p6bdf3zrqttY7/Hi0D2w+JSEwMQERO6suj5QCAyYOCEeRjHzMgZwxVwUsuQ0lNIw6X1IpdDhG5MAYgIif1xVHLMNPNo+3j7g8AeMpluGm4ZSr+50c4DEZE4mEAInJCZzV1OKuph7tMYgsc9uLW1kD25dFymLg1BhGJhAGIyAl90Tr8NWVQiN0tOjhlsKWmijoD9l/g1hhEJA4GICInIwgCvmwd/rpljP3tfyd3k2L2CMtmx18dKxe5GiJyVQxARE7mtLoOBZV6yN2kSBVp64tfMnuUJQB9e0LDHeKJSBQMQERO5uvjagCW4S9fD/sa/rKaEBcEH4UbKuoMOHyxVuxyiMgFMQAROZkdJzUAgNkjw0SupGMKNxmmD7UsivjNCbXI1RCRK2IAInIiJTUNOFWug1QC3NgaMOzVrBGW4blvjqu5KCIR9TkGICInYr37kxQdiEBvucjVXNu0IaGQu0lRWN2As5p6scshIhfDAETkRKwByN7W/mmPj8INkwYGA+AwGBH1PQYgIidR22DEgcIaAMDM4fbb/3Ml6zCYNbgREfUVBiAiJ/H96QqYzAKGhvmif5CX2OVclxuHWgLQsVItKnRNIldDRK6EAYjISXx3ynGGv6xCfBUYE6kEAOw8UyFyNUTkShiAiJxAs8mM3WerAAAz7HTxw45Yp8N/f5oBiIj6DgMQkRPILbqEOkMLgrzlGN1PKXY5nTKjdRhsz7kqGFpMIldDRK6CAYjICViHj6Y
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sns.kdeplot(Dat)\n",
"\n",
"for i in range(len(Dat)):\n",
" plt.axvline(x=Dat[i], ymax=0.05, c='r') # Plot the observations as a rugplot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Deprecated measures of asymmetry"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"Sk1 = (Mean - Mode) / S2\n",
"Sk2 = (Mean - Median) / S2\n",
"S75 = (np.quantile(Dat, 0.75) + np.quantile(Dat, 0.25) - 2*np.median(Dat)) / IQR\n",
"S90 = (np.quantile(Dat, 0.9) + np.quantile(Dat, 0.1) - 2*np.median(Dat)) / (np.quantile(Dat, 0.9) - np.quantile(Dat, 0.1))"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(np.float64(0.2276),\n",
" np.float64(0.2276),\n",
" np.float64(0.2308),\n",
" np.float64(0.3895))"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Sk1.round(4), Sk2.round(4), S75.round(4), S90.round(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Kurtosis"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"m4 = ((Dat-Mean)**4).mean()\n",
"\n",
"k4 = N**2*((N+1)*m4 - 3*(N-1)*m2**2) / ((N-1)*(N-2)*(N-3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Coefficient of kurtosis based on central moments"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"g2_byhand = m4/m2**2\n",
"g2 = kurtosis(Dat, fisher=False)"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(np.float64(4.8753), np.float64(4.8753))"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g2_byhand.round(4), g2.round(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Coefficient of kurtosis based on $k$-statistics"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"G2_byhand = k4/k2**2 + 3\n",
"G2 = kurtosis(Dat, fisher=False, bias=False)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(np.float64(7.2795), np.float64(7.2795))"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"G2_byhand.round(4), G2.round(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Excess of kurtosis"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"g2_excess_byhand = g2_byhand - 3\n",
"g2_excess = kurtosis(Dat)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(np.float64(1.8753), np.float64(1.8753))"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g2_excess_byhand.round(4), g2_excess.round(4)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"G2_excess_byhand = G2_byhand - 3\n",
"G2_excess = kurtosis(Dat, bias=False)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(np.float64(4.2795), np.float64(4.2795))"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"G2_excess_byhand.round(4), G2_excess.round(4)"
]
}
],
"metadata": {
"kernelspec": {
2024-11-20 18:06:44 +07:00
"display_name": ".venv",
2024-10-16 10:21:36 +07:00
"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",
2024-11-19 23:03:48 +07:00
"version": "3.13.0"
2024-10-16 10:21:36 +07:00
}
},
"nbformat": 4,
"nbformat_minor": 2
}
