{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normal_cdf(x, loc=0.0, scale=1.0):\n",
    "    \"\"\"\n",
    "    Compute the CDF of the normal distribution using only NumPy\n",
    "    (no SciPy, no math.erf).\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : float or np.ndarray\n",
    "        Points at which to evaluate the CDF.\n",
    "    loc : float\n",
    "        Mean of the normal distribution.\n",
    "    scale : float\n",
    "        Standard deviation of the normal distribution.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    np.ndarray\n",
    "        The CDF values.\n",
    "    \"\"\"\n",
    "    # Standardize\n",
    "    z = (x - loc) / scale\n",
    "    t = 1.0 / (1.0 + 0.2316419 * np.abs(z))\n",
    "    \n",
    "    # Coefficients for the approximation\n",
    "    a1, a2, a3, a4, a5 = 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429\n",
    "    poly = (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t\n",
    "    pdf = np.exp(-0.5 * z**2) / np.sqrt(2 * np.pi)\n",
    "    cdf = 1.0 - pdf * poly\n",
    "\n",
    "    # Use symmetry for negative z\n",
    "    cdf = np.where(z >= 0, cdf, 1.0 - cdf)\n",
    "    return cdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task 1. Probing the Central Limit Theorem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting with a Gaussian random distribution (`np.random.normal`), with `mean=0` and `std=1`, find out by brute force (i.e. use a Monte Carlo Method, i.e. generate 1000000 random numbers and test it) what the probability is that a number is between -1 and 1 (This should be approximately 68%)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's generate 1000000 random numbers from a normal distribution\n",
    "# with mean 0 and standard deviation 1\n",
    "sample = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the probability of having a number between -1 and 1\n",
    "we can write a boolean expression that checks if the number is\n",
    "between -1 and 1. This can be done using numpy's `np.logical_and`\n",
    "function. This function takes two boolean arrays and returns\n",
    "a new boolean array where the i-th element is True if both\n",
    "the i-th elements of the input arrays are True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = ... # Create a boolean mask where the condition is met\n",
    "\n",
    "# The probability of having a number between -1 and 1 is the\n",
    "# number of elements in the mask that are True divided by the\n",
    "# total number of elements in the sample.\n",
    "p = ...\n",
    "print(\"Probability of having a number between -1 and 1: \", p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create a low number N (start with 4) of such random samples. Repeat this `M=1000000` times (i.e. make a 2D array, let's call that x2DArray). Compute the mean (let's call that `xBarArray=np.mean(x2DArray,axis=0)`). Draw a histogram of the distribution of these means. Find the mean of these means (`np.mean(xBarArray)`) and the std. of these means. How does that compare to the expectation from the Central Limit Theorem? What is the probability P of one such mean (one such value in  xBarArray) to be within -1 to +1?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's create a 2D array of random normal numbers\n",
    "# with 4 rows and 1000000 columns\n",
    "\n",
    "N = 4\n",
    "M = 1000000\n",
    "x2DArray = ...\n",
    "\n",
    "# To compute the mean along the rows we can use the np.mean:\n",
    "xBarArray = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can plot the histogram of the mean values\n",
    "... # Plot histogram of xBarArray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean and std of the mean values\n",
    "mean_of_means = ...\n",
    "std_of_means = ...\n",
    "\n",
    "print(\"Computed mean of the means: \", mean_of_means)\n",
    "print(\"Computed std of the means: \", std_of_means)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the central limit theorem, the mean of the means\n",
    "should be equal to the mean of the original distribution and\n",
    "the standard deviation of the means should be equal to the\n",
    "standard deviation of the original distribution divided by\n",
    "the square root of the number of samples in the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Mean of the original distribution: \", ...)\n",
    "print(\"Std of the original distribution: \", ...)\n",
    "print(\"Std of the means according to the CLT: \", ...)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Probability of the means being between -1 and 1:\n",
    "mask = ...\n",
    "p = ...\n",
    "print(\"Probability of having a mean between -1 and 1 for N = 4: \", p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now find that probability P as a function of N (from N=2 to N=100) and make a plot of P(N). Again compare this to what you expect from the Central Limit Theorem. Optional: plot a line on top with the expected analytic result. \n",
    "\n",
    "**Hint**: Analytical probability can be computed as follows:\n",
    "\n",
    "$$P(-1 < X < 1) = P(X < 1) - P(X < -1) = F(1) - F(-1)$$\n",
    "\n",
    "where $F$ is the cumulative distribution function of the normal distribution.\n",
    "According to the central limit theorem, the standard deviation in this case\n",
    "is $1/\\sqrt{N}$. $F$ can be computed using `normal_cdf` function with `loc=0` and `scale=1/np.sqrt(N)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now, let's compute the probability of having a mean between -1 and 1\n",
    "# as a function of the number of samples N. For each N from 2 to 100\n",
    "# we will (1) create the 2D array of random normal numbers, (2) compute\n",
    "# the mean along the rows, (3) compute the probability of having a mean\n",
    "# between -1 and 1.\n",
    "\n",
    "# Note: This may take about a minute to compute.\n",
    "\n",
    "N_range = ... # N from 2 to 100 with step size 5\n",
    "p_values = []\n",
    "p_values_analytical = []\n",
    "for N in N_range:\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's plot the probability of having a mean between -1 and 1\n",
    "# as a function of N.\n",
    "\n",
    "plt.plot(...)\n",
    "plt.plot(...)\n",
    "plt.xlabel(\"N (number of samples in the mean)\")\n",
    "plt.ylabel(\"Probability of having a mean between -1 and 1\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compute the std of each set of N numbers (i.e. take `xStdArray = np.std(x2DArray,axis=0)`), and find the mean of this std (`np.mean(xStdArray)` ). How does this compare to what you expect? Plot the result as a function of N (from N=2 to N=100). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's repeat the same trick with the standard deviation\n",
    "# and analyze how the standard deviation behaves as a function\n",
    "# of the number of samples N.\n",
    "N = 4\n",
    "M = 1000000\n",
    "x2DArray = ...\n",
    "\n",
    "# To compute the mean along the rows we can use the np.mean:\n",
    "xStdArray = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we can plot the histogram of the std values\n",
    "plt.hist(...)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean and std of the std values\n",
    "print(\"Computed mean of the std: \", ...)\n",
    "print(\"Computed std of the std: \", ...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can study how the mean of the standard deviation\n",
    "# behaves as a function of the number of samples N.\n",
    "\n",
    "N_range = ...\n",
    "mean_values_of_std = []\n",
    "for N in N_range:\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's plot the probability of having a mean between -1 and 1\n",
    "# as a function of N.\n",
    "\n",
    "plt.plot(...)\n",
    "plt.xlabel(\"N (number of samples in the mean)\")\n",
    "plt.ylabel(\"Mean of the standard deviation\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task 2. Who's taller?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a certain populations, the height of men is described by a normal distribution of average 180cm and std 10cm, whereas the height of women is on average 170cm with a std of 10cm. Find the probability of one randomly chosen man to be taller than one randomly chosen woman. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "men_heights = ...\n",
    "women_heights = ...\n",
    "\n",
    "# To calculate the probability that a randomly selected man\n",
    "# is taller than a randomly selected women, we need to calculate\n",
    "# the probability that the difference between the heights is positive.\n",
    "\n",
    "difference = ...\n",
    "p = ...\n",
    "\n",
    "print(\"Probability of a man being taller than woman: \", p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now find the probability P that the mean height of N=3 randomly chosen men is larger than the mean of 3 randomly chosen women."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's calculate the probability of 3 randomly selected\n",
    "# men being taller on average than 3 randomly selected women.\n",
    "\n",
    "N = 3\n",
    "M = 100000\n",
    "men_heights = ...\n",
    "women_heights = ...\n",
    "\n",
    "mean_three_men_height = ...\n",
    "mean_three_women_height = ...\n",
    "\n",
    "difference = ...\n",
    "p = ...\n",
    "\n",
    "print(\"Probability of 3 radomly selected men being taller\")\n",
    "print(\"on average than 3 randomly selected women: \", p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find how this probability depends on N, plot P(N) (optionally: as a log plot of 1-P). How large does N need to be so that you find with 99.7% probability, that the mean height of men is bigger than the mean height of women?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_range = ...\n",
    "p_range = []\n",
    "for N in N_range:\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's plot the probability of having a mean between -1 and 1\n",
    "# as a function of N.\n",
    "\n",
    "plt.plot(...)\n",
    "plt.xlabel(\"N (number of samples in the mean)\")\n",
    "plt.ylabel(\"P of N men being taller N women on average\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finding the number of samples needed to have a probability\n",
    "# of 99.7% that the mean of the men's heights is greater than\n",
    "# the mean of women's heights.\n",
    "\n",
    "p_target = 0.997\n",
    "p_range = np.arange(p_range) # Convert to numpy array\n",
    "index = ...\n",
    "N_target = N_range[index]\n",
    "\n",
    "print(\"Number of samples needed to have a probability of 99.7%: \", N_target)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}