exercises.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "generate_data <- function(N,beta,sd){\n",
    "    x<- rnorm(N,0,1)\n",
    "    eps<- rnorm(N,0,sd)\n",
    "    y<- beta*x + eps\n",
    "    list(x=x,y=y)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_regression <- function(y,x){\n",
    "    coef <- calc_coef(y,x)\n",
    "    se <- calc_se(y,x,coef)\n",
    "    list(coef=coef,se=se)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc_coef <- function(y,x) cov(x,y)/(var(x))\n",
    "\n",
    "calc_se <- function(y,x,coef){\n",
    "    n <- length(y)\n",
    "    eps <- y-x*coef\n",
    "    e_sd <- mean(eps^2)\n",
    "    se <- sqrt(e_sd /(n*var(x)))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "eval_model <- function(coef, se, beta, conf = 1.96){\n",
    "    up <- coef + se*conf\n",
    "    down <- coef - se*conf\n",
    "    (beta > down) & (beta < up)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate <- function(N, beta, sd){\n",
    "    d <- generate_data(N,beta,sd)\n",
    "    m <- run_regression(d$y, d$x)\n",
    "    eval_model(m$coef, m$se, beta)\n",
    "}\n",
    "# 1:M gives list of int from 1 - M, list is filled with True /False from \n",
    "# eval_model at each element\n",
    "avg_simulations <- function(M,N,beta,sd){\n",
    "    inside <- sapply(1:M, function(x){\n",
    "        simulate(N,beta,sd)\n",
    "    })\n",
    "    sum(inside)/M\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAC5VBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8hISEiIiIjIyMkJCQmJiYnJycoKCgpKSkq\nKiorKyssLCwtLS0uLi4vLy8xMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9\nPT0/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09Q\nUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFi\nYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0\ndHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWG\nhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaYmJiZ\nmZmampqbm5ucnJydnZ2enp6fn5+goKCjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKyt\nra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi6urq7u7u8vLy9vb2+vr6/v7/A\nwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDS0tLT\n09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl\n5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb3\n9/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///8VPnEvAAAACXBIWXMAABJ0AAASdAHeZh94AAAg\nAElEQVR4nO3de5zddX3n8QNUEBFdUdcLZRdv1VoWWazX6tYuuDohNgQwV8IGugKSQK0BQQmm\nJJRbhCAtVMCCUKTaYF0CxF1cIYESQ0JpgIQQJgnJx8Rkcs9cfn/vuQzJiXmfmfP+5Uvmd+a8\n3n+Mk3N+55XvI5ync82klDHGDniloT4AY8NhQGIswYDEWIIBibEEAxJjCQYkxhIMSIwlGJAY\nS7AEkLo2mtuy3X1Ew23uTdjamSy1qTdha3ey1MbehK3uhKmUrXSpXb2DXrIpJaRNYa5ru/uI\nhvttlrC1M1lqQ5aw1Z0sFVnCVm/CVMpWulR3NuglG4EkWkByBqQAkm4ByRmQAki6BSRnQAog\n6RaQnAEpgKRbQHIGpACSbgHJGZACSLoFJGdACiDpFpCcASmApFtAcgakAJJuAckZkAJIugUk\nZ0AKIOkWkJwBKYCkW0ByBqQAkm4ByRmQAki6BSRnQAog6RaQnAEpgKRbQHIGpHAgbblh4pir\n1tVeXztz3OhrNmXZ1zvKGw2kxgOSl2oDSDOmrVh93QW9lVd3/8WMzpVXXJZlkx6MynMFSA0H\nJC81/CHFiOXlt0ojF1deX9bxm/INHSuz05/a5xog7TcgeanhD+nxUX3llxfeV3l9aUdXlvWM\nfGR3x81TzpnZCaTGA5KXGv6QHjq78vLy2yovt4/92+7uu0f+ZNP4G5ctmz5+a/mm1Q+U9/IW\nczt2u49ouO1Zuta27mSprVnCVm+y1JasJ12rL2EqZStdqicb9JKtTUOatBdStuS8kWfefd7c\n6uvbR88rv5x/UnkLB04wNmzXu+e1QSAtqL1rd3//L7d2d3/lidqr59+T8Rap0XiL5KWG/1uk\nDSNeyLLNpy2tvN7zWPk9widHbl45pzvLdoye/9o1fIy03/gYyUsN/4+RsllTV3ROv6Qvm1d+\nj27KzFg68dasa8zstZ0zJ+0EUsMByUu1AaRtsyeMnVm+/Noryu/Ifev0cbeX3xotv+LMcTNe\n3XMJkPYbkLxUG0BqYkDab0DyUkACkhyQvBSQgCTXkpBWz//5CqMFpACSbrU3pHveWyodNaP5\nFpACSLrV1pB+dWSpstubbgEpgKRbbQ3pnKqj0h813QJSAEm32hrSF2qQ3tp0C0gBJN1qa0hj\na5A+3HQLSAEk3WprSPOOqEK6rukWkAJIutXWkOLmo0ulw89f33QLSAEk3WpvSLHs7h8sMlpA\nCiDpVptDMgekAJJuAckZkAJIugUkZ84zdvUvlwxw77rnn193oKfZMyDlGpCsDRGk7xxVKn3s\n/za694fHlkrH/jDFkSoDUq4BydrQQLqu+on19z6v76193v2Ihw/6sQYdkHK2gOSs+Wfsu2pf\n671S3/vl2r1fPujHGnRAytkCkrOmn7Ev16SUJui7P+J+S0WiYw0+IOVsAclZ08/Y9UfVqFyi\n7/6T2r1/crCPNfiAlLMFJGfNP2PPrUo58nF9720l8+9vpDrWoANSzhaQnDX/jH2l8p3lR/9d\no7svKt97+JQkZwog5RyQrA3V15F+fs1t/9b43gW3377ggI/z2oCUa0Cyxnc2OANSzhaQnAEp\ngKRbCZ/8qzekawHJaqVLASlnKxmkG99ROvTzTyaKAclrpUsBKWcrFaRbq5/PPf6lNDUgea10\nKSDlbKWC9Pu1r4z8dZoakLxWuhSQcrYSQVo18LfEuAOS10qXAlLOViJI695YgzQ1TQ5IXitd\nCkg5W6netRtXdXREw7+A4w1IXitdCkg5W6kgvfTxiqMbEtWA5LXSpYCUs5Xs09/r751x89Op\nYkDyWulSQMrZ4jsbnAEpgKRbQHIGpACSbgHJWVEgvfw7rQM6yT4DUs4WkJwVAtKqi95Wetd3\n1ta3Epynf0DK2QKSs0JAOn2/r9gBKdeAZG24QXq09qXvQ5+tayU5UXVAytkCkrMiQJrT/91Y\nD9S1kpyoOiDlbAHJWREg3dUPqe6nSwIp14BkbbhBevHtVUfvf7WuleRE1QEpZwtIzooAKe57\nc9nR2x+tbyU4T/+AlLMFJGeFgBRLvnveNS/u0zrg0+wZkHK2gOSsGJD2b6VLASlnC0jOgBRA\n0i0gOQNSAEm3Wh/SHSccftzUlxvf/3pBWnfjhw4/fvqa3CkgAUluiCB9v/rJ4FPWN7zg9YJ0\nefU3npQ7BSQgyQ0NpDVvq3158t6GV7xOkJa9ofYbP5Y3BSQgyQ0NpCf7v87/rYZXvE6QftL/\nG9+SNwUkIMkNDaQl/c/nGQ2veJ0gPdT/G9+RNwUkIMkN0cdIJ9Z+fNETDS94nSCtObb6G791\nWd4UkIAkN0SQfnVM5fk8q/EFr9cnG+ZW/nXLI+7MnQISkOSG6tPfL3znjIvmD3D/6/Z1pCXT\nzrj4yfwpIAFJji/IeikgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSOwBI\na6765IfPWFDfApLVSpcCUs5WESCt/7PK38078v/UtYBktdKlgJSzVQRIf1/76Qcn1bWAZLXS\npYCUs1UESJNrkA7p3NsCktVKlwJSzlYRIJ3b/485rt7bApLVSpcCUs5WESD9sAbp03UtIFmt\ndCkg5WwVAVJ8peLo6LqfqAUkr5UuBaScrUJAWnfzFz9x3jP1LSBZrXQpIOVsFQLS/i0gWa10\nKSDlbAHJGZACSLoFJGdACiDpFpCcASmApFtAcgakAJJuAckZkAJIugUkZ0AKIOkWkJwBKYCk\nW0ByBqQAkm4ByRmQAki6BSRnQAog6RaQnAEpgKRbQHIGpACSbgHJGZACSJWtv/WMjumr6lut\nAGn+Oaf8rydzt9oN0pIpp07455yhdd8fe/o+zw81IMX6Uyt/g+59L9S1WgDSnMqhj/jHvK02\ng/Tomyt/XJfn6vQ/P14c+Cogxc21v9M9rq5VfEj/9qbqod/Z2fjyAddmkD5c+2/8WJ7OTbXH\njh/4KiDFiNof1LvrWsWH1P/Dukr/krPVXpB+3f+n9d08nY7aY98z8FVAilNrf1DH1LWKD+nW\n/qfGT3K22gvSk/1/Wrnetzul9ti3D3wVkOLy2h/Uf69rFR/SgtqhjxjkPfeGay9Ir76j9sf1\n0zyd/ufHKQNfBaRY9aHKn9NRC+taxYcUF1b/687M22ovSHFn9U/rz3N1Vn2w8tg3Lxz4KiBF\nLJv8vmNHPF7fagFI62488Z0fvyN3q80gxY8/+84/unJ1o0sH3rLJHzhuxBODXAQk1WoBSAe4\ndoN0YOMLsjlbQHIGpACSbgHJGZACSLoFJGdACiDpFpCcASmApFtAcgakSAxpyyZz23a6j2i4\nLVnC1u5kqa4sYasnWWpTlrDVmzCVspUu1ZMNeklXSkg73O3uth/SaLuyhK2eZKmdWcJWb7LU\njixhqy9hKmUrXaq3iad2Ski8a7ffeNfOS/GuHZDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcC\nEpDkgOSlgDTEkFbddPFfP5UqdtAg/ePUqfcbLSBZrXSp9oG06LjKXxf9fqLaQYK0rvrTIk5b\n13QLSFYrXap9IH2q+jcg3/R0mtpBgjSr9neYZzXdApLVSpdqG0jP9v9wi5lpcgcJ0idqh/7j\npltAslrpUm0DaWE/pMvS5A4SpD+sHfojTbeAZLXSpdoG0uqja8/JH6bJHSRIo2uHHt10C0hW\nK12qbSDFNdWn5Gea/rh94B0kSE9Vf47u0U1/shFIXitdqn0grb/hP5XeMvH5RLWD9envhz/x\nhjd88pHmW0CyWulS7QOp/OTf0opfkF3t/IQoIHmtdKm2gsR3NjgDkjMg5WwByRmQAki6Newh\n/fSa7/0y72OX3njZnWvrbwBSAEm3hjmktV8slUqHX5rvwXdVvtTw4aV1twApgKRbwxzStNoX\nsIzvl927Z2pfsvvTupuAFEDSrWEO6fgapDPyPPZvao895Lm9NwEpgKRbwxzS22oYvpDnsZf1\nfzdW3b/OAKQAkm4Nc0gn1yycn+exd9Qee+Qre28CUgBJt4Y5pH+q/VufS/I8ds3Hqg/+dt1N\nQAog6dYwhxR3HlcqnTw/32OXjDis9JZv139XI5ACSLo13CFFvPqb/I/tXLx+n18DKYCkW8Mf\nEt/Z4AxIOVtAcgakAJJuAckZkAJIugUkZ0AKIOkWkJwBKYCkW0ByBqQAkm4ByRmQAki6BSRn\nQAog6RaQnAEpgKRbQHIGpACSbgHJGZACSLoFJGdACiDpFpCcASmApFtAcjYMIC26496Bfkwv\nkHK2gOSs5SGtP//wUumttzS+AEg5W0By1vKQrq3+pd8jHm54AZBytoDkrOUh/UHt51CMa3gB\nkHK2gOSs5SH9h9J+P6tv3wEpZwtIzloe0kdrkM5peAGQcraA5KzlIf1d7Z/y/lXDC4CUswUk\nZy0PKa58U6n03nsa3w+knC0gOWt9SLH8Xx4Z6J98A1LOFpCcDQNIgwxIOVtAcgakAJJuOU/+\nJ+/9f+sb3wskLwWkNoW07M9KpdIfL2p4P5C8FJDaFNIp1c+cnrim0f1A8lJAak9IC/r/taAH\nGl0AJC8FpPaE9ON+SDc1ugBIXgpI7Qnp8dIg/xwrkLwUkNoTUnyu6ugPG341D0heCkhtCunZ\nT5cdnbCg4f1A8lJAalNIEb+4Y966xvcCyUsBqW0hDTwgeSkgAUkOSF4KSECSA5KXGqi17l8f\nedlpHehh9g5IOVtAcnaQID38kVLpjZcO8H2Nv9tKcJ7+ASlnC0jODg6kf39X9SsNs5pvpThQ\nbUDK2QKSs4MD6du1L33/x+ZbKQ5UG5BytoDk7OBAmtD/TSRNf5gEpFwDkrXWg/SXNUdHN/1B\nEpByDUjWWg/SwjdVIV3YfCvFgWoDUs4WkJwdpM/a3VH5OY4dA/2Mkt9pJThP/4CUswUkZwfr\n60jP33XTL5zWgR5m74CUswUkZ3xnQwBJt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3\ngOQMSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkD\nUgBJt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFEDS\nLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jO\ngBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIIUD\nacsNE8dcta72+tqZ40Zfs2nf24CkBiQv1QaQZkxbsfq6C3orr+7+ixmdK6+4bJ/bgCQHJC81\n/CHFiOXlt0AjF1deX9bxm/INHSvrbwOSHJC81PCH9PiovvLLC++rvL60oyvLekY+Un8bkOSA\n5KWGP6SHzq68vPy2ysvtY/+2u/vukT+pu23+SeUtHDjB2LDdno9vBoU0qfKyBilbct7IM+8+\nb27dbYu/Vt6S3eZ6et1HNFx3lrBVzGPt7kuXyhK2UqYKeqxs0Et2NQtpQe3duPv7f7m1u/sr\nT/zObbxrt/94185LDf937TaMeCHLNp+2tPJ6z2Plhz05cnP9bUCSA5KXGv6QsllTV3ROv6Qv\nmzc3y6bMjKUTb917G5AaDUheqg0gbZs9YezM8uXXXpFlq791+rjbu/feBqRGA5KXagNITQxI\n+w1IXgpIQJIDkpcCUqEhLTrrgydMWdF8q/7J/+q1H3/fqT/PeywgeSkgFRnSoreWyvtoZ9Ot\n+if/VyuPLd2T81hA8lJAKjKkL1ctlL7TdKvuyf+z2mPf/Wq+YwHJSwGpyJDeWcNwStOtuif/\nt2uPLT2Z71hA8lJAKjKkd9csfKnpVt2T/6p+SIvyHQtIXgpIRYZ0Vs3C9U236p78v6g99gPr\n8x0LSF4KSEWGtOz3KxY+t67pVv2T/+LKY9/4v3MeC0heCkhFhhTLp33+S9evbb61z5P/R6M+\nO/npvMcCkpcCUqEhuS2+IOsMSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIz\nIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk\n3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJt4Dk\nDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IA\nSbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0g\nOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAU\nQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQL\nSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMg\nBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTd\nApIzIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQM\nSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJ\nt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFEDSLSA5\nA1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jOgBSJ\nIe3uMdfb5z6icSpL2EqX6kl4rJ6CHitLmErZSpfqG7zVnRISb5H2G2+RvBRvkYAkByQvBSQg\nyQHJSwEJSHJA8lJAApIckLwUkIAkByQvBSQgyQHJSwEJSHJA8lJAApIckLwUkIAkByQvBSQg\nyQHJSwEJSHJA8lJAApIckLwUkIAkByQvBSQgyQHJSwEJSHJA8lJAApIckLwUkIAkByQvBSQg\nyQHJSwEJSHJA8lJAApIckLwUkIAkByQvBSQgyQHJSwEJSHJA8lJAApIckLwUkIAkByQvBSQg\nyQHJSw03SLuBlGZA8lLDDdI7Lnq6eJA6B7oTSNaA5Cw/pM8dWvrotWuKBOmVi99ees9Vaxve\nDyRrQHJ2AB8jrZ3z2UMOO/VH2wsD6cxSZVMb3g8ka0BydmCfbFj9vZNLb5m8sBiQflF1VDrs\n2UYXAMkakJwd6Gftlo4pP3c/9VQRIN1Sg1R6oNEFQLIGJGcHBOnVG04oHfalB+b+18MeKgCk\nu/ohzWt0AZCsAclZfki7ftzxe6UPzap8umH3qe8vAKQXj6k6Or7hZxuAZA1IzvJDOqb05km/\n7H/9p4cUAFLce1TZ0TEPN7wfSNaA5Cw/pM/8YOue11f9oAiQ4pkrJ1/9fOO7gWQNSM74FqGc\nLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jO\ngBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQ\ndAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKS\nMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgB\nJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFIWCtP6WL3/ma0vd\nxGsDkjUgOWstSKNL5b1lgdvoH5CsAclZS0G6u1Tdp91G/4BkDUjOWgrSeTVIh652I7UByRqQ\nnLUUpMk1SIe84kZqA5I1IDlrKUi31yCd6Db6ByRrQHLWUpDW/7eKozc+6jb6ByRrQHLWUpBi\n9ZUnf3DUr9zEawOSNSA5ay1IBzYgWQOSMyDlbAHJGZACSLoFJGdACiDpFpCcASmApFtAcgak\nAJJuAckZkAJIugUkZ0AKIOkWkJwBKYCkW0ByBqQAkm4ByRmQAki6BSRnQAog6RaQnAEpgKRb\nQHIGpACSbgHJGZACSLoFJGdACiDpFpCcASmApFtAcgakAJJuAckZkAJIugUkZ0AKIOkWkJwB\nKYCkW0ByBqQAkm4ByRmQAki6BSRnQAog6RaQnAEpHEhbbpg45qp1tddf+e7YMy99Nsu+3lHe\naCA1HpC8VBtAmjFtxerrLuitvNp37pxtO/9hdFc26cGoPFeA1HBA8lLDH1KMWF5+qzRycRVM\nx3NZtrFjWXb6U/tcA6T9BiQvNfwhPT6qr/zywvuqv/ir2V077pm8a3fHzVPOmdkJpMYDkpca\n/pAeOrvy8vLbqr/YcEFHx4QXs03jb1y2bPr4reVbll1d3vM7zO3udh/RcLuyhK2eZKmdWcJW\nb7LUjixhqy9hKmUrXao3G/yapiFN2gupe8qcTdvuH1tDuH30vPLL+SeVt3DgBGPDdr17XhsE\n0oLau3b3V15fNKLi75y5tXvOv6f8YntnebHR3NYd7iMarivbmay1eVey1KYsXeu3PclSG7Pu\ndK3ehKmUrXSp7mzQSzY1C2nDiBeybPNpSyuvP92xrfxywtyVc7qzbMfo+a9dw8dI+42PkbzU\n8P8YKZs1dUXn9Ev6snlzs20T5mzZ9cCoNV1jZq/tnDlpJ5AaDkheqg0gbZs9YezM8uXXXpFl\nK68ae9Y3n8my5VecOW7Gq3suAdJ+A5KXagNITQxI+w1IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5\nKSABSQ5IXgpIQJIDkpcCEpDkgOSlgAQkOSB5KSABSQ5IXgpIlW3bam7nbvcRDbc9S9jqTpba\nliVs9SZLbc0StvoSplK20qV6s0Ev2ZYU0hZzO3a7j2i4MqR0re5kqa1ZwlZPstSWLGGrL2Eq\nZStdqicb9JKtKSHxrt1+4107L8W7dkCSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KX\nAhKQ5IDkpYAEJDkgeSkgAUkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgARpBVxcAAAfKSURB\nVEkOSF4KSECSA5KXAhKQ5IDkpYAEJDkgeSkgAUkOSF6qDSBtuWHimKvW1V5/5btjz7z02X1v\nA5IakLxUG0CaMW3F6usu6K282nfunG07/2F0V/1tQJIDkpca/pBixPLyW6CRi6tgOp7Lso0d\ny+pvA5IckLzU8If0+Ki+8ssL76v+4q9md+24Z/Kuuts2LCxvzSZz23a6j2i4LVnC1u5kqa4s\nYasnWWpTlrDVly7V25uulfBYPdmgl3Q1C+mhsysvL7+t+osNF3R0THix/rb5J5W3cOAEY8N2\nez6+GRTSpMrLGqTuKXM2bbt/7Ma625bfXN6L283t6nYf0XA7s4StnmSpHVnCVm+y1PYsYasv\nYSplK12qNxv8mmYhLai9G3d/5fVFI3aUX54zt/62yvgYab/xMZKXGv4fI20Y8UKWbT5taeX1\npzu2lV9OmFt/G5DkgOSlhj+kbNbUFZ3TL+nL5s3Ntk2Ys2XXA6PW7LkNSI0GJC/VBpC2zZ4w\ndmb58muvyLKVV40965vP7L0NSI0GJC/VBpCaGJD2G5C8FJCAJAckLwUkIMkByUsBCUhyQPJS\nQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJSQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJS\nQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJSQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJS\nQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJSQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJS\nQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJSQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJS\nQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJSQAKSHJC8FJCAJAckLwUkIMkByUsBCUhyQPJS\nQAKSHJC8FJCAJAckLwWkZiAtn3LyyVNerLsBSFYLSFYrXapgkF46vlTe8Sv23gIkqwUkq5Uu\nVTBIF5Sq+9reW4BktYBktdKlCgbpxBqkE/beAiSrBSSrlS5VMEgfq0H6L3tvAZLVApLVSpcq\nGKSLapAu3HsLkKwWkKxWulTBIK38QMXRB1buvQVIVgtIVitdqmCQ4qVvfPpT33ip7gYgWS0g\nWa10qaJB2m9AslpAslrpUkDK2QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQ\ndAtIzoAUQNItIDkDUgBJt4DkDEgBJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKS\nMyAFkHQLSM6AFEDSLSA5A1IASbeA5AxIASTdApIzIAWQdAtIzoAUQNItIDkDUgBJt4DkDEgB\nJN0CkjMgBZB0C0jOgBRA0i0gOQNSAEm3gOQMSAEk3QKSMyAFkHQLSM6AFEDSLSA5A1IASbeA\n5AxIASTdApIzIAWQdAtIzoAUiSEN5VZfPW+oj6C2+eoHhvoIclffOdQnkLv5pqE+gdxdV/c1\nf3GLQ3rmpBuG+ghq60765lAfQe6ks4f6BHL/44tDfQK5SScBaYgHJGtAGvIByRqQnAFpyAck\na0BijFUHJMYSDEiMJRiQGEuwFobU+ZenVf5nyw0Tx1y1bqgPs2cbrht3xqXLCnesVTPGfPVb\nzxXuWOU90vFE4Y719Y7yRlvHal1Ij02YXYU0Y9qK1ddd0DvUx3ltF09bvub6sTsKdqzuiTd2\nrpl91vaCHau8344f9UTh/iNOejAq3+flHKt1IT26/okKpBixvPz/HCMXD/Vx+tc1c1WWre94\nvmDH2vRP28tvwjuWF+xY5c36wfgnivYfMTv9qer/OMdqXUhZVoX0+KjKJ/svvG+oD1O/507b\nWMBjdd36td2FO9bjk3eUIRXsWLs7bp5yzsxO61gtD+mh6tcYL79tiM9Sv67z7yzesXr/vOOy\n3xTuWFsm/DorQyrYsTaNv3HZsunjtzrHan1IkyqvFua/QXmvnHdrXxGP9cys87YU7Vjf+15W\nhVSwY1W2ffQ851gtD2lB7c3v/UN9mD1bPObBrIDHKq/3zJ8V7Fi/ntBVhVSwY9V2/j3OsVoe\n0oYRL2TZ5tOWDvVhXtuzX/3Xyv8U7FiLzt2ZZX1jf1awY107asyYMSPOmFmwY62c051lO0bP\nd47VupA2xrzTInZks6au6Jx+ifHtha/rdp37o8rfnCzasbaMv2bV2ttGrS3Ysboqf1bj5m0u\n2rHGzF7bOXPSTudYrQvpf1a+aNbxz9m22RPGztw4+PUHZ4urp+r4WcGOla28cvQZ31icFe1Y\nlZXftSvasZZfcea4Ga9ax2pdSIwVaEBiLMGAxFiCAYmxBAMSYwkGJMYSDEiMJRiQGEswIDGW\nYEBiLMGAxFiCAamFN7d0ffnlj0pzhvogDEgtvYlHvZxtfs/ni/Jd0+08ILXyNh03Irvg6JVD\nfQwGpBbfvNLlh94+1IdgGZBafeeWvjDUR2CVAaml1/uZQ/5z11AfgmVAavFd/3sPvuXcoT4E\ny4DU2vv3I6dlt5QeGupjMCC19Ho+8f7tWe8nj/3tUB+EAamV9zelR8svl75hwlAfhAGJsRQD\nEmMJBiTGEgxIjCUYkBhLMCAxlmBAYizBgMRYggGJsQQDEmMJBiTGEgxIjCXY/wcP2Md9ngAe\ntgAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "check_N <- function(M,beta,sd){\n",
    "    x <- seq(4,50,2)\n",
    "    y <- sapply(x, function(N){\n",
    "        avg_simulations(M, N, beta, sd)\n",
    "    })\n",
    "    qplot(x,y)\n",
    "}\n",
    "check_N(400,0.1,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "\t<li>TRUE</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\item TRUE\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. TRUE\n",
       "2. TRUE\n",
       "3. TRUE\n",
       "4. TRUE\n",
       "5. TRUE\n",
       "6. TRUE\n",
       "7. TRUE\n",
       "8. TRUE\n",
       "9. TRUE\n",
       "10. TRUE\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0</li>\n",
       "\t<li>1</li>\n",
       "\t<li>1</li>\n",
       "\t<li>1</li>\n",
       "\t<li>1</li>\n",
       "\t<li>1</li>\n",
       "\t<li>1</li>\n",
       "\t<li>1</li>\n",
       "\t<li>1</li>\n",
       "\t<li>1</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0\n",
       "2. 1\n",
       "3. 1\n",
       "4. 1\n",
       "5. 1\n",
       "6. 1\n",
       "7. 1\n",
       "8. 1\n",
       "9. 1\n",
       "10. 1\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1] 0 1 1 1 1 1 1 1 1 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "0.953333333333333"
      ],
      "text/latex": [
       "0.953333333333333"
      ],
      "text/markdown": [
       "0.953333333333333"
      ],
      "text/plain": [
       "[1] 0.9533333"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# some tryouts\n",
    "sapply(1:M, function(x) simulate(100,.5,.5))\n",
    "       \n",
    "M <- 10\n",
    "output <- rep(0,M)\n",
    "for (i in 1:M){\n",
    "    output[i]<-simulate(100,.5,.5)\n",
    "} \n",
    "output\n",
    "\n",
    "       \n",
    "       \n",
    "run <- function(M, N, beta, sd){\n",
    "    res <- sapply(1:M, function(x) simulate(100,.5,.5))\n",
    "    sum(res)/M\n",
    "}\n",
    "    \n",
    "run(150,5,3,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#some useful stuff\n",
    "foo <- function(d){\n",
    "    nrow(d)\n",
    "    #d will be a vector in summarize\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## Writing Functions\n",
    "\n",
    "The goal of this exercise is to replicate the process from the slides, but with multivariate data. You should think of this as a practise of two things: \n",
    "\n",
    "1. Performing linear algebra in R. This is straight forward, but takes some practise. Especially worth noting how to use vectors and scalars together, as you see how to modify the functions we originally wrote for scalars, such that they work for vectors.\n",
    "\n",
    "2. Writing small functions and testing them. Make sure each part works before trying to combine them into a whole!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "── \u001b[1mAttaching packages\u001b[22m ─────────────────────────────────────── tidyverse 1.2.1 ──\n",
      "\u001b[32m✔\u001b[39m \u001b[34mggplot2\u001b[39m 3.2.1     \u001b[32m✔\u001b[39m \u001b[34mpurrr  \u001b[39m 0.3.2\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtibble \u001b[39m 2.1.3     \u001b[32m✔\u001b[39m \u001b[34mdplyr  \u001b[39m 0.8.3\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtidyr  \u001b[39m 0.8.3     \u001b[32m✔\u001b[39m \u001b[34mstringr\u001b[39m 1.4.0\n",
      "\u001b[32m✔\u001b[39m \u001b[34mreadr  \u001b[39m 1.3.1     \u001b[32m✔\u001b[39m \u001b[34mforcats\u001b[39m 0.4.0\n",
      "── \u001b[1mConflicts\u001b[22m ────────────────────────────────────────── tidyverse_conflicts() ──\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mfilter()\u001b[39m masks \u001b[34mstats\u001b[39m::filter()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mlag()\u001b[39m    masks \u001b[34mstats\u001b[39m::lag()\n"
     ]
    }
   ],
   "source": [
    "library(tidyverse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Write a function that generates data for regressions, \n",
    "\n",
    "# y should be generated according to: \n",
    "# y <- X %*% beta + eps\n",
    "# X should be several columns of independent random normals (aka a multivariate random \n",
    "# normal matrix with diagonal covariance)\n",
    "# eps should be such that:\n",
    "# eps <- rnorm(1, 0, sd)\n",
    "\n",
    "# beta should be a vector\n",
    "# params should be a DATAFRAME/TIBBLE with two columns: \"mean\" and \"sd\" to generate the features (X) \n",
    "# sd should be a scalar: the standard deviation of the normally distributed error term in the DGP\n",
    "\n",
    "# hint: use cbind to bind vectors into the columns of a matrix or the library MASS with the function mvrnorm\n",
    "\n",
    "\n",
    "\n",
    "generate_data <- function(N, beta,params, sd){\n",
    "    x<- cbind(rnorm(N,params$mean[1],params$sd[1]),rnorm(N,params$mean[2],params$sd[1]),rnorm(N,params$mean[3],params$sd[1]))\n",
    "    eps<- rnorm(1,0,sd)\n",
    "    y<- x %*% beta + eps\n",
    "    list(x=x,y=y)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<dl>\n",
       "\t<dt>$x</dt>\n",
       "\t\t<dd><table>\n",
       "<caption>A matrix: 3 × 3 of type dbl</caption>\n",
       "<tbody>\n",
       "\t<tr><td>-0.1101786</td><td>-0.006410601</td><td>-0.02091049</td></tr>\n",
       "\t<tr><td>-0.0074748</td><td> 0.155288347</td><td>-0.31735853</td></tr>\n",
       "\t<tr><td>-0.0153277</td><td>-0.133405247</td><td> 0.17756909</td></tr>\n",
       "</tbody>\n",
       "</table>\n",
       "</dd>\n",
       "\t<dt>$y</dt>\n",
       "\t\t<dd><table>\n",
       "<caption>A matrix: 3 × 1 of type dbl</caption>\n",
       "<tbody>\n",
       "\t<tr><td>-1.5132096</td></tr>\n",
       "\t<tr><td>-2.0084973</td></tr>\n",
       "\t<tr><td>-0.9105734</td></tr>\n",
       "</tbody>\n",
       "</table>\n",
       "</dd>\n",
       "</dl>\n"
      ],
      "text/latex": [
       "\\begin{description}\n",
       "\\item[\\$x] A matrix: 3 × 3 of type dbl\n",
       "\\begin{tabular}{lll}\n",
       "\t -0.1101786 & -0.006410601 & -0.02091049\\\\\n",
       "\t -0.0074748 &  0.155288347 & -0.31735853\\\\\n",
       "\t -0.0153277 & -0.133405247 &  0.17756909\\\\\n",
       "\\end{tabular}\n",
       "\n",
       "\\item[\\$y] A matrix: 3 × 1 of type dbl\n",
       "\\begin{tabular}{l}\n",
       "\t -1.5132096\\\\\n",
       "\t -2.0084973\\\\\n",
       "\t -0.9105734\\\\\n",
       "\\end{tabular}\n",
       "\n",
       "\\end{description}\n"
      ],
      "text/markdown": [
       "$x\n",
       ":   \n",
       "A matrix: 3 × 3 of type dbl\n",
       "\n",
       "| -0.1101786 | -0.006410601 | -0.02091049 |\n",
       "| -0.0074748 |  0.155288347 | -0.31735853 |\n",
       "| -0.0153277 | -0.133405247 |  0.17756909 |\n",
       "\n",
       "\n",
       "$y\n",
       ":   \n",
       "A matrix: 3 × 1 of type dbl\n",
       "\n",
       "| -1.5132096 |\n",
       "| -2.0084973 |\n",
       "| -0.9105734 |\n",
       "\n",
       "\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "$x\n",
       "           [,1]         [,2]        [,3]\n",
       "[1,] -0.1101786 -0.006410601 -0.02091049\n",
       "[2,] -0.0074748  0.155288347 -0.31735853\n",
       "[3,] -0.0153277 -0.133405247  0.17756909\n",
       "\n",
       "$y\n",
       "           [,1]\n",
       "[1,] -1.5132096\n",
       "[2,] -2.0084973\n",
       "[3,] -0.9105734\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>-0.110178617033584</li>\n",
       "\t<li>-0.00747479998080793</li>\n",
       "\t<li>-0.0153276991439281</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item -0.110178617033584\n",
       "\\item -0.00747479998080793\n",
       "\\item -0.0153276991439281\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. -0.110178617033584\n",
       "2. -0.00747479998080793\n",
       "3. -0.0153276991439281\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] -0.1101786 -0.0074748 -0.0153277"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "d<-data.frame(list(mean=c(0,0,0), sd=c(.2, .5, .3)))\n",
    "dat <- generate_data(3,c(2,3,4),d,1)\n",
    "#a<-(c(rnorm(3,0,1),rnorm(3,0,1),rnorm(3,0,1)))\n",
    "dat\n",
    "dat$x[,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_regression <- function(y,x){\n",
    "    coef <- calc_coef(y,x)\n",
    "    se <- calc_se(y,x,coef)\n",
    "    list(coef=coef,se=se)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc_coef <- function(y,x) cov(x,y)/(var(x))\n",
    "\n",
    "calc_se <- function(y,x,coef){\n",
    "    n <- length(y)\n",
    "    eps <- y-x*coef\n",
    "    e_sd <- mean(eps^2)\n",
    "    se <- sqrt(e_sd /(n*var(x)))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "ename": "ERROR",
     "evalue": "Error in cov(x, y)/(var(x)): non-conformable arrays\n",
     "output_type": "error",
     "traceback": [
      "Error in cov(x, y)/(var(x)): non-conformable arrays\nTraceback:\n",
      "1. calc_coef(dat$y, dat$x)"
     ]
    }
   ],
   "source": [
    "calc_coef(dat$y,dat$x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 326,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write other helper functions, as in the slides, and eventually a function \"avg_simulations\" that takes\n",
    "# the same parameters as generate_data, plus an \"M\" parameter that controls how many simulations are run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tests\n",
    "\n",
    "Make sure the following assertions pass:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a <- avg_simulations(1000, \n",
    "                20, \n",
    "                c(1,2,1), \n",
    "                data.frame(list(mean=c(0,0,0), sd=c(.2, .5, .3))), \n",
    "                .5)\n",
    "\n",
    "stopifnot(round(a, 1) == .9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a <- avg_simulations(1000, \n",
    "                500, \n",
    "                c(1,2,1), \n",
    "                data.frame(list(mean=c(0,0,0), sd=c(.2, .5, .3))), \n",
    "                .5)\n",
    "\n",
    "stopifnot(a > .92)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  },
  "name": "exercises.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}