Statistics course weekly lab

{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Grade: /100 pts\n", "\n", "# Assignment 04: Confidence Intervals & The Bootstrap\n", "\n", "Once you are finished, ensure to complete the following steps.\n", "1. Restart your kernel by clicking 'Kernel' > 'Restart & Run All'.\n", "2. Fix any errors which result from this.\n", "3. Repeat steps 1. and 2. until your notebook runs without errors.\n", "4. Submit your completed notebook to OWL by the deadline." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Add the necessary imports for this homework \n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import sklearn.model_selection\n", "import sklearn.linear_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 1: /10pts \n", "In this question, you will construct a confidence interval for the sample mean, not using the normal distribution, but the t-distribution (see end of lecture 4.3), which is more accurate for small sample sizes. \n", "\n", "The $100(1-\\alpha)\\%$ confidence interval is \n", "\n", "$$ \\bar{x} \\pm t_{1-\\alpha/2, n-1} \\dfrac{\\hat{\\sigma}}{\\sqrt{n}} $$\n", "\n", "Where $ t_{1-\\alpha/2, n-1}$ is the appropiorate quantile of a Student's t distribution with $n-1$ degrees of freedom. \n", "Write a function called `confidence_interval` which takes as it's argument an array of data called `data` and returns two things:\n", "\n", "* An estimated mean of `data`, and \n", "\n", "* The lower and upper bounds of the 95% confidence interval for the mean of `data`. Ensure these are returned in a numpy array of shape (2,)\n", "\n", "To get the appropirate quantiles for the t-distribution, you can use `scipy.stats.t`, which implements some statistical functions for the t-distribution. Take a look at the documentation for `scipy.stats.t`, especially the `ppf` method.\n", "\n", "Here is the documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.t.html\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def confidence_interval(data):\n", "\n", " # Note, np.std divides by n and not n-1\n", " # Force it to apply the correct formula by ussing ddof=1\n", " # Alternaively, you can use scipy.stats.sem to compute\n", " #The standard error\n", " \n", " return estimated_mean, bounds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 2: /15pts\n", "\n", "The \"95% confidence interval\" is named so because the long term relative frequency of these estimators containing the true estimand is 95%. That is to say **if I construct 100 95% confidence intervals for the sample mean again and again from the same data generating mechanism, 95 of these intervals I construct will contain the true population mean**.\n", "\n", "Write a function called `ci_simulation` that runs some simulations to show this is the case. From a standard normal distirbution, sample 25 observations and construct a confidence interval. Do this 20 times and plot the intervals using `matplotlib.pyplot.errorbar`. Color the bar red if the confidence interval does not capture the true mean and blue if it does. If you are unfamilliar with `matplotlib.pyplot.errorbar`, I highly suggest reading Matplotlib's excellent documentation which has some examples at the bottom of the webpage.\n", "\n", "If you are unfamilliar with how to sample random numbers, I suggest you look at `numpy.random.normal`. Try searching for the documentation for that function yourself if you need to.\n", "\n", "Here is the documentation for `matplotlib.pyplot.errorbar`: https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.errorbar.html" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def ci_simulation():\n", " # Set the random seed to always get the same random numbers. \n", " # This is for Reproducibility. \n", " np.random.seed(4)\n", " \n", " # Create the figure.\n", " fig, ax = plt.subplots(dpi = 120)\n", "\n", " # If the interval crosses this line, it should be blue, else red.\n", " ax.axhline(0, color = 'k')\n", "\n", " # Do the following 20 times\n", " for i in range(20):\n", "\n", " #Draw 25 observations from a standard normal\n", "\n", " # Compute what we need for the CI, namely the mean and the bounds\n", "\n", " \n", " # color should be blue if it crosses the black line\n", " color = 'blue'\n", " if (min(bounds)>0)|(max(bounds)<0):\n", " # but in the case it does not, turn it red\n", " color = 'red'\n", "\n", " # Need to get the length of the interval from bounds\n", " interval_len = 1.0/2*(bounds[1] - bounds[0])\n", " ax.errorbar(i, mu, yerr=interval_len, color = color, fmt = 'o')\n", "\n", " # This function does not have to return anything\n", " return None\n", "\n", "ci_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 3: /8pts\n", "\n", "If you haven't changed the random seed from 4 and if you implemented the solution correctly, you should one red interval.\n", "\n", "Answer the following below in no more than 3 sentences:\n", "\n", "a) How many red intervals did we expect to see? What is your justifiation for this?\n", "\n", "Changing the random seed might affect how many red intervals you see. Try changing the random seed in your function to 3. This will yield two red intervals (which is different than what you should expect to see). \n", "\n", "b) Why does the simulation sometimes deviate from the predicted results?\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "a) \n", "\n", "b) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 4: /10pts\n", "\n", "Load in the `hockey_draftees_2005.csv` data into pandas. It contains data from hockey players drafted in 2005, including their rank, weight (wt - in pounds) and height (ht - in inches). \n", "\n", "Fit a linear model of weight (`wt`) explained by height (`ht`) using a linear regression model from sklearn, as done in the lab. Call your fitted model `model`. \n", "Make a scatter plot of the height (x-axis) against weight (y-axis). \n", "Add the predicted values for 66-80 inches. \n", "\n", "Calculate the residuals from the fit, and report the r-squared for this model." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df=\n", "\n", "# Make it and fit the model \n", "model =\n", "\n", "# Make the scatter plot \n", "\n", "\n", "#Generate and plot the predicted values\n", "\n", "# Calculate residuals, R2 and print it \n", "\n", "print('Training rsquared is ',model.rsquared)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 5: /15pts\n", "### Bootstrap confidence intervals on parameters\n", "\n", "How confident can we be about the relation between height and weight? \n", "To judge this we need confidence intervals let's use the bootstrap.\n", "\n", "Modify the function `BootstrapCoef` from lab 04 - part 2 to conduct a boostrap analysis for this regression model; \n", "\n", "* `data`, which is a dataframe having columns 'weight' and 'height'\n", "* `numboot` which is an integer denoting how many bootstrap replications to perform.\n", "\n", "The function should return `theta`, a numpy array of regression coefficients of size (numboot, 2)\n", "\n", "You can use `pd.DataFrame.sample` with `replace = True` to perform the resampling. `bootstrap` should return:\n", "params: a numpy array of size [numboot,numParams] of bootstraped parameter values. The parameters are the intercept value and the slope from the linear regression. \n", "Tip: Note that the intercept can be retrieved from model.intercept_, whereas all the other regression coefficients are stored in model.coef_. \n", "\n", "Here is the documentation for `pd.DataFrame.sample`: https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.sample.html\n", "\n", "Then call the function to get 100 boostrap samples for your linear regression model of wt explained by height. \n", "Make a joint scatter plot of the parameter value for the intercept and for the slope. \n", "Written answer: What do you notice? Why do you think the estimate for the intercept has such a high negative correlation with the slope? " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Write a Bootstrap function that records the fitted models \n", "def BootstrapCoef(data,numboot=1000):\n", "\n", " \n", " \n", " \n", " return theta\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Intercept and slope regressors are highly colinear - as we did not subtract the mean of ht before thre regression. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 6: /12pts\n", "\n", "Plot the bootstrap estimates for the slope as a histogram. Use your samples to compute a 95% confidence interval. Note that the CI should be constructed around the sample estimate of the slope. How can you interpret this confindence interval? " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculate bootstrap interval\n", "\n", "\n", "\n", "print('My confidence interval is between', boot_ci[0], ' and ', boot_ci[1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Written answer: The interval contains the true slope parameter with a probability of 95%. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 7: /15pts\n", "Modify the function `BootstrapPred` from lab04 to bootstrap your fit and generate a predict from each of these bootstrapped models. \n", "\n", "Draw again a scatter plot of height against weight and plot the predictions from the 20 fitted bootstrap models for the height ranging from 60 to 88. \n", "\n", "Written answer: Where are we most uncertain in our prediction about the weight of a player? How does the negative correlation between slope and intercept play a role here? Why is the spread of the prediction in the mean weight so much lower than the variability of our intercept parameter? " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The greatest uncertainty we have for the very small and very tall players. All lines agree in their prediction for the players of middle weight. For the lines to cross in the middle of the graph, a large slope needs to have a small intercept and a small slope needs to have a large intercept. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 8: /15pts\n", "\n", "Now, let's see how well our model performs out of sample. Load in the `hockey_draftees_test.csv` file into a dataframe. \n", "Use your fitted `model` to make predictions. \n", "\n", "Make a scatter plot of the test data and superimpose the prediction of the model. \n", "To evaluate this prediction, calculate the r-square value for the out of sample (oos) data. Statsmodels doesn't provide a function to compute r-squared on new data. You will have to write one yourself or find one that performs the computation for you." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df_test=pd.read_csv('hockey_draftees_test.csv')\n", "\n", "# Make the scatter plot \n", "\n", "# Generate and plot the predicted values\n", "\n", "# Now do the prediction for the test data and compute R2\n", "\n", "\n", "print('Out of sample rsquared is ', rsquared_oos)\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }