Lab13-Exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimating the parameters of a distribution\n",
    "\n",
    "## Exercise 1\n",
    "Linked Bernoulli distributions\n",
    "\n",
    "Consider two coins. The probability of getting \"heads\" is $p$ for the first coin and $2p$ for the second. We then toss the first coin 5 times and the second 10 times, as simulated in the code below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X: [0 1 0 0 0]\n",
      "Y: [1 1 1 1 1 1 1 0 1 1]\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import bernoulli\n",
    "p = 0.3\n",
    "X = bernoulli.rvs(p, size=5, random_state=1)\n",
    "Y = bernoulli.rvs(2*p, size=10, random_state=2)\n",
    "print('X:', X) # 1 = heads; 0 = tails\n",
    "print('Y:', Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Plot the log-likelihood of the data as a function of $\\hat{p}$. Note that the data consists of both `X` and `Y`, so the likelihood function becomes $L(\\hat{p} | X,Y)$\n",
    "1. Experimentally determine the MLE estimation for $p$ corresponding to the observations in `X` and `Y`.\n",
    "1. Analytically determine the MLE estimation for $p$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "Poisson distribution\n",
    "\n",
    "A call centre keeps track of the number of phone calls received every day. In order to accurately plan the resources, the number of calls for the last 100 days are modelled as a random variable $X$ following a Poisson distribution of an unknown parameter $\\lambda$.\n",
    "\n",
    "We will simulate this for $\\lambda = 4$ using the `poisson.rvs` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 2, 3, 4, 4, 6, 1, 5, 3, 5])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import poisson\n",
    "lambda_ = 4\n",
    "X = poisson.rvs(lambda_, size=100, random_state=1)\n",
    "X[:10] # Calls received in the first 10 days"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Plot the histogram of the data.\n",
    "1. Plot the log-likelihood of the data as a function of $\\hat{\\lambda}$.\n",
    "1. Experimentally determine the MLE estimation for $\\lambda$ corresponding to the observations in `X`.\n",
    "1. Analytically determine the MLE estimation for $\\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "Uniform distribution\n",
    "\n",
    "Consider a hashing function that returns a number in the interval $[-w, w]$ for any file. Any value in that interval is equally likely to appear so the hash values are following a uniform distribution $U(-w,w)$. We can simulate the hashes for 100 files using the `uniform.rvs` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.65955991,  4.40648987, -9.9977125 ])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import uniform\n",
    "w = 10\n",
    "X = uniform.rvs(-w, 2*w, size=100, random_state=1)\n",
    "X[:3] # The first 3 hashes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Plot the histogram of the data.\n",
    "1. Experimentally determine the MLE estimation for $w$ given the observations in `X`.\n",
    "1. Analytically determine the MLE estimation for $w$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4\n",
    "Exponential distribution\n",
    "\n",
    "Seismologists are tracking the time interval between consecutive major earthquakes. They noticed that it follows an exponential distribution $Exp(\\lambda)$.\n",
    "\n",
    "To simulate the observed intervals between 100 earthquakes that occur on average once per year, we can use the `expon.rvs` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5.39605837e-01, 1.27412525e+00, 1.14381359e-04])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import expon \n",
    "lambda_ = 1 # Once per year, on average\n",
    "X = expon.rvs(scale=1/lambda_, size=100, random_state=1)\n",
    "X[:3] # The first 3 intervals "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Plot the histogram of the data.\n",
    "1. Experimentally determine the MLE estimation for $\\lambda$ corresponding to the observations in `X`.\n",
    "1. Analytically determine the MLE estimation for $\\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5\n",
    "Variance of a Gaussian distribution\n",
    "\n",
    "Consider $X$ a random variable representing size of pollen grains following a normal (Gaussian) distribution with known mean 0 and a variance $\\sigma^2$, formally written as $X \\sim N(0, \\sigma^2)$. Load the observations for this variable as a `numpy` array, by running the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_openml\n",
    "pollen = fetch_openml('pollen', version=1, as_frame=False, parser='auto')\n",
    "X = pollen.data[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Plot the histogram of the data.\n",
    "1. Experimentally find the value of $\\hat{\\sigma}^2_\\text{MLE}$ by testing candidates in the interval $[1, 10]$. Note that since the dataset is quite large, calculating the likelihood of the data can quickly result in an underflow on most systems. Try using the log-likelihood instead.\n",
    "1. Analytically find the estimator $\\hat{\\sigma}^2_\\text{MLE}$ and apply the resulting formula to the dataset."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}