Question 2- Logistic Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homewok 2: Logistic Regression\n",
    "\n",
    "Welcome to the course **AI and Deep learning**!\n",
    "\n",
    "In practice, categorical data are more commonly analyzed, including but not limited to fraud detection, face recognition in Xiamen University and object  detection. As a fundation for categorical data analysis, logistic regression is of great importance in both statistics and deep learning. In the second homework, we would use simulation again to review logistic regression in more details. Computation graph is implicitly used to obtain the derivatives and update the estimated parameters. Hope you enjoy the second homework!  \n",
    "\n",
    "**Learning Goal**: In this homework, we are going to conduct a simple simulation study based on a logistic regression. After this homework, you will know:\n",
    "    * How to estimate the model parameters by gradient-based methods. \n",
    "    * How Newton-Raphson method works. \n",
    "    * How (batch) gradient descent method works. \n",
    "    * How important the learning rate is for a (batch) gradient descent method.\n",
    "    * How to make inference for the parameters in a logistic regression model.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of content\n",
    "* [1 - Packages](#1)\n",
    "* [2 - Generate a training dataset](#2)\n",
    "* [3 - Parameter estimation](#3) \n",
    "    * [3.1 - Newton-Raphson method](#3.1)\n",
    "    * [3.2 - Gradient descent method](#3.2)\n",
    "* [4 - Statistical inference](#4)\n",
    "    * [4.1 - Bias, variance and MSE](#4.1)\n",
    "    * [4.2 - Coverage rate](#4.2)\n",
    "* [5 - Closer check for the two estimation methods](#5)\n",
    "* [6 - Additional homework](#6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='1'></a>\n",
    "## 1- Packages\n",
    "\n",
    "In order to finish a task, we need commands from certain **Python** packages. Again, one of the commonly used package is **numpy**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt # for plots\n",
    "import pandas as pd # for printing tables\n",
    "import scipy.stats # for quantile calculation when we consider coverage rates\n",
    "import time # for timing the code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='2'></a>\n",
    "## 2 - Generate a training dataset\n",
    "\n",
    "First, we generate a training dataset from a pre-specified logistic regression model. **In order to guarantee that our simulation results are reproducible, we need to control the random seed.** That is, after controlling the seed, others can generate the **SAME** random variables as we did, so our simulation results can be reproduced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following logistic regression model \n",
    "$$\n",
    "y^{(i)}\\sim\\mbox{Bernoulli}\\{\\pi(x^{(i)})\\},\\\\\n",
    "\\mathrm{logit}\\{\\pi(x^{(i)})\\} = b_0 +w_{00}x^{(1i)}+w_{01}x^{(2i)},\n",
    "$$\n",
    "where $\\mbox{Bernoulli}(p)$ is a Bernoulli distribution with success probability $p$, $x^{(i)} = (x^{(1i)},x^{(2i)})^T$, $b_0=-0.5$, $w_{00}=0.1$, $w_{01}=-0.1$, $x^{(ki)}\\sim N(2,2^2)$ $k=1,2$.\n",
    "Let us write a function to generate a training dataset of size $n$ with a random number $rn$. \n",
    "\n",
    "First, we need a sigmoid function $\\sigma(x) = 1/\\{1+\\exp(-x)\\}$, and complete the following function (**DO NOT CHANGE THE EXISTING PARTS**). The following command may be useful, and check the help document for details:\n",
    "\n",
    "* `np.exp`: Exponential function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE EXISTING CODE\n",
    "# Provide your code between ``YOUR CODE BEGINS HERE (approximately 1 lines)'' and ``YOUR CODE ENDS''.\n",
    "def sigmoid(x):\n",
    "    # x: input\n",
    "    \n",
    "    ### YOUR CODE BEGINS HERE (approximately 1 lines)\n",
    "    sig = \n",
    "    ### YOUR CODE ENDS    \n",
    "\n",
    "    return sig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, compete the following function. (**DO NOT CHANGE THE EXISTING PARTS**) The following command may be useful, and check the help document for details:\n",
    "\n",
    "* `np.random.binomial`: generate normally distributed random variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE EXISTING CODE\n",
    "# Provide your code between ``YOUR CODE BEGINS HERE'' and ``YOUR CODE ENDS''.\n",
    "def train_data_generation(n, rn):\n",
    "    # n: sample size\n",
    "    # rn: random seed\n",
    "    \n",
    "    # Step 1. Set random seed to be rn\n",
    "    # Step 2. Generate x of size nX2 from a normal distribution with mean 2 and standard deviation 2\n",
    "    # Step 3. Calculate z by z=b_0 + Xw_0\n",
    "    # Step 4. Calculate a by a = \\sigma(z)\n",
    "    # Step 5. Generate 2d y of size nX1\n",
    "    \n",
    "    ### YOUR CODE BEGINS HERE (approximately 5--6 lines)\n",
    "\n",
    "    x = \n",
    "    z = \n",
    "    a = \n",
    "    y =\n",
    "    ### YOUR CODE ENDS\n",
    "    \n",
    "    \n",
    "    return x, y\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize your data, you may would like to run the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "x, y = train_data_generation(1000, 100)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "scatter = ax.scatter(x[:,0], x[:,1],  c=y[:,0])\n",
    "legend1 = ax.legend(*scatter.legend_elements(),\n",
    "                    loc=\"lower right\", title=\"Classes\")\n",
    "ax.set_title('Simulated training dataset')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='3'></a>\n",
    "## 3 - Parameter estimation\n",
    "\n",
    "The estimated parameters $\\hat{b}$ and $\\hat{w}$ are obtained by minimizing the following cost function,\n",
    "$$\n",
    "J(b,w) = -n^{-1}\\sum_{i=1}^n\\{y^{(i)}\\log a^{(i)} + (1-y^{(i)})\\log(1-a^{(i)})\\}, \\tag{1}\n",
    "$$\n",
    "where $a^{(i)} = \\sigma(z^{(i)})$ and $z^{(i)} = (x^{(i)})^Tw+b$.\n",
    "\n",
    "Different from linear regression, we do not have a 'normal equation' for a logistic regression model. To find parameters minimizing the cost function (1), we have two different gradient-based methods, including a Newton-Raphson method and a (batch) gradient descent method. We would consider both of them in this homework, and  both methods are iterative. The main difference is that a Hessian matrix is involved for the Newton-Raphson method, but only gradients are required for the gradient descent method (with the Hessian matrix replaced by a **learning rate** $\\alpha$). \n",
    "\n",
    "\n",
    "\n",
    "For both methods, denote $b^{(0)}$ and $w^{(0)}$ as the initials when estimating both parameters $b_0$ and $w_0=(w_{00},w_{01})^T$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='3.1'></a>\n",
    "### 3.1 - Newton-Raphson method\n",
    "\n",
    "Newton-Raphson method estimates model parameters based on the following steps. \n",
    "\n",
    "- Step 1. Initialize $(b^{\\{0\\}}, w^{\\{0\\}})$.\n",
    "- Step 2. Given $(b^{\\{t\\}}, w^{\\{t\\}})$ for $t=0,1,\\ldots$, obtain \n",
    "$$\\nabla J(b^{\\{t\\}},w^{\\{t\\}}) = n^{-1}\\sum_{i=1}^n(a^{(i)}-y^{(i)})\\tilde{x}^{(i)}, \\quad H(b^{\\{t\\}},w^{\\{t\\}})=\\frac{\\partial^2 J}{\\partial (b,w)^T\\partial (b,w)}(b^{\\{t\\}},w^{\\{t\\}})=n^{-1}\\sum_{i=1}^na^{(i)}(1-a^{(i)})\\tilde{x}^{(i)}(\\tilde{x}^{(i)})^T. $$\n",
    "Notice that we do not use a superscript '$\\{t\\}$' for $a^{(i)}$'s. \n",
    "- Step 3. Update the model parameters by \n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "b^{\\{t+1\\}}\\\\\n",
    "w^{\\{t+1\\}}\n",
    "\\end{pmatrix} = \\begin{pmatrix}\n",
    "b^{\\{t\\}}\\\\\n",
    "w^{\\{t\\}}\n",
    "\\end{pmatrix} - \\left\\{H(b^{\\{t\\}},w^{\\{t\\}})\\right\\}^{-1}\\nabla J(b^{\\{t\\}},w^{\\{t\\}})\n",
    "$$\n",
    "- Step 4. Go back to Step 2 until convergence.\n",
    "\n",
    "When we code the Newton-Raphson method, we usually augment the design matrix by a column of 1's as what we did for linear regression models, so that we treat the two parameters by a single one. That is, denote $X=(\\tilde{x}^{(1)},\\ldots,\\tilde{x}^{(1)})^T$, where $\\tilde{x}^{(i)} = (1, (x^{(i)})^T)^T$ for $i=1,\\ldots,n$. Then, we treat  $\\tilde{w}=(b,w^T)^T$ as the single model parameter vector. Please notice that we use $\\tilde{X}$ to denote the augmented design matrix in the slides. \n",
    "\n",
    "\n",
    "**Note:**\n",
    "* In Step 2, it is easy to show the following results: ( _check by yourself_ )\n",
    "$$\n",
    "\\nabla J(\\tilde{w}) = \\frac{\\partial J}{\\partial \\tilde{w}} = n^{-1}\\sum_{i=1}^n(a^{(i)} - y^{(i)})\\tilde{x}^{(i)},\\\\\n",
    "H(\\tilde{w}) = \\frac{\\partial^2 J(\\tilde{w})}{\\partial \\tilde{w}^T\\partial \\tilde{w}} = n^{-1}\\sum_{i=1}^na^{(i)} (1- a^{(i)})\\tilde{x}^{(i)}(\\tilde{x}^{(i)})^T,\\\\\n",
    "$$\n",
    "* Summation undermines the computation efficiency, so we consider vectorization. It can be shown that the above results are equivalent to ( _check by yourself_ )\n",
    "$$\n",
    "\\nabla J(\\tilde{w}) =n^{-1}X^T(A-Y), \\quad H(\\tilde{w}) = n^{-1}X^TWX,\n",
    "$$\n",
    "where $Y=(y^{(1)},\\ldots,y^{(n)})^T$, $A=(a^{(1)},\\ldots,a^{(n)})^T$ and $W = \\mbox{diag}((a^{(1)}(1-a^{(1)}),\\ldots, a^{(n)}(1-a^{(n)})))$. \n",
    "\n",
    "<font color='red'>Notice that multiplying a diagonal matrix is still computationally inefficient. Do you know how to further improve the computation efficiency?</font>\n",
    "* In Step 4, there are many different metrics to check the convergence. For example, if the Euclidean norm between $(b^{(t)}, w^{(t)})$ and $(b^{(t+1)}, w^{(t+1)})$ is less than a tolerance, say $10^{-4}$.   \n",
    "* <font color='red'>The maximum iteration is 1\\,000</font>, and you can modify this number by yourself. However, it is proved that Newton-Raphson method converges very fast. If you are interested, the Newton-Raphson method achieves quadratic convergence; check relevant references for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we code up the Newton-Raphson method mentioned above. <font color='red'>Your output should be a $3\\times1$ 2d array.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE EXISTING CODE\n",
    "def NR_est(x,y, tol = 10**(-8)):\n",
    "    # x: the generated feature vector of size nX2\n",
    "    # y: the generated feature vector of length n\n",
    "    # tol: tolerance for convergence criterion\n",
    "    \n",
    "    # Step 1. Obtain the sample size\n",
    "    # Step 2. Augment the feature vector with the intercept term, check the relationship between \\tilde{x} and x\n",
    "    # Step 3. Initialize the parameter for \\tilde{w}. Please notice that the size of this parameter is 3X1 2d array\n",
    "    # Step 4. Initialize a counter for iterations\n",
    "    # Step 5. Main loop. If number of iterations is less than 1000 and the norm between consecutive estimators is larger than the \n",
    "    #         tolerance, then do the following:\n",
    "    #         Step 5.1. Obtain z using the current parameter\n",
    "    #         Step 5.2. Obtain a using the current parameter\n",
    "    #         Step 5.3. Obtain \\nabla J\n",
    "    #         Step 5.4. Obtain the Hessian matrix H\n",
    "    #         Step 5.5. Update the parameter by the Newton-Raphson method\n",
    "    #         Step 5.6. add 1 to the counter\n",
    "    \n",
    "    ### YOUR CODE BEGINS HERE (approximately 13 lines)\n",
    "    n = \n",
    "    aug_x = \n",
    "\n",
    "    par = \n",
    "    \n",
    "    iter_index = \n",
    "    while :\n",
    "        z = \n",
    "        A = \n",
    "        \n",
    "        nabla_J = \n",
    "        H_J = \n",
    "        \n",
    "        par -= \n",
    "        iter_index += \n",
    "    ### YOUR CODE ENDS\n",
    "    \n",
    "    return par\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "x, y = train_data_generation(50000, 1010)\n",
    "print('Your estimators by the Newton-Raphson method are:')\n",
    "print(NR_est(x,y).flatten())\n",
    "print('The true values are\\n (-0.5, 0.1, -0.1)')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your Newton-Raphson method is correct, your estimator should be very close to the true value. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='3.2'></a>\n",
    "### 3.2 - Gradient descent method\n",
    "\n",
    "In this section, we consider the gradient descent method, which is mainly based on the following steps. \n",
    "\n",
    "- Step 1. Initialize $(b^{\\{0\\}}, w^{\\{0\\}})$.\n",
    "- Step 2. Given $(b^{\\{t\\}}, w^{\\{t\\}})$ for $t=0,1,\\ldots$, obtain \n",
    "$$\\nabla J(b^{\\{t\\}},w^{\\{t\\}}) = n^{-1}\\sum_{i=1}^n(a^{(i)}-y^{(i)})\\tilde{x}^{(i)}$$\n",
    "- Step 3. Update the model parameter by \n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "b^{\\{t+2\\}}\\\\\n",
    "w^{\\{t+1\\}}\n",
    "\\end{pmatrix} = \\begin{pmatrix}\n",
    "b^{\\{t\\}}\\\\\n",
    "w^{\\{t\\}}\n",
    "\\end{pmatrix} -\\alpha\\nabla J(b^{\\{t\\}},w^{\\{t\\}}),\n",
    "$$\n",
    "where $\\alpha$ is a pre-specified learning rate.\n",
    "- Step 4. Go back to Step 2 until convergence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note:\n",
    "* Again, we should use vectorization to calculate $\\nabla J(b^{\\{t\\}},w^{\\{t\\}})$ as before.\n",
    "* The difference between the Newton-Raphson method and the gradient descent method is that a simple learning rate is used for the latter, but a complex Hessian matrix is needed for the former.\n",
    "* The maximum iteration is 5\\,000, which is larger than the one for the Newton-Raphson method above since the gradient descent method converges much slower."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we code up the Newton-Raphson method mentioned above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE EXISTING CODE\n",
    "def BGD_est(x,y, alpha=10**(-2), tol = 10**(-8)):\n",
    "    # x: the generated feature vector of length n\n",
    "    # y: the generated feature vector of length n\n",
    "    # alpha: learning rate\n",
    "    # tol: tolerance for convergence criterion\n",
    "    \n",
    "    # Step 1. Obtain the sample size\n",
    "    # Step 2. Augment the feature vector with the intercept term, check the relationship between \\tilde{x} and x\n",
    "    # Step 3. Initialize the parameter for \\tilde{w}. Please notice that the size of this parameter is 3X1 2d array\n",
    "    # Step 4. Initialize a counter for iterations\n",
    "    # Step 5. Main loop. If number of iterations is less than 5000 and the norm between consecutive estimators is larger than the \n",
    "    #         tolerance, then do the following:\n",
    "    #         Step 5.1. Obtain z using the current parameter\n",
    "    #         Step 5.2. Obtain a using the current parameter\n",
    "    #         Step 5.3. Obtain \\nabla J\n",
    "    #         Step 5.4. Update the parameter by the gradient descent method\n",
    "    #         Step 5.5. add 1 to the counter\n",
    "    \n",
    "    ### YOUR CODE BEGINS HERE (approximately 10 lines)\n",
    "    n = \n",
    "    aug_x = \n",
    "\n",
    "    par = \n",
    "    \n",
    "    iter_index = \n",
    "    while :\n",
    "        z = \n",
    "        A = \n",
    "        \n",
    "        nabla_J = \n",
    "        \n",
    "        par -= \n",
    "        iter_index += \n",
    "    ### YOUR CODE ENDS\n",
    "    \n",
    "    return par"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "x, y = train_data_generation(50000, 1010)\n",
    "print('Your estimators by the gradient descent method are:')\n",
    "print(BGD_est(x,y).flatten())\n",
    "print('The true values are\\n (-0.5, 0.1, -0.1)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>Different from the Newton-Raphson method, we do not expect that estimators are VERY close to the truth, especially when the learning rate and the maximum number of iterations are both small, since performance of the gradient descent method mainly depends on the learning rate. Different learning rates lead to different estimation result. Try some different values by yourself. </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='4'></a>\n",
    "## 4 - Statistical inference\n",
    "Usually, we conduct 1\\,000 Monte Carlo simulations. That is, we generate the training datasets of size $n$ independently 1\\,000 times. Then, for each generated training dataset, we estimate the model parameters $(b_0,w_0)$. Then, we pool them together to check their bias, variance and MSE. If we have more than one estimation methods, we can compare when in terms of bias, variance and MSE. In this homework, we have two different methods to estimate the model parameters, and we would like to compare the statistics in this section and computation efficiency in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='4.1'></a>\n",
    "### 4.1 - Bias, variance and MSE\n",
    "We first consider the case with $n=1\\,000$, and record all the estimated parameters.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "n = 1000\n",
    "par_NR_result = np.zeros((1000,3))\n",
    "par_GD_result = np.zeros((1000,3))\n",
    "\n",
    "for i in range(1000):\n",
    "    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i\n",
    "    par_NR_result[i,:] = NR_est(x,y).flatten()\n",
    "    par_GD_result[i,:] = BGD_est(x,y).flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, check the bias of the two methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "true_par = np.array([-0.5, 0.1, -0.1])\n",
    "bias_NR = np.mean(par_NR_result,axis = 0) - true_par \n",
    "bias_GD = np.mean(par_GD_result,axis = 0) - true_par\n",
    "\n",
    "bias_all = np.concatenate((bias_NR, bias_GD)).reshape((2,3))\n",
    "df=pd.DataFrame(bias_all, index=['Newton-Raphson', 'Gradient descent'],columns=['b','w0','w1'])\n",
    "df.style.set_caption('Bias of the two methods')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, check the variance of the two methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "var_NR = np.var(par_NR_result,axis = 0)\n",
    "var_GD = np.var(par_GD_result,axis = 0)\n",
    "\n",
    "var_all = np.concatenate((var_NR, var_GD)).reshape((2,3))\n",
    "df=pd.DataFrame(var_all, index=['Newton-Raphson', 'Gradient descent'],columns=['b','w0','w1'])\n",
    "df.style.set_caption('Variance of the two methods')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, check the MSE of the two methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "MSE_NR = bias_NR**2 + var_NR\n",
    "MSE_GD = bias_GD**2 + var_GD\n",
    "\n",
    "MSE_all = np.concatenate((MSE_NR, MSE_GD)).reshape((2,3))\n",
    "df=pd.DataFrame(MSE_all, index=['Newton-Raphson', 'Gradient descent'],columns=['b','w0','w1'])\n",
    "df.style.set_caption('MSE of the two methods')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your code is correct, we can observe that the bias of the gradient method is slightly larger than the Newton-Raphson method, and their variance is kind of comparable. Again, try different learning rates and maximum number of iteration for the gradient descent method to observe the difference. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='4.2'></a>\n",
    "### 4.2 - Coverage rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For logistic regression models, we do not have a 'normal equation', and we should ues Taylor expansion to derive the variance of the estimators and use the Slutsky theorem to derive its limiting distribution. In order to obtain meaning results, we need to show consistency of our estimators. All of these materials will be covered in an advanced (statistics or) probability course, and check relevant materials for more details.  \n",
    "\n",
    "It can be show that the estimators obtained by both methods are consistent, but the convergence of the Newton-Raphson method is much faster. Thus, we only consider the estimators obtained by the Newton-Raphson method in this section.\n",
    "\n",
    "It can be shown that a variance estimator for the estimated parameter is (check by yourself if you are interested)\n",
    "$$\\left\\{n H(\\hat{\\tilde{w}})\\right\\}^{-1},\\tag{2}$$\n",
    "where $\\hat{\\tilde{w}}$ is the estimator obtained by the Newton-Raphson method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need a function to estimate the covariance matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE EXISTNG CODE\n",
    "def var_est(x, y, est_par):\n",
    "    # x: the generated feature vector of length n\n",
    "    # y: the generated feature vector of length n\n",
    "    # est_par: the estimated model parameter\n",
    "\n",
    "    # Step 1. Obtain the sample size\n",
    "    # Step 2. Augment the feature vector with the intercept term, check the relationship between \\tilde{x} and x\n",
    "    # Step 3. Obtain z using the estimated parameter\n",
    "    # Step 4. Obtain a using the estimated parameter\n",
    "    # Step 5. Obtain the Hessian matrix H\n",
    "    # Step 4. Obtain the estimated covariance matrix by (2) \n",
    "    \n",
    "    ### YOUR CODE BEGINS HERE (approximately 6 lines)\n",
    "    n = \n",
    "    aug_x = \n",
    "    z = \n",
    "    A = \n",
    "    H_J = \n",
    "    var_par = \n",
    "    ### YOUR CODE ENDS\n",
    "    \n",
    "    return var_par\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By law of large numbers, the variance estimator does not change too much when the sample size is large. Thus, if your function is correct, your variance estimator should approximately matches the sample variance. That is, the following two matrices should be almost the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "n = 1000\n",
    "x, y = train_data_generation(n=n, rn=i)\n",
    "est_par = NR_est(x,y)\n",
    "var_par = var_est(x, y, est_par)\n",
    "print('Your estimated covariance matrix is:')\n",
    "print(var_par)\n",
    "print('\\nThe sample covariance matrix is:')\n",
    "print(np.cov(np.transpose(par_NR_result)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be shown that under our model setup, the following asymptotic result holds\n",
    "$$\n",
    "nH(\\hat{\\tilde{w}})(\\hat{\\tilde{w}} - \\tilde{w}_0)^T\\to N((0,0,0)^T,I_3)\n",
    "$$\n",
    "in distribution as $n\\to\\infty$, where $I_3$ is a $3\\times 3$ identity matrix. By this result, we can also check the coverage rates of two-sided 95\\% confidence intervals \n",
    "$$\n",
    "(\\hat{b}_0 - q_{0.975}\\hat{\\sigma}_b,\\hat{b}_0 - q_{0.025}\\hat{\\sigma}_b)\\\\\n",
    "(\\hat{w}_0 - q_{0.975}\\hat{\\sigma}_{w0},\\hat{w}_{0} - q_{0.025}\\hat{\\sigma}_{w0}),\\\\\n",
    "(\\hat{w}_1 - q_{0.975}\\hat{\\sigma}_{w1},\\hat{w}_1 - q_{0.025}\\hat{\\sigma}_{w1}),\n",
    "$$\n",
    "where $q_\\alpha$ is the $\\alpha$-th quantile of a standard normal distribution, and $\\hat{\\sigma}_b$, $\\hat{\\sigma}_{w0}$ and $\\hat{\\sigma}_{w1}$ are the estimated standard error of the estimators. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your results are very close to the desired value, you can move to the coverage rates. For each simulated training dataset, we need to check whether the constructed confidence interval covers the true value or not. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE EXISTING CODE\n",
    "def cr_indicator(x,y,alpha):\n",
    "    # x: the generated feature vector of length n\n",
    "    # y: the generated feature vector of length n\n",
    "    # alpha: significance level. For example, alpha=0.05 corresponds to a 95% two-sided confidence interval.\n",
    "\n",
    "    # Step 1. Obtain (1-\\alpha/2)th quantile of a standard normal distribution.\n",
    "    # Step 2. Obtain the model estimator\n",
    "    # Step 3. Obtain the estimated covariance matrix for the estimated parameter\n",
    "    # Step 4. Obtain the estimated variances for every parameter\n",
    "    # Step 5. Generate the indicator whether the true value of b_0 lies in the confidence interval\n",
    "    # Step 6. Generate the indicator whether the true value of w_00 lies in the confidence interval\n",
    "    # Step 7. Generate the indicator whether the true value of w_00 lies in the confidence interval\n",
    "    \n",
    "    ### YOUR CODE BEGINS HERE (approximately 7 lines)\n",
    "    quan_normal = \n",
    "    \n",
    "    est_par =\n",
    "    var_par = \n",
    "    est_se = \n",
    "    \n",
    "    ind_b =\n",
    "    ind_w1 = \n",
    "    ind_w2 = \n",
    "    ### YOUR CODE ENDS    \n",
    "\n",
    "    \n",
    "    return np.array([ind_b, ind_w1, ind_w2])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "cr_est = np.ones((1000,3))-10 # Usually, we should have an NA matrix or a NULL matrix\n",
    "for i in range(1000):\n",
    "    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i\n",
    "    cr_est[i,:]  = cr_indicator(x, y, alpha = 0.05).flatten()\n",
    "print('Your coverage rates for the two parameters are:')    \n",
    "print(np.mean(cr_est, axis = 0)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**If your results are close to $1-\\alpha = 0.95$, then your code is good.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='5'></a>\n",
    "## 5 - Closer check for the two estimation methods\n",
    "We have mentioned that the convergence of the Newton-Raphson method is much faster than the gradient descent method, and we compare the computation efficiency of the two methods in this section. We still consider the case for $n=1\\,000$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we have a closer look at the estimators for both methods in each iteration. Fill in the following blank space to finish the following task:\n",
    "\n",
    "* Generate a training dataset with sample size $n=1\\,000$ with random seed to be 1234\n",
    "* Set initials to be (0,0,0) for the three parameters.\n",
    "* ONLY update your estimator for TEN iterations, and the learning rate for the gradient descent method is set to 0.01.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "n = 1000\n",
    "rn = 1234\n",
    "alpha = 10**(-2)\n",
    "iter_NR = np.ones((100,3))-10\n",
    "iter_GD = np.ones((100,3))-10\n",
    "\n",
    "x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i\n",
    "aug_x = np.concatenate((np.ones((n,1)),x), axis = 1)\n",
    "\n",
    "par_NR = np.zeros((3,1))\n",
    "par_GD = np.zeros((3,1))\n",
    "iter_index = 0\n",
    "iter_NR[iter_index,:] = par_NR[:,0]\n",
    "iter_GD[iter_index,:] = par_GD[:,0]\n",
    "\n",
    "while iter_index < 99:\n",
    "    iter_index += 1\n",
    "    z_NR = aug_x @ par_NR\n",
    "    A_NR = sigmoid(z_NR)        \n",
    "    nabla_J_NR = np.transpose(aug_x) @ (A_NR-y) / n\n",
    "    H_J_NR = np.transpose(aug_x*(A_NR*(1-A_NR))) @ aug_x /n\n",
    "    par_NR -= np.linalg.inv(H_J_NR) @ nabla_J_NR\n",
    "\n",
    "    z_GD = aug_x @ par_GD\n",
    "    A_GD = sigmoid(z_GD)        \n",
    "    nabla_J_GD = np.transpose(aug_x) @ (A_GD-y) / n\n",
    "    par_GD -= alpha * nabla_J_GD \n",
    "                          \n",
    "    iter_NR[iter_index,:] = par_NR[:,0]\n",
    "    iter_GD[iter_index,:] = par_GD[:,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We visualize the estimators for $b_0$ based on the two methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "plt.scatter(np.arange(100),iter_NR[:,0], c = 'red',  label='Newton-Raphson')\n",
    "plt.scatter(np.arange(100),iter_GD[:,0], c = 'blue', label='Gradient descent')\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.axhline(y=-0.5, color='black', linestyle='-')\n",
    "plt.ylim(-1,1)\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('Estimators for b0')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We conclude that the Newton-Raphson method converges very quickly, but it takes more iterations for the gradient descent method to converge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLEASE DO NOT CHANGE THE FOLLOWING CODE\n",
    "T1 = time.time()\n",
    "for i in range(100):\n",
    "    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i\n",
    "    par_est =  NR_est(x,y)\n",
    "T2 = time.time()\n",
    "print('The Newton-Raphson method takes about %s seconds' % ((T2 - T1)))\n",
    "\n",
    "T1 = time.time()\n",
    "for i in range(100):\n",
    "    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i\n",
    "    par_est = BGD_est(x,y)\n",
    "T2 = time.time()\n",
    "print('The (Batch) gradient descent method takes about %s seconds' % (T2 - T1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By comparison, we observe that the Newton-Raphson method is indeed much faster than the gradient descent method. Besides, the Newton-Raphson method generates more accurate estimators than the gradient descent method. <font color='red'>However, the gradient descent method is more commonly used than the Newton-Raphson method, especially in deep learning algorithms, since it is impossible to obtain the Hessian matrix when we have billions of parameters. </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='6'></a>\n",
    "## 6 - Additional homework\n",
    "When you finish the above basic function, we consider the following additional homeworks. \n",
    "\n",
    "* Usually, when we conduct simulation studies, we may consider more than one setup. For simplicity, we only consider different setups for the sample size. Consider $n\\in \\{500, 1\\,000, 1\\,500, 2\\,000\\}$ and use figures to report the following aspects. Plese notice that for all figures, sample sizes index the x-axis. \n",
    "   * In one figure, show the coverage rate for $b_0$ as the y-axis and comment your results.\n",
    "   * In one figure, show the coverage rate for $w_0$ as the y-axis and comment your results.\n",
    "   * In one figure, show the coverage rate for $w_1$ as the y-axis and comment your results.\n",
    "   * In one figure, show the computation efficiency of the two estimation methods, and different colors should be used for different methods. Comment your results. \n",
    "   * For the case $n=1\\,000$, use a figure to show the estimated parameters for $(w_{00},w_{01})$ in 10 iterations. You may set $x$-axis as the estimator for $w_{00}$ and $y$-axis for that of $w_{01}$. Then, use arrows to show how they are updated by different methods. Please use different colors to denote different methods."
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}