{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lecture 34: A Look Ahead; Examples of Regression Example, Sampling from a Finite Population\n",
"\n",
"\n",
"## Stat 110, Prof. Joe Blitzstein, Harvard University\n",
"\n",
"----"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Top-10 List\n",
"\n",
"1. Conditioning ... is the soul of statistics!\n",
"1. Symmetry ... is powerful but dangerous\n",
"1. Random variables and their distributions\n",
"1. Stories (proofs, backgrounds of the distributions covered)\n",
"1. Linearity\n",
"1. Indicator random variables\n",
"1. LOTUS\n",
"1. Law of Large Numbers\n",
"1. Central Limit Theorem\n",
"1. Markov Chains\n",
"\n",
"Items 1 through 4 deal with the Big Picture<sup>™</sup> questions: _What is randomness? How do we think about uncertainty?_\n",
"\n",
"Items 5 through 7 are for computing expected values (mean, variance & standard deviation).\n",
"\n",
"Items 8 through 10 are important for understanding long-run behavior."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Where to go from here?\n",
"\n",
"Some topics to study from here on out:\n",
"\n",
"* Statistical inference (we have data, need to estimate parameters or make predictions)\n",
"* Regress & linear models\n",
"* Finance\n",
"* Computational biology\n",
"* Stochastic processes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Advice\n",
"\n",
"* Learn R\n",
"* Learn C\n",
"* Read Mostly Harmless Econometrics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Ex. A Simple Linear Regression\n",
"\n",
"You've seen this before:\n",
"\n",
"\\begin{align}\n",
" Y &= \\beta_0 + \\beta_1 \\, X + \\epsilon\n",
"\\end{align}\n",
"\n",
"* We want to use $X$ to predict $Y$\n",
"* $\\beta_j$ are linear coeffiecients, with $\\beta_0$ being the value of $Y$ when $x=0$ (default value)\n",
"* $\\epsilon$ error term (since $X$ is not perfect)\n",
"* a common assumption is $\\mathbb{E}(\\epsilon | X) = 0$ (centered at 0, $\\epsilon$'s distribution may or may not be normal)\n",
"\n",
"So how would we solve for $\\beta_1$?\n",
"\n",
"We can start by treating $Cov$ as an _operator_!\n",
"\n",
"\\begin{align}\n",
" Cov(Y, X) &= Cov\\left( (\\beta_0 + \\beta_1 \\, X + \\epsilon), X \\right) \\\\\n",
" &= Cov(\\beta_0, X) + Cov\\left( (\\beta_1 \\, X), X\\right) + Cov(\\epsilon, X) \\\\\n",
" \\\\\n",
" \\text{now } Cov(\\beta_0, X) &= 0 &\\quad \\text{ since } Cov \\text{ of constant with anything is } 0 \\\\\n",
" \\\\\n",
" \\text{and } Cov\\left( (\\beta_1 \\, X), X\\right) &= \\beta_1 \\, Cov(X, X) &\\quad \\text{by definition of }Var \\\\\n",
" &= \\beta_1 \\, Var(X) \\\\\n",
" \\\\\n",
" \\text{and since } \\mathbb{E}(\\epsilon) &= \\mathbb{E}\\left( \\mathbb{E}(\\epsilon|X) \\right) = \\mathbb{E}(0) = 0 \\\\\n",
" \\text{and further } \\mathbb{E}(\\epsilon \\, X) &= \\mathbb{E}\\left( \\mathbb{E}(\\epsilon \\, X | X) \\right) &\\quad \\text{ by Adam's Law} \\\\\n",
" &= \\mathbb{E}\\left( X \\mathbb{E}(\\epsilon | X) \\right) &\\quad \\text{ since } X \\text{ is known, we can pull it out} \\\\\n",
" &= \\mathbb{E}(0) \\\\\n",
" &= 0 \\\\\n",
" \\text{so }Cov(\\epsilon, X) &= \\mathbb{E}(\\epsilon \\, X) - \\mathbb{E}(\\epsilon) \\, \\mathbb{E}(X) \\\\\n",
" &= 0 - 0 = 0 \\\\\\\\\n",
" \\Rightarrow \\beta_1 &= \\frac{Cov(X,Y)}{Var(X)} &\\quad \\text{(population version)} \n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculate $\\beta_1$ with $Cov(X,Y)$ and $Var(X)$\n",
"\n",
"Here we calculate $\\beta_1 = \\frac{Cov(X,Y)}{Var(X)} \\,$ using [numpy.cov](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.cov.html):"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"0.64383561643835618"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"\n",
"X = np.array([95, 85, 80, 70, 60])\n",
"Y = np.array([85, 95, 70, 65, 70])\n",
"\n",
"# numpy.cov(X, Y) returns the matrix\n",
"# [ Cov(X,X), Cov(X,Y)]\n",
"# [ Cov(X,Y), Cov(Y,Y)]\n",
"covM = np.cov(X,Y)\n",
"\n",
"beta_1 = covM[0,1]/covM[0,0]\n",
"beta_1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculate $\\beta_1$ via sklearn LinearRegression API\n",
"\n",
"For comparison's sake, we also obtain $\\beta_1$ via [`sklearn.linear_model.LinearRegression`](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"0.64383561643835607"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sklearn import linear_model\n",
"\n",
"regr = linear_model.LinearRegression()\n",
"regr.fit(np.matrix(X).T, Y).coef_[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Ex. Sampling from a Finite Population\n",
"\n",
"Here's the set-up:\n",
"\n",
"* We have a finite population , size $N$.\n",
"* Let $Y_1, Y_2, ..., Y_N$ be some value of interest (height, weight, opinion.\n",
"* Each person in the population can be uniquely identified.\n",
"* $Y_j$ are fixed, non-random values, but nevertheless _unknown_.\n",
"* Using some sampling scheme to obtain a sample of size $n$.\n",
"* Using this sample, we want to infer to the sum of $Y$ (or perhaps the average).\n",
"* Assume that the _inclusion probability_ of person $j$ ending up in our sample is $p_j$ (assume that the true value is known). \n",
"* Our sample data takes the form $(X_1, Z_1), (X_2, Z_2), ..., (X_n, Z_n)$, where\n",
" * $Z_j$ is the ID of the $j^{th}$ person in our sample\n",
" * $X_j = Y_j$\n",
"\n",
"### The difference between $X_j$ and $Y_j$\n",
"\n",
"It is important to understand the difference between $X_j$ and $Y_j$:\n",
"\n",
"* $Y_j$ is a **_fixed, non-random value_**.\n",
"* $X_j$, due to our random sampling (person $j$ was randomly selected from the population with probability $p_j$), is a **_random variable_**.\n",
"\n",
"### How do we get an unbiased estimator for the total?\n",
"\n",
"Let $t_y$ be the true population total $\\sum_{1}^{N} Y_i$. How can we use random sampling of this finite population to find $\\hat{t_y}$? \n",
"\n",
"The claim is that $\\sum_{j=1}^{n} \\frac{X_j}{P_{Z_j}}$ is an unbiased estimator for $t_y$ we are looking for.\n",
"\n",
"\\begin{align}\n",
" t_y &= \\sum_{j=1}^{n} \\frac{X_j}{P_{Z_j}} \\\\\n",
" &= \\sum_{j=1}^{N} \\frac{I_j \\, Y_j}{P_j} &\\quad \\text{where } I_j = 1 \\text{ if person } j \\text{ included in sample}\\\\\\\\\n",
" \\mathbb{E}(t_y) &= \\mathbb{E}\\left( \\sum_{j=1}^{N} \\frac{I_j \\, Y_j}{P_j} \\right) &\\quad \\text{ find expected value to get } \\hat{t_y} \\\\\n",
" &= \\sum_{j=1}^{N} \\frac{P_j \\, Y_j}{P_j} &\\quad \\text{ by linearity} \\\\\n",
" &= \\boxed{ \\sum_{j=1}^{N} Y_j }\n",
"\\end{align}\n",
"\n",
"This is known as the [Horvitz-Thompson Estimator](https://en.wikipedia.org/wiki/Horvitz%E2%80%93Thompson_estimator), or alternately inverse probability weighting.\n",
"\n",
"_But is an unbiased estimator good?_"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Basu's Circus Elephants\n",
"\n",
"Statistics is not easy, and it requires a lot of effort to keep your eyes open and question whether or not a tentative method is really going to yield a proper answer. Here is an anecdote to illustrate an example of when blindly applying an Horvitz-Thompson estimate ends in disaster.\n",
"\n",
"> The circus owner is planning to ship his 50 adult elephants and so he needs a rough estimate of the total\n",
"> weight of the elephants. As weighing an elephant is a cumbersome process, the owner wants to estimate\n",
"> the total weight by weighing just one elephant. Which elephant should he weigh?<p/>So the owner looks\n",
"> back on his records and discovers a list of the elephants' weights taken 3 years ago. He finds that 3 years\n",
"> ago Sambo the middle-sized elephant was the average (in weight) elephant in his herd. He checks with\n",
"> the elephant trainer who reassures him (the owner) that Sambo may still be considered to be the average\n",
"> elephant in the herd. Therefore, the owner plans to weigh Sambo and take 50y (where y is the present\n",
"> weight of Sambo) as an estimate of the total weight of the 50 elephants.<p/>But the circus statistician is\n",
"> horrified when he learns of the owner's proposed sampling plan. \"How can you get an unbiased estimate\n",
"> of Y this way?\", protests the statistician.<p/>So, together they work out a compromise sampling plan. With\n",
"> the help of a table of random numbers they devise a plan that allots a selection probability of 99/100 to\n",
"> Sambo and equal selection probabilities of 1/4900 to each of the other 49 elephants. Naturally, Sambo is\n",
"> selected and the owner is happy. <p/>\"How are you going to estimate Y?\", asks the statistician.<p/>\"Why? The\n",
"> estimate ought to be 50y of course,\" says the owner.<p/>\"Oh! No! That cannot possibly be right,\" says the\n",
"> statistician, \"I recently read an article in the Annals of Mathematical Statistics where it is proved that\n",
"> the Horvitz-Thompson estimator is the unique hyperadmissible estimator in the class of all generalized\n",
"> polynomial unbiased estimators.\"<p/>\"What is the Horvitz-Thompson estimate in this case?\" asks the owner,\n",
"> duly impressed.<p/>\"Since the selection probability for Sambo in our plan was 99/100,\" says the statistician,\n",
"> \"the proper estimate of Y is 100y/99 and not 50y.\"<p/>\"And, how would you have estimated Y\", inquires\n",
"> the incredulous owner, \"if our sampling plan made us select, say, the big elephant Jumbo?\"<p/>\"According to\n",
"> what I understand of the Horvitz-Thompson estimation method,\" says the unhappy statistician, \"the proper\n",
"> estimate of Y would then have been 4900y, where y is Jumbo's weight\".<p/>That is how the statistician lost\n",
"> his circus job (and perhaps became a teacher of statistics)\n",
"\n",
"----"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"View [Lecture 34: A Look Ahead | Statistics 110](http://bit.ly/2Npa7ze) on YouTube."
]
}
],
"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.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}