Exercise16.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/choiitea/Choi_DSPN_S24/blob/main/Exercise16.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "614c6a4a"
      },
      "source": [
        "# Exercise 16: Model selection"
      ],
      "id": "614c6a4a"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7e861f62"
      },
      "source": [
        "This homework is designed to give you practice implementing model selection techniques, including best subset selection and forward/backward stepwise selection.\n",
        "\n",
        "You won't need to load in any data for this homework, we will be simulating our own."
      ],
      "id": "7e861f62"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b065256e"
      },
      "source": [
        "---\n",
        "## 1. Best subset selection (4 points)"
      ],
      "id": "b065256e"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bfe412c5"
      },
      "source": [
        "In this question, we will first generate simulated data, and then use it to perform best subset selection.\n",
        "\n",
        "a) Use `rnorm()` to generate a dataset including a predictor $X$ of length $n = 100$ and a noise vector $\\epsilon$ of length $n = 100$. Generate data for a response variable $Y$ of length $n = 100$ according to the model\n",
        "\n",
        "$Y = \\beta_0 + \\beta_1X + \\beta_2X^2 + \\beta_3X^3 + \\epsilon$\n",
        "\n",
        "where $\\beta_0$, $\\beta_1$, $\\beta_2$, and $\\beta_3$ are constants of your choice.\n"
      ],
      "id": "bfe412c5"
    },
    {
      "cell_type": "code",
      "execution_count": 79,
      "metadata": {
        "id": "82a01047",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        },
        "outputId": "86b8c3b7-eac8-48e7-9aad-b5371afc7df5"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<table class=\"dataframe\">\n",
              "<caption>A data.frame: 6 × 2</caption>\n",
              "<thead>\n",
              "\t<tr><th></th><th scope=col>Y</th><th scope=col>X</th></tr>\n",
              "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
              "</thead>\n",
              "<tbody>\n",
              "\t<tr><th scope=row>1</th><td>-2.156658</td><td>-1.1066080</td></tr>\n",
              "\t<tr><th scope=row>2</th><td>-1.407943</td><td>-0.1120660</td></tr>\n",
              "\t<tr><th scope=row>3</th><td>-4.816910</td><td>-1.7375028</td></tr>\n",
              "\t<tr><th scope=row>4</th><td>-1.137927</td><td> 0.2295856</td></tr>\n",
              "\t<tr><th scope=row>5</th><td> 5.590512</td><td> 1.2597499</td></tr>\n",
              "\t<tr><th scope=row>6</th><td>-1.213036</td><td> 0.4484742</td></tr>\n",
              "</tbody>\n",
              "</table>\n"
            ],
            "text/markdown": "\nA data.frame: 6 × 2\n\n| <!--/--> | Y &lt;dbl&gt; | X &lt;dbl&gt; |\n|---|---|---|\n| 1 | -2.156658 | -1.1066080 |\n| 2 | -1.407943 | -0.1120660 |\n| 3 | -4.816910 | -1.7375028 |\n| 4 | -1.137927 |  0.2295856 |\n| 5 |  5.590512 |  1.2597499 |\n| 6 | -1.213036 |  0.4484742 |\n\n",
            "text/latex": "A data.frame: 6 × 2\n\\begin{tabular}{r|ll}\n  & Y & X\\\\\n  & <dbl> & <dbl>\\\\\n\\hline\n\t1 & -2.156658 & -1.1066080\\\\\n\t2 & -1.407943 & -0.1120660\\\\\n\t3 & -4.816910 & -1.7375028\\\\\n\t4 & -1.137927 &  0.2295856\\\\\n\t5 &  5.590512 &  1.2597499\\\\\n\t6 & -1.213036 &  0.4484742\\\\\n\\end{tabular}\n",
            "text/plain": [
              "  Y         X         \n",
              "1 -2.156658 -1.1066080\n",
              "2 -1.407943 -0.1120660\n",
              "3 -4.816910 -1.7375028\n",
              "4 -1.137927  0.2295856\n",
              "5  5.590512  1.2597499\n",
              "6 -1.213036  0.4484742"
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "# random seed 1997\n",
        "set.seed(1997)\n",
        "# parameters for generating data\n",
        "n <- 100\n",
        "random_coef <- runif(4, min = -3, max = 3)\n",
        "#random_coef\n",
        "# -1.6716605569236 0.789223083760589 1.90374576766044 1.54804862150922\n",
        "# simplifying for my convenience\n",
        "beta0 <- -1.67\n",
        "beta1 <- 0.79\n",
        "beta2 <- 1.9\n",
        "beta3 <- 1.55\n",
        "\n",
        "# Generate data\n",
        "X <- rnorm(n)\n",
        "e <- rnorm(n, mean = 0, sd = 0.5)\n",
        "Y <- beta0 + beta1*X + beta2*X^2 + beta3*X^3 + e\n",
        "\n",
        "data <- data.frame(Y = Y, X = X)\n",
        "\n",
        "head(data)"
      ],
      "id": "82a01047"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "19065c60"
      },
      "source": [
        "b) Use `regsubsets()` to perform best subset selection to determine the best model that contains the predictor variables $X$, $X^2$, . . . , $X^{10}$. Print the model summary."
      ],
      "id": "19065c60"
    },
    {
      "cell_type": "code",
      "execution_count": 80,
      "metadata": {
        "id": "2d0ea51f",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 464
        },
        "outputId": "e19e89a9-98d3-45d4-bb64-5fc2c0f43083"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "Subset selection object\n",
              "Call: regsubsets.formula(Y ~ ., data)\n",
              "10 Variables  (and intercept)\n",
              "    Forced in Forced out\n",
              "X       FALSE      FALSE\n",
              "X2      FALSE      FALSE\n",
              "X3      FALSE      FALSE\n",
              "X4      FALSE      FALSE\n",
              "X5      FALSE      FALSE\n",
              "X6      FALSE      FALSE\n",
              "X7      FALSE      FALSE\n",
              "X8      FALSE      FALSE\n",
              "X9      FALSE      FALSE\n",
              "X10     FALSE      FALSE\n",
              "1 subsets of each size up to 8\n",
              "Selection Algorithm: exhaustive\n",
              "         X   X2  X3  X4  X5  X6  X7  X8  X9  X10\n",
              "1  ( 1 ) \" \" \" \" \"*\" \" \" \" \" \" \" \" \" \" \" \" \" \" \"\n",
              "2  ( 1 ) \" \" \"*\" \"*\" \" \" \" \" \" \" \" \" \" \" \" \" \" \"\n",
              "3  ( 1 ) \"*\" \"*\" \"*\" \" \" \" \" \" \" \" \" \" \" \" \" \" \"\n",
              "4  ( 1 ) \"*\" \"*\" \"*\" \" \" \"*\" \" \" \" \" \" \" \" \" \" \"\n",
              "5  ( 1 ) \"*\" \"*\" \"*\" \"*\" \" \" \"*\" \" \" \" \" \" \" \" \"\n",
              "6  ( 1 ) \"*\" \"*\" \"*\" \" \" \"*\" \" \" \"*\" \" \" \"*\" \" \"\n",
              "7  ( 1 ) \"*\" \"*\" \"*\" \"*\" \"*\" \" \" \"*\" \" \" \"*\" \" \"\n",
              "8  ( 1 ) \"*\" \"*\" \"*\" \"*\" \"*\" \"*\" \" \" \"*\" \" \" \"*\""
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "#install.packages(\"leaps\")\n",
        "library(leaps)\n",
        "#?regsubsets\n",
        "\n",
        "# I guess I need to create the predictor variables X^2 to X^10\n",
        "data$X2 <- data$X^2\n",
        "data$X3 <- data$X^3\n",
        "data$X4 <- data$X^4\n",
        "data$X5 <- data$X^5\n",
        "data$X6 <- data$X^6\n",
        "data$X7 <- data$X^7\n",
        "data$X8 <- data$X^8\n",
        "data$X9 <- data$X^9\n",
        "data$X10 <- data$X^10\n",
        "\n",
        "regfit.full = regsubsets(Y~., data) #use every variable\n",
        "summary(regfit.full)"
      ],
      "id": "2d0ea51f"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "e584b852",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# Note: if your model summary doesn't show up when viewing your notebook on github, copy and paste the output below.\n",
        "# Otherwise you can either delete this cell or leave it blank.\n",
        "\n",
        "```\n",
        "# paste here\n",
        "\n",
        "```"
      ],
      "id": "e584b852"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3528231b"
      },
      "source": [
        "c) Plot the Mallows' Cp, Bayesian Information Criteria (BIC) and the adjusted coefficient of determination, $R^2$, for each model tested. Which is the best model? Report the coefficients of the best model."
      ],
      "id": "3528231b"
    },
    {
      "cell_type": "code",
      "execution_count": 81,
      "metadata": {
        "id": "6a817913",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 437
        },
        "outputId": "f8e020fc-4e50-4e3d-e067-9389ffc12fe1"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "plot without title"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAC91BMVEUAAAAAv8QBAQECAgID\nAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQV\nFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYn\nJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5\nOTk6Ojo7Ozs8PDw9PT0+Pj5AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlLS0tNTU1O\nTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19g\nYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFy\ncnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OE\nhISFhYWGhoaHh4eIiIiKioqLi4uMjIyNjY2Ojo6Pj4+QkJCSkpKTk5OUlJSVlZWWlpaXl5eY\nmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamq\nqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8\nvLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3O\nzs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g\n4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy\n8vLz8/P09PT19fX29vb39/f4dm34+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///8T+uNrAAAACXBI\nWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3de2BU9Z338RMIqAiBVcELYG3rttpK8XGxjygg\ninXtAiL1BgrecGvxqRa1aNUHqU+hKkXduo/Wy3YttV23WgQLIvqsworFG16QKgIqlxBCuIRL\nwjW/P55zvknOTMLMkJnPLyQzvN9/zDmZmXxzkpwXczI5GQJHRHJBS28AUSEEJCIPAYnIQ0Ai\n8hCQiDwEJCIPAYnIQ0Ai8tD+IH1574jLf7bEuS1TRo+YWLbvsq53X0nfyxluy6G5c/zOe2Wu\n53GeP9/Wvn1zPW/fy5m2r3kxKO0H0u7Rv1q1Zupl292945evvn/s3n2Wdb2xNf2MTb62tbaq\ndX7n7a30O2/ber/zdmf40ubSlg1+5+3c7nfe5gz7y1y/H8pn+4G06bnwy7Rq8LLyIcvCb8GF\nixov6+8HpDggaRUmpKjKR27Y9cbwmnDtxj82XtbfB0hxQNIqVEh7Lxp8x3o3+6po/c7HGi/D\ni6o1YfPKt6ZtXfqbcmnjGr/ztmTY9FyqKPU7b3OF33nry/zO27TB77zyDPvLy81qQWr/j0gr\nP5h8/ZbZV0erIaBGy/Ci8uOw11ZXpG1N+ptyqWyV33nrS/3OW5vhS5FL5Wv9zvO9fevK/M4r\nzbC/zGlOClpNefp776Uz36w9lHu28bL+LhzaxXFop1WYh3bvjtnhXM3ImRVDloaf4tCPGi/r\n7wekOCBpFSakLVf+8svSx4aXusk3L191z7iafZZ1ASkOSFqFCcl9PuHiS25dFO4eU0eNnLRh\n32VdQIoDklaBQmpiQIoDkhaQ0gUkKSAlApLHgKQFpGYJSEDSApIFJCBpAckCEpC0gGQBCUha\nQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0g\nWUACkhaQLCABSQtIFpCApAUkC0hA0gKS1dyQqi0/H6Q2IGkBqVlqZkjV9fn5MFFA0gJSs9S8\nkKqBJAekREDyKAlIWkBqlpoVUjWQ9ICUCEhAyjkgJQISkHIOSIkOVkj8jOQhICUCEs/a5RyQ\nEh20kPg9kh6QEh28kDizQQ5IiQ5iSI5z7cSAlAhIHgOSFpCaJSABSQtIFpCApAUkC0hA0gKS\nBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckC\nEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJ\nSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQk\nLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIW\nkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtI\nFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQL\nSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUk\nIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQ\ntIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGR5gjSv\ndFPa1qa/KZfWr/Y7b2OZ33nla/zO27DO77x1Gb5VuVRR7nfe2gz7yxw/e2tzxCMSj0haPCJZ\nQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwg\nAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCA\npAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDS\nApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkB\nyQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBk\nAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKA\nBCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUAC\nkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJ\nC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQF\nJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKS\nBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckC\nEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJ\nSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJ2h+kivuvuOT2T8Iv/5TR\nIyaW7busC0hxQNIqUEg/Gb9szQMjq9y945evvn/s3n2WdQEpDkhahQmpctKXzq0b/Gn5kGXh\nt+DCRY2X9fcDUhyQtAoTkrVk6IY3hteEKzf+sfGy/i5AigOSVuFCqvzRU272VdHanY81XoYX\nFQvCXl1RmraV6W/KpdVf+p23ZpXfeata+/b5/n6s9jtvZYbte6kZJYjtH9LK6x+pcbOvjlZD\nQI2W4cWebWHzN+9O24b0N+XS1jK/83Zu8juvcp3feTsyfGlzadN6v/O2b/E7b2OG/SWfH5EW\njZgRXr5Zeyj3bONl/b04tIvj0E6rQA/tFl/+drSoGLI0/BSHftR4WX83IMUBSaswIe0c80x5\nWJWbfPPyVfeMq9lnWReQ4oCkVZiQFg22ZrptU0eNnBR+Cxov6wJSHJC0ChNSUwNSHJC0gJQu\nIEkBKRGQPAYkLSA1S0ACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCy\ngAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlA\nApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCAB\nSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICk\nBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNIC\nkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJ\nAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQB\nCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAE\nJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZyZCq\n/jo/150eSHFA0sp/SA+WBEH7H+3IaQyQ4oCklfeQ/hR8ddxd/YMbcxoDpDggaeU9pAFf3RJe\nXtN+cy5jgBQHJK28h9Tp7uhyYTA/lzFAigOSVt5DCh6NLtcEs3MZA6Q4IGnlP6TfRJelwaxc\nxgApDkhaQEoXkKSAlCg/IN2xIOzFYGq0yHYMkOKApJX/kJLLdgyQ4oCklfeQJiSX7RggxQFJ\nK2dIE4Juu2rXrg36pbzHGSdlvTVrx/fq1Onkmz61+VElZ/8p5R051w5IWq0HUtvi521le+f2\n/VLeY+qkbDdm/pHFQ37xi0vbd5oZQbrjN7959K6vBA+mumcypNLaXXRBRbYfDkhJAUkrd0jt\nzhpiK9Pan97Pz7as7XrUX6Plkh5d1oaQ7KmDyhNKqlLcNQnSjC5TbXlSt/ey/ohAigOSlnBo\nN6nd2mjl/PP79QsXz3y3Q0mfZ5x7IXggeit42A7tBvR//buHdb9v1/junQYtc+7UU6P3GHrU\nPjdYdwS/r115Zcq6ekhuXPDXFB89AenTjj1esZX/Pva4rL/WQIoDkpYAaWnbiMyqtk/27efc\nH4JhM2deEITHZKM7fuE2H3dOjUEa1POcd1YOC86buOq1zv+UDKnRDda3j9yTPL8W0l0pT/5J\nQPpx+yV1a+8V/3w/n+w+ASkOSFoCpKrzTgmXv+xQeUY/5yaduzMc1m5kuAMeP8SNLfm89smG\nQcEi5+YFZ4V3HNmxAaSGN0TVFJ/TYH4tpP7tUm1hAtJJl8WrQ0/OtMGpAlIckLQUSL8LFjr3\nrctdBKm2ngPCiznBnW2j0w0MUoTks+DW8PLWoDIZUsMborYGQxvMf7G0dM3Ca4IbUn30BKSO\nv4hXJx6aaYNTBaQ4IGkpkLaV3OAWBrMM0ua7e3UuLq59InxMcF60MEgnhCsrgsnh5fhgYzKk\nhjdE1bQb0GC+1S71n+wlQZqceJeSTBucKiDFAUlLgeSu+buqscftMUhnF//s9Q8+7B6uub39\n25wQfZezheR6d9qePH/qrFmz521M/dETkHpdFK8O7J1pg1MFpDggaUmQXgumd73FRZCWBmPC\n63Yf1i+8fKDdjM7RWykgndYretcz0kC6J6h9Itu9cdKC+Gek1CUg/bT47bq1GcFdmTY4VUCK\nA5KWBKnm66cH7xukj4OJ4XUPBX2d+1uH8e5foj8PSgHp3K41zpV1SANp/bEdZkbLRd27VTQZ\n0tojuz4TPdm3fWqHbln/ShZIcUDSkiCFjyHRsVoIadfxPabPv2XgwJJXN/c9cbvbe2aPjakg\nPRhMXvvuOaekgeQWdgsG3XPv8OJj3nVNhuTmdQ2OOnfogE7BsQub9BknB6Q4IGlpkJa3meIM\nknvrzMOP+efNM7oecW0Q/X70o/ajUkHaOa7HoafOGFuSBpJbf3uvjiWn3l3usoDk1t7x7eKg\nXe97OEVICUhaeX/2t1WztSanMUCKA5JWYUDKNSDFAUkLSOkCkhSQEgHJY0DSAlKzBCQgaQHJ\nAhKQtIBkNYS0bY1z2596YFmaO6cPSHFA0ioESEuOnux2nx4EXd7NdgyQ4oCkVQiQLvrOZ+7p\n4JHPzvpBtmOAFAckrUKAdPQ054b1cm7a8dmOAVIckLQKAdIhr7o9R/zUuTmHZDvmvzfsSNv6\n9DflUuVav/O2Z9j0XNpU5nfeto1+520s9ztvy2a/8yoy7C8v+9vxfdcA0vGPuznBq849cVy2\nY4AUByStQoB07bG3n3DiHlfWm5+Rco9DO61COLRb0zfousC5S7u8n+0YIMUBSasQIIWfRfTq\nyW+tzXoMkOKApFUYkHINSHFA0sp7SOMmPPDongz3zBiQ4oCklfeQgradjk35il1NCUhxQNLK\ne0jbnFua8xggxQFJywuk6upqH9uSRQ1+RgpO/NH03HYrIMUBScsDpOraUt7Wp/b1Up9yruyQ\nnvajTJ+xtbfsmdy70yHfnLQ3vs+J2Wx2A0gPX9Q1aD9w8rvZv24DkOKApKVDqq7OIKnPFUuj\nwm/6LwYd/We7pg7ST7rPLF0zrfPd8X2+yGazG7/4yQcPhZiOviKbEVFAigOSlkdIqSTVs3F7\nT3jyx99PvqbX7dHlS7OS7pNFKZ7+Xv2/u/GfMecekLRkSNVNgzSj45b32q5IuuaKv3+78X2y\nqBGZ1b+//qSg5Pz/k+0YIMUBScsnpBSSYiTfv9q5036WdE3FZcVfu+LRsuia9h2jfp3NZjeA\ndN03gq4X/uqtHH6bBKQ4IGk1N6Q2xVFvL28737mHj9mV/PhT8dyt3zr03+OfkbLacRs+a9f1\n9o+zeedEQIoDklZzQ7r8w6iq8UGXLl06Bf/R+EDu5o679UO7f7/6q8HRF//6o+zHACkOSFoH\n5meknd0mrAj7waD4mi8utmfpnm271cuTDSuevLJn0I0/o8g9IGkdmGftph0SvS6+e63Np3UH\ncqV7evd+YcUX00+8ID60W7ori81O8fzc8ju78qxd7gFJy98vZFP/HqkWUv+6/zH55Fvqfv06\n3FWM+8bhh3xzfGX8C9lgSaoBaWpIpmbxv47oGXQ4/4EsJlhAigOSlo9ThDKc2dBcNYA0rGvQ\n5n/c9nIOp64CKQ5IWnl/0mpYj9HTynIbA6Q4IGkVAiTn1s987PHZOexZQIoDklYhQNp7S/vo\nZ6yO92U9BkhxQNIqBEj3BcOe+MvMR/8x+G22Y4AUByStQoD0rXG1y+v/IdsxQIoDklYhQDr0\nldrlix2yHQOkOCBpFQKkjjNql3/ulO0YIMUBSasQIPU/Z2e0qDp/YLZjgBQHJK1CgPRim6/8\n8N6fj+neNuvXWAZSHJC0CgGSe/7k6Onv77yY9RggxQFJqyAgObd6YQ4vWAykpICkVSCQcgxI\ncUDSyntIJyWX7RggxQFJK+8h9Usu2zFAigOSlh9IRR62JKs4tAOSViuFdKAlAQlIWq0SUlEL\nQuJnJD8BScsXpAMsiZ+RgKTVGiEVpYU0vO7lGEbXvjBD1398M7zyhAl2W4/Jub1+vpXy0G7L\np9mOAVIckLS8QUopqXTp0ueDl5YuLXN9rlq5cuXCSzotbwgph9fPt1JCmntktmOAFAckLQ+Q\niorSS3LurdpXB6p9OaHdnf6lIaQcXtLOaghp5sgB4XFd35Ku2Y4BUhyQtHRIRVlA2tvlV80A\n6ZmgXc+g+2HBOVmfbAekOCBpeYSUWlIypMpbD2l0aJfD6+dbDSD1uaDSFX+4+6GBWe9bQIoD\nkpYMqaiJkCIyQe+XXKqfkbLfZxtAKpnpXPEHzt2c9eMbkOKApOUTUkpJ9ZBCMu8cYQ8937g5\nutx7+IO+Du0Om+Vc59edm9c92zFAigOSlgqpqGg/kpIP7Z7usDi8HNon+s9e/18wzxek036w\n051yp3PTO2Y7BkhxQNISIRVlBcld1HuHc4s7XbFg8VPHjXS5vX6+1QDS08Egd3fxmIk9zsp2\nDJDigKTlF1IKSQ0grTv6pvDy/YuOO/yUKdFLdefy+vlWw6e/n5nstn0vCI5/K9sxQIoDklbe\n/xlF3NKPs35cA1IiIGkVDqRcAlIckLTyHhInrfoJSFp5Dylof3q/vmfUle0YIMUBSSvvIV3b\nuf35T+a6UwEpDkhaeQ/JbZ92ftsOl8/ak8sYIMUBSSv/IYWt+uW3guNueT/7MUCKA5JWQUAK\nWzi2a3DqlGzHACkOSFqFAslVTe8TtMl2DJDigKRVIJAW3vB3wdfvWZ7tGCDFAUmrECCtfeCU\noOSa12qyHwOkOCBp5T2kXc8Padf2vN9ty2kMkOKApJX3kLoVnzXps6q6sh0DpDggaeU9pKBt\nuyAu2zFAigOSVt5DmpBctmOAFAckrbyHJAWkOCBpASldQJICUiIgeQxIWkBqloAEJC0gWUAC\nkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJ\nC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQF\nJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKS\nBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckC\nEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJ\nSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUka7+QVt0yNFpsmTJ6xMSyfZd1ASkOSFoFCun1\nUVMN0r3jl6++f+zefZZ1ASkOSFoFCumVdQsiSOVDloXfggsXNV7W3w1IcUDSKlBIzhmkN4bX\nhJc3/rHxsv5OQIoDklZhQ5p9VbR652ONl+HFnm1h8zfvTtuG9Dfl0tYyv/N2bvI7r3Kd33k7\nMnxpc2nTer/ztm/xO29jhv2lACBdHa2GgBotw4uKBWGvrihN28r0N+XS6i/9zluzyu+8Va19\n+3x/P1b7nbcyw/a91IwSxJoI6c3aQ7lnGy/r78ShXRyHdlqFfWhXMWRp+CkO/ajxsv5OQIoD\nklaBQtpQPmdoeXmVm3zz8lX3jKvZZ1kXkOKApFWgkK4dHDXdbZs6auSk8FvQeFkXkOKApFWg\nkJoYkOKApAWkdAFJCkiJgOQxIGkBqVkCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQ\ngKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA\n0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBp\nAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSA\nZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCy\ngAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlA\nApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCAB\nSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICk\nBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNIC\nkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJ\nAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQB\nCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAE\nJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZniDN\nK6tMW4abcqlitd95m9f5nbe+1O+8jeV+55Wv9TtvQ4Xfeesy7C8v+9lbmyMekXhE0uIRyQIS\nkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlI\nWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQt\nIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQ\nLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gW\nkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtI\nQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQg\naQHJAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0\ngGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQFJAtIQNICkgUkIGkByQISkLSAZAEJSFpA\nsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKSBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZ\nQAKSFpAsIAFJC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkKz8\ng1RdXe11HpC0gGTlG6Tq2jxOBJIWkKw8g1Rd7V0SkLSAZOUXpGogyQGpWcpXSP4kAUkLSFZe\nQaoGkh6QmiUgAUkLSBaQgKQFJCuvIPEzkoeA1CzlKyR/M4GkBSQrvyA5/46AJAYkK88g1VHy\nORBIWkCy8g4SJ62KAalZAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJC0gWkICkBSQLSEDS\nApIFpL1FfucBSQtI6WrlkIqAJAUkC0hFniW1dkiez7DyfMZW5lMpgeQxz5CKPEPyfC6g53N0\nfZ/0e4DnAcljfiEVFXmV1Mp3VN9/hnKg5xUgpC1TRo+YWBa/CSQr33fUFp6X9HebqQcWIKR7\nxy9fff/YvfVv5iukoiKvkva3I7TwvP3uqK19XuFBKh+yLHxUunBR/dv5DsmPpLzfUVv7vMKD\n9MbwmvDyxj/Gb+cnpCIg5dW8woM0+6ro8s7HwovKj8NeW12RtjXpb8qlslX+ZiUgFfkYl7wj\nMK8Z5s3xuOd7LldIV0eXBqlqTdi88q1pW5f+plzauMbbqCRHRT7mJe8IPuZtzZ95Xsbtd97L\nHvd8z+UI6c3aQ7tn69/Oy0O7ogb5mOj3SOfgm3fw/R6pYshS5zYP/aj+7QKA5EPSgd6xCm7e\nwXdmw+Sbl4k9iJkAAAfMSURBVK+6Z1xN/Zv5CKnIP6QDvWO1/DzfpwgddOfabZs6auSkxOmP\n+QipLs9/RuF9x+Kk1UQFCKlRQIpr7SetAqlZAhKQtIBkAQlIWkCygAQkLSBZQAKSFpAsIAFJ\nC0gWkICkBSQLSEDSApIFJCBpAckCEpC0gGQBCUhaQLKABCQtIFlAApIWkCwgAUkLSBaQgKQF\nJAtIQNICkgUkIGkByQISkLSAZAEJSFpAsoAEJC0gWUACkhaQLCABSQtIFpCApAUkC0hA0gKS\nBSQgaQHJAhKQtIBkAQlIWkCygAQkLSBZviB9/FnaPkp/Uy797X2/85Yu9jtvyQd+532a4Uub\nSx9/6HfeJ0v8zlucYX+Z5WdvbY48QVr7RfqWZbgthz552++8z5f7nbfkHb/zPl/hd95Hi/zO\n8719H7yf4UY/e2tz5AnSAazi3Zbegsyt/aCltyBzK5e09BZkbvlnLb0FOQUk3wFJC0gHKCBp\nAalZApLvgKQFJKKDNyAReQhIRB4CEpGH8g1Sxf1XXHL7Jy29FRmbO3hBS29Chl68btiNC1t6\nI9K38ucjL719cUtvRQ7lG6SfjF+25oGRVS29GRnaeOXwVgxp7qi3yv48ZltLb0a6asY8vG3H\n7y72fPrjgSjPIFVO+tK5dYM/bentyNDkJ65sxZDGvNLSW5CxTYOXOLdhcOs+5EhZnkGylgz1\nfAKzz964rqoVQ1o/+JX/9YNbWvFvkm6bWln1++t2tvRmZF8eQqr80VMtvQnp2zLqPdeKIX0y\n+I6VlY9d5vkvWzxWMXbw4FH5+CvZ/IO08vpHalp6G9L34IOudUNa5Nyey1vt3/XsvunhTdue\nHdmKjzjSlXeQFo2Y0dKbkKH3RlW2akjlg5eGl2OfbentSNe7Q6Lnka55oaW3I/vyDdLiy99u\n6U3I1H3DR4wYMeSSSS29HenaOyr8Z2jnJa+39Hak653B0ROKo4DU3O0c80x5WKt9+rsy2ror\n5mxu6e1I27Mj3yt/aFSr/fptG/Xwlp1/Gr6mpbcj+/IM0qLB1syW3o6MteJDO7f3t1cOu/3L\nlt6K9H0+ceRlP23l58+nLM8gEbXOgETkISAReQhIRB4CEpGHgETkISAReQhIRB4CUpObEPSt\nPVu2z6AM9+p3UhPH7b7y8A4ro5Xvta//bwFO6tb47wfOaDAtMbvJH4UOUEBqchOC4FFb8QNp\nZjDyBftL1WeDKbXXvBHc2vhOUxuctQek1huQmtyEw75/hD12+IH0b0HdqaO7ju5VuzIm+Fvm\ndwFS6w1ITW5C8Olho6OVCNKpp0ZrQ49ybkD/1797WPf7do3v3mnQsnAXP/md/ocfMWpjeOt/\nnVfS4bQnwpV+A2b0PLNuyl8GdDrslCk1blAQtsKuui2wVyPZ3vls5575boeSPs/E72OHdonr\n4tkGKZ6/5rqvHHrMRa34714PgoDU5CYEOyYGr7nGkAb1POedlcOC8yaueq3zP4W7eM+T7nv+\ntjaDnZtbfPaMOT8MHnDu3N4n/7ruPNvn21zw57njgtvcJxOCx9+q/ZHo0+CH0eLp4Gn3h2DY\nzJkXBDPr3yeClLguMTuClJjf99jHX532naNb7UuaHAwBqclNCKp2fOPbu/aBFCxybl5wVvjm\nyI7hLh78Z7g2IvjCnfb30Z49pKQqvMtz9UNO/kqE58L2691Twbz6Kwd2if7Xu3OPqHKTzg1v\n3txupKt7nwhS4rrE7AhSPH9zcHu48tmk1Qfwi0GNAlKTCyG5l4Jf7gMp1OM+s+cJbg0qXb9D\nQ2qhkufKgpuqwv5veNw26JBddTNW1z74PBE+viRBmhZMc+7zNj+uf7vnAFf3Poln7aLr4tkR\npMT8XUedMHdv83/+lCkgNbkIkrv48M8bQzohXFkRTA4vxwcbXb+vRTfMCh59L6jrOTeoe/2M\nhcG90eIvwWPJkHYcGU6cGHwYPvDc3atzcXHQz9W9TwQpcV08O4KUNH/+14Kjhk/bfaC+EJQi\nIDU5g7SqZIg7PSOkE6Mb/hL85r3gmgVWee1drLeCidHixeDxZEjupjaf13y9b7hydvHPXv/g\nw+796sYapMR18exaSPF8t+eVW78dnO75f0WmbAJSkzNIbkow/cwQ0mn2lPUZKSB1jI6yngxe\nqAhG179nAlJpcH20eCyY3QDS4mDSvOBJ55YGY8K3dh/WLwlS0nXx7AhS0nzrkeDfmuOzpqYF\npCZXC2n3d044K4R0btca58o6pIAUzAnXLmxb6v5nl+g58N/euTsJkuvVPRpyweGbG0ByZ/W5\nqfNW5z62B6yHgr5JkJKuS8yOnmyI5799aZmLfky7/8B8HShVQGpytZDc/DZtQkgPBpPXvnvO\nKftCOrPnN/917vjgcuf+q33v3750V/urkh+R3Ittz58+64bo3g0gPRUcET0Lsev4HtPn3zJw\nYMmrW2NISdclZkeQ4vmlJb2fePkPZ3XOx9dVLJiA1OTqILmrgxDSznE9Dj11xtiSfSD9w5lv\n9+9wxHVbwjfnfa+k/Tfv290AkpvTv+Ohp4VHcQ0hbesSvBMt3zrz8GP+efOMrkd8kvgZKXFd\nYrb9Qjae//6wo9t3H9bK/0vQAg9IRB4CEpGHgETkISAReQhIRB4CEpGHgETkISAReQhIRB4C\nEpGHgETkISAReej/A/i2V8aLMq+mAAAAAElFTkSuQmCC"
          },
          "metadata": {
            "image/png": {
              "width": 420,
              "height": 420
            }
          }
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "library(tidyverse)\n",
        "reg.summary = summary(regfit.full)\n",
        "numvar = 1:length(reg.summary$rss)# Make a vector that lists the number of variables in sequence, from 1 to 19.\n",
        "allfalse = rep(FALSE,length(reg.summary$rss))# Starting point for an indicator that marks the best model choice for each metric.\n",
        "\n",
        "# mallow's cp\n",
        "cp.df <- data.frame(numvar = numvar, cp <- reg.summary$cp, mincp <- allfalse)\n",
        "cp.df$mincp[which.min(reg.summary$cp)] <- TRUE\n",
        "ggplot(cp.df,aes(x=numvar,y=cp,shape=mincp,col=mincp)) +\n",
        "    geom_point(size=3) + theme_light() +\n",
        "    labs(x = \"Number of Variables\", y = \"Mallow's CP\", color=\"Maximum CP\", shape=\"Maximum CP\")"
      ],
      "id": "6a817913"
    },
    {
      "cell_type": "code",
      "source": [
        "# #bic\n",
        "bic.df <- data.frame(numvar = numvar,bic <- reg.summary$bic, minbic <- allfalse)\n",
        "bic.df$minbic[which.min(reg.summary$bic)] <- TRUE\n",
        "ggplot(bic.df,aes(x=numvar,y=bic,shape=minbic,col=minbic)) +\n",
        "    geom_point(size=3) + theme_light() +\n",
        "    labs(x = \"Number of Variables\", y = \"BIC\", color=\"Minimum BIC\", shape=\"Minimum BIC\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 437
        },
        "id": "xeSPT4ghmwP4",
        "outputId": "64e0191d-046a-4d02-8920-f372733fa8cd"
      },
      "id": "xeSPT4ghmwP4",
      "execution_count": 82,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "plot without title"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAC7lBMVEUAAAAAv8QBAQECAgID\nAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQV\nFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8hISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj5AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlLS0tNTU1OTk5P\nT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBh\nYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJz\nc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1/f3+AgICBgYGCgoKDg4OFhYWGhoaH\nh4eIiIiKioqLi4uMjIyNjY2Ojo6Pj4+QkJCSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqb\nm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKyt\nra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/\nv7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR\n0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj\n4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT1\n9fX29vb39/f4dm34+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///9SOdBcAAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAgAElEQVR4nO3de4DVdYH38d8g3lDQx7wVumvbRbdcfHysrbyUYj21MBppCoLa\nas+WlWa2aBckgzTDetaeNLO27LLuZoUVeWnXUMwrlKmrCAwXgWG4OQgMF5nff8/vnDln5jvw\nPcPMfD5f5sL7/cf8zpzf8OEXzatzmDlMWU5EcllfXwDRYAhIRIaARGQISESGgERkCEhEhoBE\nZAhIRIYckJ5+sHb3d3GuN91vHzTvPdDvL/A+96B578GuLvAZw+drkhyQ5rxa+9x6w37YK+vM\ng+4L3NRkHnRf4Pbl5sFXt5sHl9Q+tfZJ8+9lC0jmPSDJASkSkNSAFAQkV0BSA1KSgGTeA5Ic\nkCIBSQ1IQUByBSQ1ICUJSOY9IMkBKRKQ1IAUBCRXQFIDUpKAZN4DkhyQIgFJDUhBQHIFJDUg\nJQlI5j0gyQEpEpDUgBQEJFdAUgNSkoBk3gOSHJAiAUkNSEFAcgUkNSAlCUjmPSDJASkSkNSA\nFAQkV0BSA1KSgGTeA5IckCIBSQ1IQUByBSQ1ICXJAemRtVtqtrr2qV61psk86L7AV1aaB90X\nuHGZeXDdRvPg4tqnGp8wfL4myQHp4dWbatZU+1SvWr3KPOi+wHUrzYPuC9ywzDy4doN5cHHt\nUyseN3y+JomnduY9ntrJ7bVP7YAUBCQ5IO1cS1uG36I9IKkBKUkpIbW0+CUBSQ1ISUoIqQVI\njoAUtJdDMkoCkhqQkgQk8x6Q5IDUOSBZAlIQkFwBSQ1ISQKSeQ9IckDqHJAsASlor4TE95Es\nASlo74TEKxscASloL4WU81o7PSAFAckVkNSAlCQgmfeAJAekSEBSA1IQkFwBSQ1ISQKSeQ9I\nckCKBCQ1IAUByRWQ1ICUJCCZ94AkB6RIQFIDUhCQXAFJDUhJApJ5D0hyQIoEJDUgBQHJFZDU\ngJQkIJn3gCQHpEhAUgNSEJBcAUkNSEkCknkPSHJAigQkNSAFAckVkNSAlCQgmfeAJAekSEBS\nA1IQkFwBSQ1ISQKSeQ9IckCKBCQ1IAUByRWQ1ICUJCCZ94AkB6RIQFIDUhCQXAFJDUhJApJ5\nD0hyQIoEJDUgBQHJFZDUgJQkIJn3gCQHpEhAUgNSEJBcAUkNSEkCknkPSHJAigQkNSAFAckV\nkNSAlCQgmfeAJAekSEBSA1IQkFwBSQ1ISQKSeQ9IckCKBCQ1IAUByRWQ1ICUJCCZ94AkB6RI\nQFIDUhCQXAFJDUhJApJ5D0hyQIoEJDUgBQHJFZDUgJQkIJn3gCQHpEhAUgNSEJBcAUkNSEkC\nknkPSHJAigQkNSAFAckVkNSAlCQgmfeAJAekSEBSA1IQkFwBSQ1ISQKSeQ9IckCKBCQ1IAUB\nyRWQ1ICUJCCZ94AkB6RIQFIDUhCQXAFJDUhJApJ5D0hyQIoEJDUgBQHJFZDUgJQkIJn3gCQH\npEhAUgNSEJBcAUkNSEkCknkPSHJAigQkNSAFAckVkNSAlCQgmfeAJAekSEBSA1IQkFwBSQ1I\nSQKSeQ9IckCKBCQ1IAUByRWQ1ICUJCCZ94AkB6RIQFIDUhCQXAFJDUhJApJ5D0hyQIoEJDUg\nBQHJFZDUgJQkIJn3gCQ3eCGtuemi8699Ic833Hzx+KmNHcdKQAoCktzghfTZyQuWf33C5vyG\nyQtfvumKHe3HSkAKApLcoIXUPK34j7ZqzItNYxcUj0bnzKseq+eBFAQkuUELqdzz9WvnjGst\nbnzqZ9Vj9RSQgoAkN6ghNX/yznzWJaVbX7yteizevPpS0UMr1tdsZe1TvWplF79X7wbNe00v\nmwfdF7h2iXlw1VrzYEPtU0seSwVBbTeQZtfX1z9XHJd+/NbWfNalpbsKSJVj8Wbj4qKHGptr\n1sWpXtW40j1o3luz3DzovsD1S82DTevNg4trn1o2UCFtbGho2JLn88bPLN75Y9tTururx+oH\n8dQuiKd2coP3qd2zF5avf83Y+cVncv1fqsfqaSAFAUlu0ELaevlPmoo259OvWrjs+qtb24+V\ngBQEJLlBC2nemHL35htnTJowbW3efqwEpCAgyQ1aSLsLSEFAkgNSJCCpASkISK6ApAakJAHJ\nvAckOSBFApIakIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUCkhqQgoDk\nCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCCgOQKSGpAShKQzHtAkgNS\nJCCpASkISK6ApAakJAHJvAckOSBFApIakIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQk\nAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCC\ngOQKSGpAShKQzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIakIKA5ApIakBKEpDMe0CS\nA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQG\npCQBybwHJDkgRQKSGpCCgOQKSGpAShKQzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIa\nkIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7\nQJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCCgOQKSGpAShKQzHtAkgNSJCCpASkISK6A\npAakJAHJvAckOSBFApIakIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUC\nkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCCgOQKSGpAShKQ\nzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIakIKA5ApIakBKEpDMe0CSA1IkIKkBKQhI\nroCkBqQkAcm8ByS5vRbS7MbmmnVxqlc1rnQPmvfWLDcPui9w/VLzYNN68+Di2qeWPWb4fE2S\nA9IjzTtqtq72qV61fq150H2BG1eZB90XuHW5eXDDVvPgktqnVg/qRySe2gXx1E5ur31qB6Qg\nIMkBKRKQ1IAUBCRXQFIDUpKAZN4DkhyQIgFJDUhBQHIFJDUgJQlI5j0gyQEpEpDUgBQEJFdA\nUgNSkoBk3gOSHJAiAUkNSEFAcgUkNSAlCUjmPSDJASkSkNSAFAQkV0BSA1KSgGTeA5IckCIB\nSQ1IQUByBSQ1ICUJSOY9IMkBKRKQ1IAUBCRXQFIDUpKAZN4DkhyQIgFJDUhBQHIFJDUgJQlI\n5j0gyQEpEpDUgBQEJFdAUgNSkoBk3gOSHJAiAUkNSEFAcgUkNSAlCUjmPSDJASkSkNSAFAQk\nV0BSA1KSgGTeA5IckCIBSQ1IQUByBSQ1ICUJSOY9IMkBKRKQ1IAUBCRXQFIDUpKAZN4DkhyQ\nIgFJDUhBQHIFJDUgJQlI5j0gyQEpEpDUgBQEJFdAUgNSkoBk3gOSHJAiAUkNSEFAcgUkNSAl\nCUjmPSDJASkSkNSAFAQkV0BSA1KSgGTeA5IckCIBSQ1IQUByBSQ1ICUJSOY9IMkBKRKQ1IAU\nBCRXQFIDUpKAZN4DkhyQIgFJDUhBQHIFJDUgJQlI5j0gyQEpEpDUgBQEJFdAUgNSkoBk3gOS\nHJAiAUkNSEFAcgUkNSAlCUjmPSDJASkSkNSAFAQkV0BSA1KSgGTeA5IckCIBSQ1IQUByBSQ1\nICUJSOY9IMkBKRKQ1IAUBCRXQFIDUpKAZN4DkhyQIgFJDUhBQHIFJDUgJQlI5j0gyQEpEpDU\ngBQEJFdAUgNSkoBk3gOSHJAiAUkNSEFAcgUkNSAlCUjmPSDJASkSkNSAFAQkV0BSA1KSgGTe\nA5IckCIBSQ1IQUByBSQ1ICUJSOY9IMkBKRKQ1IAUBCRXQFIDUpKAZN4DkhyQIgFJDUhBQHIF\nJDUgJQlI5j0gyQEpEpDUgBQEJFdAUgNSkroFackN4y/8wvN5vuHmi8dPbew4VgJSEJDkBi2k\n7Rd/Y9nyGRdsym+YvPDlm67Y0X6sBKQgIMkNWkjrf74pz5eNWdA0dkHxaHTOvOqxeh5IQUCS\nG7SQSjXf+oltc8a1Frc+9bPqsXoOSEFAkhvEkHZ8eMx1q/NZl5Ruf/G26rF4s+7Jov9c3Fiz\nZbVP9aplS92D5r3l/f0CV3bx31avenmlebCh9qkFj6ajoLUbSLPr6+ufK45L/zz94xtmXVq6\nq4BUORZvtq0venjtlpqtrn2qV61pMg+6L/CVleZB9wVuXGYeXLfRPLi49qnGJ1KD6G27gbSx\noaFhS/nWjo/e+8e2p3R3V4/VD+KpXRBP7eQG7VO7py8vLLVOuHfN2PnFZ3L9X6rH6nkgBQFJ\nbtBC2jDxa0tW3DZuRT79qoXLrr+6tf1YCUhBQJIbtJDyhinnnX/NvOKJ3oxJE6at7ThWAlIQ\nkOQGL6TdBKQgIMn1GtKU7Mhtbbf+MTstz991fHiy83vdaEpWNPz4y0qXc0/27dJdKyefOHz4\nCVe+GPloIJn3gCTXe0j7DL2nfGPTIfsVkGZMC092fq8bTcmuu/32Wybu+4Y1VUgPv27o2K9+\n9aP7Db93148GknkPSHK9h7TvqWPLN+7a7x2nyZcxJSt/02p6dmsF0sojDn+sdM/zxxy6cpeP\nBpJ5D0hywlO7afuWP8U/8IHTKk/tzjj96bNGHHlBY/W9P7zzwJE3bps8cvjoBXl+0kmlD64/\nfJcTlbUypIezqRVI12U/bjvz4M2rdvm9gWTeA5KcAGn+Pl8vjsv2+d67K5BG/9U772/896EX\nV9479synlp6bnT112UOH/EMIaacTlbUypDuyeyqQ3va612r/3kAy7wFJToC0+ey3F8evDWt+\nVxVS9nBxx+iR1ffm5fns7NTirgkHd4LU+URlbVZT06KfHn3WjjZIrUPP7OL3BpJ5D0hyCqQf\nZY/n+d9emLdDOqh04uJ9Ku+VkLyUXVO8vSZrDiF1PlFZK3fWorwN0qtZfRe/N5DMe0CSUyBt\nHPGJ/PHstx2Qjiud+McseG9RNr14OzlbF0LqfKKy9s377//td88e9sPKI9K+Z3TxewPJvAck\nOQVS/rH/sfmKN7zmgVT+O1Lr2Qetb/s70qjhm2r/3kAy7wFJToL0UPbLIz6Xdw/SySeWTr6r\na0j5ddmcNkjXZzPazsw5ftd/zQEk8x6Q5CRIrX/zjuxP3YR01hGted44rGtI2/9+yJI2SKtf\nP6z8ndh5I49cs8vvDSTzHpDkJEjF40bpCVu3IN2STV/59Jlvrwnputtv/860U7JPVl/Z8PiR\n2ejrbxg39Oind/29gWTeA5KcBmnhkJvzbkLaevUxB5w084oRtSAVDTn6zB+3tr/WbvW1Jx48\n4qQvx/4rBpJ5D0hyvPo7EpDUgBQEJFdAUgNSkoBk3gOSHJAiAUkNSEFAcgUkNSAlCUjmPSDJ\nASkSkNSAFAQkV0BSA1KSgGTeA5IckCIBSQ1IQUByBSQ1ICUphLSi7WejPLrra8S7DkhBQJIb\n6JBmHtr275aOP3JuzzaAFAQkuQEO6cWDj3mwfOOR179hba0PjwakICDJDXBIn9nv+cqtuUO/\n0qMNIAUBSW6AQzr+gvab9Sf0aANIQUCSG+CQDv5q+82pB/RoA0hBQJIb6JCmt9+cMqJHG0AK\nApLcAId04ofbb75vVI82gBQEJDkLpJaWFse1dLsOSP88tHqRM7Mv9WgDSEFAkjNAamkreu6U\nth9FfGeeN+5/bPnH4p9yRduZ16aPGr7/W6ftaP+YN3X/ojsgrXzdET8pzW6aMSzyY7u6CkhB\nQJLTIbW0dCHplIvml2rO86+OPuoX5XsqkD478t4Vy+865MvtH7O4+xcdfEN29hHZ4WfVnzE8\ne/3j3f/1pYAUBCQ5GVJLS1eSqmzyHcd97zMfCu858drS29/9NviYbhe+RGjldW8bmu076npe\nIiQEJLk9BWnmwRvm7rMouOeiNz+588d0u51etNr6amtPF4DUKSDJ7SlIH7o0z0/+QnDPmguG\nvvGi7zSW7tnv4FL/0v2L5tXf5j0gyaWGNGRoqScX7vNwnn/r6G3h48+an1/ztwf8a/vfkXrw\nZ98B6dth3R/IgdQpIMmlhnThM6U2T84OPfTQ4dm/7fxE7qqDt0tP7bKwHm0AKQhIcnvmqd3W\nI6csKvrI6PZ7Fp9X/ird3fu8KkG6J6xHG0AKApKcEVLsbAXJXfuX/5t6aMiLlSdyK14bNepX\nixb/8k0fbH9qN39bty+avyOZ94Akl/r7SG2QTq+8SPuEz1W+/TouX3P1Ww7a/62Tm9u/IZs9\nHxuItiuk9XOWdvtXtwWkICDJOV4iVPuFDYkKIK2+deo9O/J/PSgb8umebQApCEhyA/xFq4te\nXzyUjW4Ydva178ju7tEGkIKAJDfAIV2y7/SH7zhq9Adb823/8/092gBSEJDkBjikN08s3vwo\n+37x9isje7QBpCAgyQ1wSAfdUrxZkM0q3n5nnx5tACkISHIDHFJ2e/GmKbu/eHs735DtdUCS\nA1IkIKkBKWhAQLru0UcfnZV9s3h7HZB6HZDkBjokXmvnCEhyAxzSlLAebQApCEhyAxxS7wNS\nEJDkgBQJSGpACgKSKyCpASlJQDLvAUkOSJGApAakoB5AqjNcSQ8CknkPSHJAigQkNSAFdR9S\nXd2elQQk8x6Q5ByQ6ur2sCQgmfeAJJcW0rjKi3cubvvBDEf87z8Wdx7X9hKEY6b35ufnlwOS\neQ9IcgZIdXU1Ja2YP/+e7Hfz5zfmp1yydOnSx88fvrAzpB7//PxyQDLvAUkuLaSiJ9p+OlDb\njxPaPvz/dobU4x9pVw5I5j0gyemQ6uq6lNQJ0o5DvwEkQ0BS63+Q6noAqfma/Xd6atfjn59f\nDkjmPSDJGSHFJVUhlchko36Xx/6O1NM/diCZ94Akp0Kqq9uNpCqkgsxTh5Ufet5yVentjtLP\nLeGpXe8CktrAhVQi88NhzxZv608p/f+C/Wc2G0i9DUhq/Q1SXd3uJHX6YsOHR23J82eHX/To\ns3e+YULem5+fXw5I5j0gye1ZSKuOurJ4+6cPv+Ggt9+8Je/Nz88vByTzHpDk+GcUkYCkBqQg\nILkCkhqQkgQk8x6Q5IAUCUhqQAoCkisgqQEpSUAy7wFJDkiRgKQGpCAguQKSGpCSBCTzHpDk\ngBQJSGpACgKSKyCpASlJQDLvAUkOSJGApAakICC5ApIakJIEJPMekOSAFAlIakAKApIrIKkB\nKUlAMu8BSQ5IkYCkBqQgILkCkhqQkgQk8x6Q5IAUCUhqQAoCkisgqQEpSUAy7wFJDkiRgKQG\npCAguQKSGpCSBCTzHpDk9lpIsxuba9bFqV7VuNI9aN5bs9w86L7A9UvNg03rzYOLa59a9pjh\n8zVJPCKZ93hEkttrH5GAFAQkOSBFApIakIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQk\nAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCC\ngOQKSGpAShKQzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIakIKA5ApIakBKEpDMe0CS\nA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQG\npCQBybwHJDkgRQKSGpCCgOQKSGpAShKQzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIa\nkIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7\nQJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCCgOQKSGpAShKQzHtAkgNSJCCpASkISK6A\npAakJAHJvAckOSBFApIakIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUC\nkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCCgOQKSGpAShKQ\nzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIakIKA5ApIakBKEpDMe0CSA1IkIKkBKQhI\nroCkBqQkAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQGpCQBybwHJDkg\nRQKSGpCCgOQKSGpAShKQzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIakIKA5ApIakBK\nEpDMe0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7QJIb1JAeGPNo\nnm+4+eLxUxs7jpWAFAQkucEMad3EcQWkGyYvfPmmK3a0HysBKQhIcoMZ0vQ7Jj6aN41dUDwa\nnTOveqyeBFIQkOQGMaQ5l20uIM0Z11rc/tTPqsf2s0DqCEhygxfShklz8wLSrEtK73zxtuqx\neLPqwaLfz19Ss4bap3pVg33QvLe4v1/gkkXmvcWLzYNdXODzcxJakNoNpNn19fXP5bfckpch\nXVq6q4BUOVY/iEekIB6R5AbjI9LGhoaGLXMnNZch/bHtKd3d1WP1g4AUBCS5wQip3I3jxo8f\nP/b8aWvGzi8+k+v/Uj1WzwMpCEhygxZSc1PRRfe9kk+/auGy669ubT9WAlIQkOQGLaRyxVO7\nfOOMSROmre04VgJSEJDkBjekLgJSEJDkgBQJSGpACgKSKyCpASlJQDLvAUkOSJGApAakICC5\nApIakJIEJPMekOSAFAlIakAKApIrIKkBKUlAMu8BSQ5IkYCkBqQgILkCkhqQkgQk8x6Q5IAU\nCUhqQAoCkisgqQEpSUAy7wFJDkiRgKQGpCAguQKSGpCSBCTzHpDkgBQJSGpACgKSKyCpASlJ\nQDLvAUkOSJGApAakICC5ApIakJIEJPMekOSAFAlIakAKApIrIKkBKUlAMu8BSQ5IkYCkBqQg\nILkCkhqQkgQk8x6Q5IAUCUhqQAoCkisgqQEpSUAy7wFJDkiRgKQGpCAguQKSGpCSBCTzHpDk\ngBQJSGpACgKSKyCpASlJQDLvAUkOSJHMnwYtRd5FIKkBqdSAgtRSybkJJDUglRpIkFpaEkgC\nkhqQSgHJuFUKSHJAimT8NGhpSSEJSGpAKgUk31Q5IMkBKRKQ1IAUBCRDQLIEpCQNIEh8scES\nkJIEJONWKSDJASkS35BVA1IQkEzxEiE5ICVpgEHiRatyQEoSkMx7QJIDUiQgqQEpCEiugKQG\npCQBybwHJDkgRQKSGpCCgOQKSGpAShKQzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIa\nkIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7\nQJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCCgOQKSGpAShKQzHtAkgNSJCCpASkISK6A\npAakJAHJvAckOSBFApIakIKA5ApIakBKEpDMe0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUC\nkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCEiugKQGpCQBybwHJDkgRQKSGpCCgOQKSGpAShKQ\nzHtAkgNSJCCpASkISK6ApAakJAHJvAckOSBFApIakIIGN6TZK9bXbGXtU71qZRe/V+8GzXtN\nL5sH3Re4dol5cNVa82BD7VNLHjN8viaJRyTzHo9IcnvtIxKQgoAkB6RIQFIDUhCQXAFJDUhJ\nApJ5D0hyQIoEJDUgBQHJFZDUgJQkIJn3gCQHpEhAUgNSEJBcAUkNSEkCknkPSHJAigQkNSAF\nAckVkNSAlCQgmfeAJAekSEBSA1IQkFwBSQ1ISQKSeQ9IckCKBCQ1IAUByRWQ1ICUJCCZ94Ak\nB6RIQFIDUhCQXAFJDUhJApJ5D0hyQIoEJDUgBQHJFZDUgJQkIJn3gCQHpEhAUgNSEJBcAUkN\nSEkCknkPSHJAigQkNSAFAckVkNSAlCQgmfeAJAekSEBSA1IQkFwBSQ1ISQKSeQ9IckCKBCS1\nvQ1SS6laJ4HkCkhq/RxSS0tXkoDkCkhq/RtSS0uXkoDkCkhq/RpSS0vXkoDkCkhqQEoSkMx7\nQFICUq2ApAakICC5ApIakJIEJPMekJSAVCsgqfV3SF18+7R3a3z5OxaQ1Po3pK6/f2ofBJIr\nIKlZIe3mAaS3mzXPAckVkNTSQHJK4kWrkYCk1p8hBV8bMEoCUiQgqQEpCEiugKQGpCQBybwH\nJDkgRQKSmvkCrZ/0SRwBKRaQ1KwX6P60B1I1IJn3+jUk/+d9AkdAigUktTSQvJJcY20BKRKQ\n1Po5JH6KUFtAMu/1Z0hpvsgGpFJAMu8BSQ5IkYCkxlO7ICC52tsgpfpqNZDMAcm854WU7Ps+\nrsEcSG0BybxnhZTs+z6uuVJAKgUk854TUopnYv5v+wCpFJDMe4kg+T75+/UPPykHpEhAEgLS\nzgHJFZDUgJQkIJn3+vvfkYCUJiCZ9/r5V+2AlCggeee8n/Qpvu0DpCQByTmW5IvV5q9XAylJ\nQDJupXgmlru/7QOkJO31kOqMW0m+ONCff/hJOSCVApJPUpKvVgMpDEiu3JDqgKQGpFJ7OaS6\nOqMkIHkCUiQgqQEpCEiuvJDq6lJJck0CqVNAcjVAILkWcyB1CkiurJDq6hJJcu2VAlIQkFyl\ngeT9yp1trByQgoDkygmpLgGk/v1ThEoBKUl7MaS6uhSSgCQHpEhAUgNSEJBc+SDV1SWRBCQ5\nIEUaMJBckoAkB6RI/RdSpX79L2RLASkISK6ApAakJAHJvAckOSBFApIakIKA5ApIakBKEpDM\ne0CSA1IkIKkBKQhIroCkBqQkAcm8ByQ5IEUCkhqQgoDkCkhqQEoSkMx7QJIDUiQgqQEpCK7D\nq6AAAAiWSURBVEiugKQGpCR1C9KnxxSdl+cbbr54/NTGjmMlIAUBSW7wQrp0ZlNT05o8v2Hy\nwpdvumJH+7ESkIKAJDd4IX3kifKhaeyC4tHonHnVY/U0kIKAJDdoIW0b880rPzZtWT5nXGvx\n3qd+Vj1WzwMpCEhygxbS+onfeOGF6ye+OuuS0ntfvK16LN60bit6pHlHzdbVPtWr1q81D7ov\ncOMq86D7ArcuNw9u2GoeXFL71OqBCml2fX39c+Vbm867b9alpRsFpMqxeLPqwaLfz19Ss4ba\np3pVg33QvLe4v1/gkkXmvcWLzYNdXODzc9Jy6H27gbSxoaFhS9vNT/74j21P6e6uHqsfxFO7\nIJ7ayQ3ap3YN3yr+qDaf9/s1Y+cXn8n1f6keq+eBFAQkuUELqXn8jBXLpl26JZ9+1cJl11/d\n2n6sBKQgIMkNWkj5gi999KIbVhZP9GZMmjBtbcexEpCCgCQ3eCHtJiAFAUkOSJGApAakICC5\nApIakJIEJPMekOSAFAlIakAKGuSQnnupZn+pfapXPfuMedB9gc//yTzovsAX55oHn3vRPPhU\n7VPPDGpIKxfXbkEX53rTM382D7ov8L+fMg+6L/Clx8yDixrMg490cW6V4fM1SQ5Ie7JF8/v6\nCnbTqnm7/5g+beMjfX0Fu+v+vr6A3gQkc0CSA9IeCEhqQEoSkMwBSQ5IRHtrQCIyBCQiQ0Ai\nMjSwIK256aLzr32hr6+i6x4Y82hfX0JX/fqycz/1eF9fRBct/cqEj177bF9fRc8bWJA+O3nB\n8q9P2NzXl9FV6yaO68+QHpj0ROMvLt/Y15dRs9bLv7Vxy4/Oa+7r6+hxAwpS87Qleb5qzIt9\nfR1dNf2Oif0Z0uUP9vUVdN36Mc/n+dox/fxZR6QBBanc8/Vrd/9Bfdacyzb3Z0irxzz46Y98\n7vm+vowu+vyM5s0/vmxrX19GjxtwkJo/eWdfX0IXbZg0N+/PkF4Yc93S5tsucP/TDGNrrhgz\nZtJLfX0VPW+gQVr68Vtbd/9RfdYtt+T9HNK8PH/twgf6+jpqtv3Kb63fePeE/vykI94AgzRv\n/My+voSumjupuX9DahpTeonVFXfv9gP7qqfHlr6U9LFf9fV19LiBBenZC/vtP+wqd+O48ePH\njz1/Wl9fR812TCr+h2jr+X/o6+uo2VNjSl9RnASkpG29/CdNRf33y9/Npcu76L5X+vo6anf3\nhLlN35zUf/8EN0761oat/zHO/c/h0zegIM0bU+7evr6OruvPT+3yHT+YeO61XfxQhD6vYeqE\nC/75z319FT1vQEEi6q8BicgQkIgMAYnIEJCIDAGJyBCQiAwBicgQkLrblOzdbS+XPWV0Fx91\n2vHdnNs+8aBhS0s33r9f9cfwHn/kzv964F2d1jq2u/270J4KSN1tSpZ9p3zDA+nebMKvyv9Q\n9e7s5rZ75mTX7PxBMzq9ag9I/TggdbcpB37osPJjhwfS97PKK0e3HXVi243Ls//u+pcAqR8H\npO42JXvxwItLN0qQTjqpdKv+8Dw/4/Q/vPPAkTdumzxy+OgFxaf4CU+dftBhk0r/L07/dfaI\nYSffUdw47YyZx76nsvKbM4Yf+PabW/PRWdGi8l2fz8o/jGTTIe/N85+8c9iIU37S/mvKT+06\n7mvfLkNq319+2V8fcPSH+/M/e90LAlJ3m5JtmZo9lO8MafSxZz619Nzs7KnLHjrkH4pP8WOP\nv/Gezw8Zk+cPDH3vzPv+Kft6np816oR/qbzQ9p4hH/zFA1dnn89fmJJ994m2vxK9mP1T6fDD\n7If5T7Nz7733g9m91V9TgtRxX8d2CVLH/rtf/93f3/V3R/Xfn2iyNwSk7jYl27zlLW/btguk\nbF6ez85OLd6dcHDxKZ79e3FrfLY4P/nNpc/ssSM2Fx/y8+rICX9dwnPOfqvzO7PZ1Tvfd+im\n4u1Zh23Op51VnH5l3wl55deUIHXc17FdgtS+/0p2bXHjpWkv78E/DNo5IHW3AlL+u+xru0Aq\n9OQvlb9OcE3WnJ92QEGtUPLzxuzKzUX/r3jeNnr/bZWNl9sefO4oHl8CSHdld+V5w5DPVN8/\n9oy88ms6vmpXuq99uwSpY3/b4cc9sCP9f37qMiB1txKk/LyDGnaGdFxxY1E2vXg7OVuXn/bG\n0onfZt+Zm1X6eT56ZHXj8eyG0uE32W0hpC2vKxanZs8UDzxfPvGQoUOz0/LKrylB6rivfbsE\nKdh/+I3Z4ePucv8/uVKPAlJ3K0NaNmJs/o4uIb2pdOI32e1zs489Wq6p7UPKPZFNLR1+nX03\nhJRfOaSh9W/eXdx479Av/OHPz4w8rTJbhtRxX/t2G6T2/fy1B695W/aOTXvqT4IiAam7lSHl\nN2e/fE8B6eTyl6zfFYF0cOlZ1veyX63JLq7+yg5IK7KPlw63ZbM6QXo2mzY7+16ez88uL97b\nfmAIKbivfbsEKdgvd2v2/TT/ualbAam7tUHa/nfHnVpAOuuI1jxvHBaBlN1X3DpnnxX53x9a\n+hr4D764PYCUnziyNPLBg17pBCk/9ZQrD3k1z58rP2B9M3t3ACm4r2O79MWG9v0nP9qYl/6a\ndtOe+XOgaEDqbm2Q8oeHDCkg3ZJNX/n0mW/fFdJ7jn3rtx+YnF2Y5/+136gf/O5L+10SPiLl\nv97nA7/87SdKH90J0p3ZYaWvQmz7q2N++fDn3ve+Eb9/tR1ScF/HdglS+/6KEaPuuP+npx4y\nAH+q4iAKSN2tAim/NCsgbb36mANOmnnFiF0g/a/3PHn6sMMu21C8O/v9I/Z7643bO0HK7zv9\n4ANOLp7FdYa08dDsqdLxifccdPT/eWXmEYe90PF3pI77OrbL35Bt3//TuUftN/Lcp/fMHwPF\nAxKRISARGQISkSEgERkCEpEhIBEZAhKRISARGQISkSEgERkCEpEhIBEZ+v/5w8qcIA/SqAAA\nAABJRU5ErkJggg=="
          },
          "metadata": {
            "image/png": {
              "width": 420,
              "height": 420
            }
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# adjusted R2\n",
        "adjr2.df <- data.frame(numvar = numvar, adjr2 <- reg.summary$adjr2, maxadjr2 <- allfalse)\n",
        "adjr2.df$maxadjr2[which.max(reg.summary$adjr2)] <- TRUE\n",
        "ggplot(adjr2.df,aes(x=numvar,y=adjr2,shape=maxadjr2,col=maxadjr2)) +\n",
        "    geom_point(size=3) + theme_light() +\n",
        "    labs(x = \"Number of Variables\", y = 'Adj'~R^2, color='Maximum Adj'~R^2, shape='Maximum Adj'~R^2)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 437
        },
        "id": "iAyTYt15m0NL",
        "outputId": "465d1894-88f1-4596-f82d-24beaf06f7c1"
      },
      "id": "iAyTYt15m0NL",
      "execution_count": 83,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "plot without title"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAAAv8QBAQECAgID\nAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQV\nFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYn\nJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5\nOTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpL\nS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxd\nXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5v\nb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCSkpKTk5OUlJSV\nlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqan\np6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5\nubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrL\ny8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd\n3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v\n7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4dm34+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///8t\n+x4CAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3dC5xVdb338TWIgIiIPlkdyfQx\nS+vkJR/0Ubx0SqtTnTOAIWkM4iW1U1iYl7A0JM4DpRnGyUorTxczozRRU0F7Ei94QylRQy5y\nGS7DMMgMjNxG5n/W2ntuDHv/2T/4MjM/+XxeL9fes9aeH2sW++1es2fPJglEtMslnb0DRG+H\ngEQkCEhEgoBEJAhIRIKARCQISESCgEQkSAnpjUeL9kjxTdam6UYJ9+qR6bpZylHCL1G4W5G/\nQ+H9sUNTQlr9QrEtb63T/SlLdKPqtspGbVglGxXW6katER544W5VFj3wj+j+kI4NSJqAZAlI\n0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiW9mxIlVcObr667qbzh4+rar3MByRNQHJY6ZBmjJzU\nAmn8mIXLbhy1teUyH5A0AclhpUN6dNXMZkjVgxakj0ZDZjdfNq0GkiYgOczyPVILpKeGNqbL\ny+5qvmxaDSRNQHLYTkF66IJsee2tzZfp4m8zZ8587skVRVpeWWyLvdd1oyqXy0YtWywbtWKp\ncJTwwAt3a1HRLQ+L798d1s5BujBbppCaLtPFhvr6+mWzGoq0qbbYFntLdKPWbpaNWl8lG9Ww\nRjdq9VrdLOFuLS164PesR6Sn86d0U5ovm1ZzaqeJUzuH7RSkmkHzQqgdPKf5smk1kDQByWGl\nQ1pTPW1wdfWGMG1qCBMvX1h5/RWNLZf5gKQJSA4rHdIXy7PuDTdcF0L9pJEVE9a0XuYDkiYg\nOYyXCGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQw\nIGkCkiUgRQOSJiA5DEiaNpTJRgHJY0DSVAYkQ0CKtgdDKisTSgKSw/ZgSBtzaWaV6SQJ90r5\nFQpHhcgoIAVvkDZu1N3PymSQNgp3q2uOCtFRQArOIG1U3mXLZJKEeyX8CpUHKz4KSMEXpI3K\n+0aZDJJwt5RfYcftFpCCW0i7fN8oK5NJ2j332F2dpTxYQNpxQFLulnAUkHZ/QNrl+0ZZmU7S\nHgWp0CwgBbeQdnVUmVBS14TE90g7bk+FpHsiqqyLQto9z9oBqUhA2tVBZbtH0q7u1u75kRRP\nfxdrj4UU/7mgOdmrv4W7pfwKO2q3gBS8QQrZX6lslPDXKHS6laOEs2IkgRQcQuL3kSzxotVY\nQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIE\nJEtAigYkTUBymDdIyheaAskUkGI5g6R9yTaQLAEpli9Iyl+xCUCyBaRYriApf380C0iWgBQL\nSJqAZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBSNp781AclhviDxA1ljQOqonEHiJUK2\ngNRReYPEi1ZNAamjApImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqA\nZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqAZAlI\n0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqAZAlI0YCk\nCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqAZAlI0YCkCUgO\nA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDlMCWnVM3VFWltdbIu9xbpR\nq2plo2qWyUbVVelGrRQeeOFuLSl64KcL748dGo9ImnhEssQjUjQgaQKSw4CkCUiWgBQNSJqA\n5DAgaQKSJSBFA5ImIDkMSJqAZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAg\naQKSJSBFA5ImIDkMSJqAZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKS\nJSBFA5ImIDkMSJqAZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBF\nA5ImIDkMSJqAZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5Im\nIDkMSJqAZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkM\nSJqAZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqA\nZAlI0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqAZAlI\n0YCkCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqAZAlI0YCk\nCUgOA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqAZAlI0YCkCUgO\nA5ImIFkCUjQgaQKSw4CkCUiWgBQNSJqA5DAgaQKSJSBFA5ImIDkMSJqAZGmPhrTupvOHj6vK\nX18xYcSw76YH9qvlacOabwEkTUByWOmQxo9ZuOzGUbkjsOVL4ysXXffNEC68r7q6uqb5FkDS\nBCSHlQypetCC9FFpyOzs+tzy1emK8kXh7Ofa3gRImoDksJIhPTW0MV1edld2fU55XYpjyCNb\nyn84+qIJlc03AZImIDmsZEgPXZAtr701W75Z8ZOGhjuG3LP2vB/MnXv9eevTVauWL18+/9n1\nRapbXWyLvcW6UdXrZKPeWC4btX6VblSV8MALd2tJ0QM/fTfdz3d7pUO6MFvmIYWXLh1yzh2X\nTs1df3PYtHS58JVXXpnzdE2RVq8otsXeIt2oFatlo6oqZaNqlgtHCQ+8cLcWFz3w03bDfbxD\nKhnS0/lTuylNH65vaDhrZv7qV37btI5TO02c2jmsZEg1g+aFUDt4Tnb9rRlrQnh2SO2iyQ3p\nXWjYX5puAiRNQHJY6U9/T7x8YeX1VzSGaekZ3egJ1XPOvyXUDZ+0onLChZuabgEkTUByWOmQ\n6ieNrJiQPhLdcF0Iy7519ojb0kejBdedM2L8yuZbAEkTkBzGS4Q0AckSkKIBSROQHAYkTUCy\nBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRo\nQNIEJIcBSROQLAEpGpA0AclhQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIE\nJIcBSROQLAEpGpA0AclhQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcB\nSROQLAEpGpA0AclhQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQ\nLAEpGpA0AclhQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEp\nGpA0AclhQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0\nAclhQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0Aclh\nQNIEJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIE\nJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtA\nigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtAirbq\nmboira0utsXeYt2oVbWyUTXLZKPqqnSjVgoPvHC3lhQ98NOF98cOTQmpetbWIm2pK7bF3hLd\nqNoG2aj6KtmorW/oRtUID7xwtyqLHngekQKndqo4tXMYkDQByRKQogFJE5AcBiRNQLIEpGhA\n0gQkhwFJE5AsASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQk\nhwFJE5AsASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJ\nE5AsASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5As\nASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASka\nkDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASkakDQB\nyWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASkakDQBqQs3\n44QD3veLAuuBpAlIlvxCWr3/78Lj+7y2/QYgaQKSJb+QVv4qXRx+//YbgKQJSJb8Qsp6uc/K\n7VcCSROQLLmGtODIXxVYCyRNQLLkGdKTh95VaDWQNAHJ0k5CGpu8c0v+2heTUwve4qSjzPtS\n2ye5u82H5/TZbsrYJKvHEUOfTT944r0zC44BkiYgWdpZSHt1vyd35c39exSGNGmCeV9uSQ78\nbJsPM0jtpoxNvnnbbbdNHtWv5+Nh/aFPFB4DJE1AsrSzkPY+ZVDuyh09TigMaSc6/vjLuy9t\n/TCD1P6PTfIPQrO7/Wu4M+mVNnH7MUDSBCRLO31qN2Hv3BNmn/rUqRmkO0/s3XfAnSFMTb6f\nfZRMzp2UnX7ajBP36X/DljH99ztzQQjHHZd9xuB3bLch33PJ955PxmfXGscd0uvoKQVP7ZrO\n5g49ovi+AUkTkCztNKR5e2VkKvf6xckppN8lZ91//6eT+0M4v8/iUHvwxxtzBM485OOzlp6V\nfGJc5WP7/1tbSO025Lu0+7Jw7OGN6bXvJRXT7zr6qOKQVvX8ZPF9A5ImIFnaaUgbPvHh9PK7\nvetOSiFNOGNzCLV7V6S79t5BYVTfRXkCZyazQ3g8OSW9YUWfbSBtuyHXur7pN0g3J9PTB6T+\nR6cfL+9RCNIDK1asWPTgid0eKL5vQNIEJEs7D+k3ybMhfOgL4aSW75EOOT1dTEuu3eu20AQp\nQzI/uSpdXpXUtYW07YZctyV/SO+4Pc8NYXHytWzFwGLP2iXJUXdG9g1ImoBkaech1ff9cng2\neTAHqfbbR+/fvXv+ifBLkk9kFzlIh6VXXk+yJwTGJG+0hbTthlwn9lteXV39771qwjP575SG\nFoI06cEHH/yXfZaGSEDSBCRLOw8pXHTAhlEHv5WD9NHu35rx95f6Z5C2ntbtsOwxxgppdtOD\nTXJzeDoPaUix75Fe6n52bN+ApAlIlnYB0mPJvQddGTJI85JL0nUN+2SQvr/3fftnHxWAdHz2\nrU84qTCkryR3Ts/qf0xYkIzK1nyk6JMNX04ejuxbG0hz7889tVjgOfISA5ImIBUpg9T4vhOS\nv+UgvZKMS9f9MDk5hH/0HhP+K3moIKQzDmoMoap3QUhv9huYn3xN8szWg45I92put6KQqg84\ncnPxfWuFdFuvw/tkr8brFftSogFJE5CKlEEK1yfZuVoKact733PvE1d+7GN9/1J78hFvhq0D\n3/NGIUg3JxNXvvDxDxeEdHvys/zk19JHt28nn/vjjw8bUPzp70nJfxbft1ZI75sZ5nzgD0Da\nuYBkaVcgLex2U8hBCs8N3PfdX6q976ADv5g8mq6a02NkIUibr3hPr+PuG9W3EKRT+jQfn4/2\nXf/WNf/U85h7LutZFNKWo3q/XnTfWiEdkP43/7BngbRTAclSF3/190689LUNpAG/TxdPHvYE\nkHYmIFnq4pA+coz9c1oh/bXvr9PlE4fu/PN4QNIEpE5t/u09hto/qw2b6tyTdutiP76NByRN\nQOrUbt7rmDn2z4o9/hT/zqpwQNIEJIdtB+mxT73/U3/Ormz6z97GWUDSBCSHtYc0s0e3Q3t0\n+30ID38gsT53ASRNQHJYe0iD+80OVSd8aOnZyQGTthhnAUkTkBzWHtL/Hp0uHkz26f7lavMs\nIGkCksPaQ9r7R+liUfIvL+3ELCBpApLD2kNKsl+PWpE8uDOzgKQJSA4DkiYgWQJSNCBpApLD\ntoP0zZkzZz6QTJqZZZwFJE1Acth2kNq2zZZ1N50/fFxV/vqKCSOGfXfttuuApApIDmsPaWzb\nttkyfszCZTeOyh2BLV8aX7noum9usy4ASRWQHFbya72rBy1IH4GGzM6uzy1fna4oX9R2XQCS\nKiA5rGRITw3N3o3ystw/aTGnvC7FMeSRtusCkFQByWElQ3rogmx57a3Z8s2KnzQ03DHknjbr\nnnzooYcef3xRsV4vusVc1xzVRXfr9a65WwuLbnlod9zJO6LSIV2YLfOQwkuXDjnnjkunbrOO\nRyRVPCI5rGRIT+dP46Y0fbi+oeGsme3WAUkTkBxWMqSaQfNCqB2c+93Bt2asCeHZIbVt1wUg\nqQKSw0p/h4aJly+svP6KxjBtagijJ1TPOf+W1nX5gKQJSCW0ceNGzb6IKh1S/aSRFRPSR6Ib\nrgth2bfOHnFbQ+u6fEDSBKQdtjFfwW0D8q8muD2Eqp6HvJVbMyq/5a2Jx+7X88gJW1tuE/mH\nw6zx3t+agGRpVyFt3BiRNGDEvKy6EP7fme/6U25NE6Sv979/xfI79v92y20Wm/e8aEDSBCRL\nMkiFJDWzCVsP+8XXPtN2zdHXZMuHH2xzG1lA0gQkS7sIaWNpkO7rs+7FvV5vs2bE+59vfxtZ\nQNIEJEs6SAUktSD5zIUhHP+tNmtqzu1++IifZi+yHtCjT9aPjPsdCUiagGRp90Lq1j3r+YV7\nPRHC5Hdvafv4U3P3VR/q9auW75GEX1J7SEdNCEe1dszZrxhmAUkTkHbQDiB94aWsDWOSfv36\n7Zf8vv2J3OV9Gjrg1O6kSeGk1j7S4yOGWUDSBKQdVNL3SJvfOfb1tLPPbFmzeFjuWbope63v\n+O+Rfr6XYRaQNAFpR5UC6Y6euTeUe6zba00nciveOvbYqa8vvveIT7ec2s2zvnVj8XYA6bWb\nDbOApAlIOyz6c6Q8pNPOzX/4wSubfvw6NNRc8YF9ex45pq7lB7LJq8YdL157SH3a1NM4C0ia\ngLTjir+woZNqD+mctKN6DBw65CPdBlxmnAUkTUByWIFTuylHL88u/vHBqcZZQNIEJIcVgHT0\n7/OXPz7OOAtImoDksAKQejZ9MVOs/5oskDQByWEFIPWvyF00nnOwcRaQNAHJYQUgjU2O+dr4\n8aM+lFxjnAUkTUByWAFIjd87OHuK/aBvv2WcBSRNQHJYwR/INi5+5ukFW8N64ywgaQKSw4q/\nsuGZS/oaZwFJE5AcVgRSzc3HJMnpxllA0gQkhxWC1Dj93F5J/2tes84CkiYgOWx7SEvHH570\n+vdkun0WkDQByWHtId392e7JsTevrgaSLSBZ2gMgJQd+PXuDCCAZA5KlPQBSn+T47y0Dkjkg\nWVJAKpPsiaz2kGpvOT7p/tk/LAeSLSBZkkDqWpIKPGv33KV9kwOSO+2zgKQJSCVU1vUhhbDu\nthOTZODPrH8JQNIEpBIq62KSir2yYfZX+iV9jLOApAlIO66srJikoU1vx3B+/o0ZDvrXp9OV\nh43NbXvPxN3x/vm5ir9EqP72gcZZQNIEpB1WVhzSinnz7kkenjevKgy4YOnSpc9+fr+F20KS\nv39+Lt5pVROQLKkgFTm5ey7/7kD5txNq2O+/toUkf0u7XEDSBCRLuwqpzABpa78fAKlQQDL0\ndodUWFJbSHVX9Wx3aid///xcQNIEJEu7CKmsREgZmeTYh0Oh75GEX00uIGkCkiUdpIKSmiGl\nZGYdmHvo+cDl2XLrvjdzatcSkAy9LSGVle1AUttTu1/3fjldDh6Q/Yvh/z95HEgtAcnQ2xFS\nmQlS+Nyxm0J4eb8RM1++/eDs/bH075+fC0iagGRJCamApG0grXrX6HT5t88dvO+Hb9oUdsf7\n5+cCkiYgWdoDfo1iVwKSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUia\ngOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQw\nIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQwIGkC\nkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQwIGkCkiUg\nRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQwIGkCkiUgRQOS\nJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQwIGkCkiUgRat6Zm2R\n1lQV22JvkW5U1RuyUauXyUatXakbtUJ44IW7tbjogZ8mvD92aDwiaeIRyRKPSNGApAlIDgOS\nJiBZAlI0IGkCksOApAlIloAUDUiagOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZ\nAlI0IGkCksOApAlIloAUDUiagOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0\nIGkCksOApAlIloAUDUiagOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkC\nksOApAlIloAUDUiagOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOA\npAlIloAUDUiagOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlI\nloAUDUiagOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAU\nDUiagOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUia\ngOQwIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQw\nIGkCkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQwIGkC\nkiUgRQOSJiA5DEiagGQJSNGApAlIDgOSJiBZAlI0IGkCksOApAlIloAUDUiagOQwIGkCkqU9\nGtK6m84fPq4qf33pdyrOueblEL5anjas+RZA0gQkh5UOafyYhctuHJU7Ao2XTK7f9JthdeHC\n+6qrq2uabwEkTUByWMmQqgctSB+VhszOrq8tfzX9OyqfG85+ru1NgKQJSA4rGdJTQxvT5WV3\n5T64elLdht9evHlL+Q9HXzShsvkmQNIEJIeVDOmhC7LltbfmPqgZVV4+cn5Ye94P5s69/rz1\n6Zol8+fPf+WZtUVas6rYFnuLdKOq3pCNWr1MNmrtSt2oFVW6WcLdWlz0wE/bPXfz3V/pkC7M\nlnlIDaMnr62fUrEmt+HNYdkXv3zx4sVzn6kr0trqYlvsLdaNWlUrG1WzTDaqrko3aqXwwAt3\na0nRAz99d9zJO6KSIT2dP7Wbkl1/YdCGdHnR1PyWr/y26Sac2mni1M5hJUOqGTQvhNrBc7Lr\ns8rr0+XIqYsmN6R3oWF/aboJkDQByWGlP/098fKFlddf0RimTQ31Iyev2/zHocvrhk9aUTnh\nwk1NtwCSJiA5rHRI9ZNGVkxIvy264boQFo2rOPcbfw9hwXXnjBi/svkWQNIEJIfxEiFNQLIE\npGhA0gQkhwFJE5AsASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA\n0gQkhwFJE5AsASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQk\nhwFJE5AsASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJ\nE5AsASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5As\nASkakDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASka\nkDQByWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASkakDQB\nyWFA0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASkakDQByWFA\n0gQkS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASkakDQByWFA0gQk\nS0CKBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASkakDQByWFA0gQkS0CK\nBiRNQHIYkDQByRKQogFJE5AcBiRNQLIEpGhA0gQkhwFJE5AsASkakDQByWFA0gQkS0CKBiRN\nQHIYkDQByRKQoq16pq5Ia6uLbbG3WDdqVa1sVM0y2ai6Kt2olcIDL9ytJUUP/HTh/bFD4xFJ\nE49IlnhEigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIE\nJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtA\nigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtAigYk\nTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtAigYkTUBy\nGJA0AckSkKIBSROQHAYkTeQj5HYAAAonSURBVECyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIE\nJEtAigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtA\nigYkTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtAigYk\nTUByGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtAigYkTUBy\nGJA0AckSkKIBSROQHAYkTUCyBKRoQNIEJIcBSROQLAEpGpA0AclhQNIEJEtAigYkTUByGJA0\nAckSkKIBSROQHAYkTUCyBKRoRSFtTJP9KUAyBKSOqgMgbWxK9KcAyRCQOqrdD2njRq0kIBkC\nUkcFJE1AsgSkaAUhbdwolgQkQ0DqqICkCUiWgBQNSJqA5DAgaQKSJSBF48kGTUByGJA0AcnS\nHg1p3U3nDx9Xlb++9DsV51zz8rbr+IGsKiA5rHRI48csXHbjqNwRaLxkcv2m3wyra7su8BIh\nVUByWMmQqgctSB+BhszOrq8tfzX9Oyqf23Zd4EWrqoDksJIhPTW0MV1edlfug6sn1W347cWb\nt1kHJFVAcljJkB66IFtee2vug5pR5eUj57dd17Bly5aVs7YWaUtdsS32luhG1TbIRtVXyUZt\nfUM3qkZ44IW7VVn0wO8BkC7MlnlIDaMnr62fUrGmzbpZM2bMmPnksiJVLi22xd7rulFLdKOW\nLtbNEu7WEuUs3ajif4cP75Z7eQdUMqSn86dxU7LrLwzakC4vmtp2XeDUThWndg4rGVLNoHkh\n1A6ek12fVV6fLkdObbsuAEkVkBxW+tPfEy9fWHn9FY1h2tRQP3Lyus1/HLq8ZV0+IGkCksNK\nh1Q/aWTFhDUh3HBdCIvGVZz7jb+3rssHJE1AchhvfqIJSJaAFA1ImoDkMCBpApIlIEUDkiYg\nOQxImoBkCUjRgKQJSA4DkiYgWQJSNCBpApLDgKQJSJaAFA1ImoDkMCBpApIlIEUDkiYgOQxI\nmoBkCUjRgKQJSA4DkiYgWQJSNCBpApLDgKQJSJaAFA1ImoDkMCBpApIlIEUDkiYgOUwKacb8\nIr32SrEt9mbpRr08TzbqH3+TjZo/RzfqJeGBF+7WC0UP/IPC+2OHpoS0YXGxFi0susnck7pR\nCxfJRs19XjZq8QLdqNkv62YJd2tm8fuD8P7YoSkhdUzTO3sHClb8tLZTe7Wys/egYI9t7uw9\nkAckTUCyBKQuEJAMAamjApImIFkCEhEVCkhEgoBEJAhIRIKcQaq5ccTnr5nb2XtRqEfKZ3b2\nLmzfAxefddmznb0T27X0OxXnXPNyZ++FOGeQvj5mwfLvV2zo7N3YvjfOG9r1ID0y8rmqP11S\n39m70a7GSybXb/rNsLrO3g9tviDVTVgSwqry1zp7P7Zv4s/P63qQLnm0s/egUGvLXw1hTXmX\nPK/Y+XxByvXq4DU7vlEH99TFG7oepNXlj3717Ctf7ezd2K6rJ9Vt+O3Fb7MfJfmDVPeV2zt7\nF7Zr3cgXQ9eDNLf8m0vrbj1X+OsPmmpGlZePnN/ZeyHOHaSll97SuONbdXA33xy6JKTZIbz1\nha72Kz4NoyevrZ9S0fVOK3Ypb5BmD7+vs3dh+14cWdcVIVWXz0uXo6Z09n6064VB2ZNFF03t\n7P3Q5gzSy194vrN3oUA3DB0+fPigz0/o7P1o19aR6f90Nn9+RmfvR7tmlWfPI44EUie2+ZI7\nq9O62tPfddlOjZhW29n70b4pFS9W/3BkVzta9SMnr9v8x6HLO3s/tPmCNLs81/2dvR+F6nqn\ndmHrL8876xrhe1yIWjSu4txv/L2z90KcL0hEXTQgEQkCEpEgIBEJAhKRICARCQISkSAgEQkC\nUqyxycn5V8gOODNyq1OPKnFcw3n79l6aXflkj+b33D/qne1/neCkbaa1zi75T6HOCEixxibJ\nT3NXNJDuTyqm5n5fdUpyU37NU8lV7W80aZuX7AHJSUCKNXafzxyYe+zQQPrvpOkFpFvedXT+\nyiXJP+KfAiQnASnW2OS1fc7PrmSQjjsuuzb4HSGcftqME/fpf8OWMf33O3NBehf/4KzT9j1w\n5Bvp1r9+om/v43+eXjn19PsOGdg05c+n77fPh29qDGcmaa/nVl2d5N6T5M39PxrCnSf27jvg\nzpbPyZ3ata5rmZ2D1DJ/+cWH9nr357reb7/usQEp1thk07jksdAe0pmHfHzW0rOST4yrfGz/\nf0vv4occdcM9V3crD+GR7h+9b9p/JN8P4YxjP/ijptfW3tPt03965Irk6jB3bPKz5/LfEr2W\n/Ed28evk1+F3yVn33//p5P7mz8kgta5rnZ1Bap1/8j/97C93HPOurvbGJntuQIo1Ntmw6QP/\nvGU7SMnsEB5PTkk/rOiT3sWTP6TXhieLw/Hvz+7Zg/puSG9yd/OQDx6a4RnSY3W4PXm8eeXH\n+r2ZLs84cEOYcEa6uXbvitD0ORmk1nWtszNILfNrk2vSK/MnLOvAg0GxgBQrhRQeTr67HaRU\nT5ife57gqqQunNorpZYqubsqGb0h7SfpeduZPbc0zViWf/D5efr40gbSHckdISzq9rXmjw85\nPTR9Tuuzdtm6ltkZpNb5W95x2CO6f7eTdj0gxcoghWH7LmoP6bD0yuvJxHQ5JnkjnHp4tuHB\n5KcvJk3dHc7s3zzj2WR8dvHn5Na2kDb9r3TiuOSl9IHn20fv3717cmpo+pwMUuu6ltkZpDbz\nnzg8ecfQOxo67lBQPCDFykGq7DsonBCFdES24c/JbS8mF83MVZ2/Sa7nknHZxQPJz9pCCqO7\nLWp838nplY92/9aMv7/U/9SmsTlIretaZuchtcwPbz161T8nJ7zZUUeCdhCQYuUghZuSewem\nkI7PPWV9UgFIfbKzrF8kU2uS85s/sxXSiuTS7OLW5KFtIL2cTHg8+UUI85JL0o8a9mkLqc26\nltkZpDbzc92S/Pfu+brJHJBi5SE1HHPYKSmkMw5qDKGqdwFIybT02pC9VoT/2y97DvyX1za0\ngRSO7p8N+fS+tdtACqcMGL3/+hBeyT1g/TA5uQ2kNutaZ2dPNrTMf/6cqpB9m3ZjxxwH2mFA\nipWHFJ7o1i2FdHMyceULH//w9pAGHnLkjx8Zk3whhL/2OPaXD1/X44K2j0jhgb0+de+DX85u\nvQ2k25MDs2chtrz3Pfc+ceXHPtb3L+tbILVZ1zo7g9Qyf0XfY38+/Xen7P92e5tFvwEpVhOk\ncGGSQtp8xXt6HXffqL7bQfo/A58/rfeBF69LP3z8k317HHlDwzaQwrTT+vQ6Pj2L2xZSfb9k\nVnb53MB93/2l2vsOOnBu6/dIretaZ+d+INsy/29nvatH/7O65r+3uUcGJCJBQCISBCQiQUAi\nEgQkIkFAIhIEJCJBQCISBCQiQUAiEgQkIkFAIhL0PyLGZ2X3u2j5AAAAAElFTkSuQmCC"
          },
          "metadata": {
            "image/png": {
              "width": 420,
              "height": 420
            }
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# identify which model has the best bias-adjusted fit\n",
        "which.min(reg.summary$cp)\n",
        "which.min(reg.summary$bic)\n",
        "which.max(reg.summary$adjr2)\n",
        "\n",
        "# let's pull out the coefficients\n",
        "coef(regfit.full, 3)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 86
        },
        "id": "xl-7nOWqm19A",
        "outputId": "b62e61e7-b334-459b-c5d5-16f14baea6d3"
      },
      "id": "xl-7nOWqm19A",
      "execution_count": 84,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "3"
            ],
            "text/markdown": "3",
            "text/latex": "3",
            "text/plain": [
              "[1] 3"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "3"
            ],
            "text/markdown": "3",
            "text/latex": "3",
            "text/plain": [
              "[1] 3"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "4"
            ],
            "text/markdown": "4",
            "text/latex": "4",
            "text/plain": [
              "[1] 4"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<style>\n",
              ".dl-inline {width: auto; margin:0; padding: 0}\n",
              ".dl-inline>dt, .dl-inline>dd {float: none; width: auto; display: inline-block}\n",
              ".dl-inline>dt::after {content: \":\\0020\"; padding-right: .5ex}\n",
              ".dl-inline>dt:not(:first-of-type) {padding-left: .5ex}\n",
              "</style><dl class=dl-inline><dt>(Intercept)</dt><dd>-1.7786055643429</dd><dt>X</dt><dd>0.765488944992594</dd><dt>X2</dt><dd>1.99154367427184</dd><dt>X3</dt><dd>1.56719996937048</dd></dl>\n"
            ],
            "text/markdown": "(Intercept)\n:   -1.7786055643429X\n:   0.765488944992594X2\n:   1.99154367427184X3\n:   1.56719996937048\n\n",
            "text/latex": "\\begin{description*}\n\\item[(Intercept)] -1.7786055643429\n\\item[X] 0.765488944992594\n\\item[X2] 1.99154367427184\n\\item[X3] 1.56719996937048\n\\end{description*}\n",
            "text/plain": [
              "(Intercept)           X          X2          X3 \n",
              " -1.7786056   0.7654889   1.9915437   1.5672000 "
            ]
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e7f16d76"
      },
      "source": [
        "> *Write your response here*\n",
        "> In looking at model selected based on Mallow's Cp, BIC, and adjusted $R^2$, we first see that there is disagreement between adjusted $R^2$ and Cp and BIC. While Cp and BIC point to the model with 3 parameters, adjusted $R^2$ prefers the model with 4 parameters. It's important to note here that we do not only choose the model with smaller dimensions just because sparsity is always preferred. It should be noted that adjusted $R^2$ is choosing the model with 4 parameters because it is trying to maximize the variance explained by the model. We choose the model with lower complexity here because maximially explaining the variance within the data is not necessarily the primary goal (assuming that the primary goal is generalizable inference).\\\n",
        "The betas I chose in this 'true' model I created were: beta0 <- -1.67, beta1 <- 0.79, beta2 <- 1.9, beta3 <- 1.55. The best model chosen through bias-adjusted fit comparison produces coefficients that are very close to the 'true' values."
      ],
      "id": "e7f16d76"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "19d388b0"
      },
      "source": [
        "---\n",
        "## 2. Forward and backwards stepwise selection (3 points)"
      ],
      "id": "19d388b0"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "620dc14c"
      },
      "source": [
        "Using the same simulated data from question 1, use forward stepwise selection, and backwards stepwise selections to determine the best model. Again, for both model selection methods, plot the Mallows' Cp, Bayesian Information Criteria (BIC) and the adjusted coefficient of determination, $R^2$, for each model tested. Report the coefficients of the best model."
      ],
      "id": "620dc14c"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d734954e"
      },
      "source": [
        "a) Forward stepwise selection:"
      ],
      "id": "d734954e"
    },
    {
      "cell_type": "code",
      "execution_count": 85,
      "metadata": {
        "id": "8018ba62",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 86
        },
        "outputId": "7030db4a-60f9-4e85-a156-97cf867d39c8"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "3"
            ],
            "text/markdown": "3",
            "text/latex": "3",
            "text/plain": [
              "[1] 3"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "3"
            ],
            "text/markdown": "3",
            "text/latex": "3",
            "text/plain": [
              "[1] 3"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "4"
            ],
            "text/markdown": "4",
            "text/latex": "4",
            "text/plain": [
              "[1] 4"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<style>\n",
              ".dl-inline {width: auto; margin:0; padding: 0}\n",
              ".dl-inline>dt, .dl-inline>dd {float: none; width: auto; display: inline-block}\n",
              ".dl-inline>dt::after {content: \":\\0020\"; padding-right: .5ex}\n",
              ".dl-inline>dt:not(:first-of-type) {padding-left: .5ex}\n",
              "</style><dl class=dl-inline><dt>(Intercept)</dt><dd>-1.7786055643429</dd><dt>X</dt><dd>0.765488944992593</dd><dt>X2</dt><dd>1.99154367427184</dd><dt>X3</dt><dd>1.56719996937048</dd></dl>\n"
            ],
            "text/markdown": "(Intercept)\n:   -1.7786055643429X\n:   0.765488944992593X2\n:   1.99154367427184X3\n:   1.56719996937048\n\n",
            "text/latex": "\\begin{description*}\n\\item[(Intercept)] -1.7786055643429\n\\item[X] 0.765488944992593\n\\item[X2] 1.99154367427184\n\\item[X3] 1.56719996937048\n\\end{description*}\n",
            "text/plain": [
              "(Intercept)           X          X2          X3 \n",
              " -1.7786056   0.7654889   1.9915437   1.5672000 "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "regfit.fwd <- regsubsets(Y ~., data = data, method = \"forward\")\n",
        "#summary(regfit.fwd)\n",
        "reg.fwd.summary = summary(regfit.fwd)\n",
        "which.min(reg.fwd.summary$cp)\n",
        "which.min(reg.fwd.summary$bic)\n",
        "which.max(reg.fwd.summary$adjr2)\n",
        "\n",
        "coef(regfit.fwd, 3)"
      ],
      "id": "8018ba62"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "861cedcd",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# Note: if your model summary doesn't show up when viewing your notebook on github, copy and paste the output below.\n",
        "# Otherwise you can either delete this cell or leave it blank.\n",
        "\n",
        "```\n",
        "# paste here\n",
        "\n",
        "```"
      ],
      "id": "861cedcd"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "955e381d"
      },
      "source": [
        "b) Backwards stepwise selection:"
      ],
      "id": "955e381d"
    },
    {
      "cell_type": "code",
      "execution_count": 86,
      "metadata": {
        "id": "c9a50756",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 86
        },
        "outputId": "88227aae-de19-42da-a141-16f8fc49985b"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "6"
            ],
            "text/markdown": "6",
            "text/latex": "6",
            "text/plain": [
              "[1] 6"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "6"
            ],
            "text/markdown": "6",
            "text/latex": "6",
            "text/plain": [
              "[1] 6"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "8"
            ],
            "text/markdown": "8",
            "text/latex": "8",
            "text/plain": [
              "[1] 8"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<style>\n",
              ".dl-inline {width: auto; margin:0; padding: 0}\n",
              ".dl-inline>dt, .dl-inline>dd {float: none; width: auto; display: inline-block}\n",
              ".dl-inline>dt::after {content: \":\\0020\"; padding-right: .5ex}\n",
              ".dl-inline>dt:not(:first-of-type) {padding-left: .5ex}\n",
              "</style><dl class=dl-inline><dt>(Intercept)</dt><dd>-1.79456846310476</dd><dt>X</dt><dd>1.5061126431511</dd><dt>X2</dt><dd>1.97615390578811</dd><dt>X5</dt><dd>0.907296456592658</dd><dt>X7</dt><dd>-0.186214665706477</dd><dt>X9</dt><dd>0.0119835636003406</dd><dt>X10</dt><dd>0.000203221769327799</dd></dl>\n"
            ],
            "text/markdown": "(Intercept)\n:   -1.79456846310476X\n:   1.5061126431511X2\n:   1.97615390578811X5\n:   0.907296456592658X7\n:   -0.186214665706477X9\n:   0.0119835636003406X10\n:   0.000203221769327799\n\n",
            "text/latex": "\\begin{description*}\n\\item[(Intercept)] -1.79456846310476\n\\item[X] 1.5061126431511\n\\item[X2] 1.97615390578811\n\\item[X5] 0.907296456592658\n\\item[X7] -0.186214665706477\n\\item[X9] 0.0119835636003406\n\\item[X10] 0.000203221769327799\n\\end{description*}\n",
            "text/plain": [
              "  (Intercept)             X            X2            X5            X7 \n",
              "-1.7945684631  1.5061126432  1.9761539058  0.9072964566 -0.1862146657 \n",
              "           X9           X10 \n",
              " 0.0119835636  0.0002032218 "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "regfit.bwd <- regsubsets(Y ~., data = data, method = \"backward\")\n",
        "#summary(regfit.bwd)\n",
        "reg.bwd.summary = summary(regfit.bwd)\n",
        "which.min(reg.bwd.summary$cp)\n",
        "which.min(reg.bwd.summary$bic)\n",
        "which.max(reg.bwd.summary$adjr2)\n",
        "\n",
        "coef(regfit.bwd, 6)"
      ],
      "id": "c9a50756"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "30746619",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# Note: if your model summary doesn't show up when viewing your notebook on github, copy and paste the output below.\n",
        "# Otherwise you can either delete this cell or leave it blank.\n",
        "\n",
        "```\n",
        "# paste here\n",
        "\n",
        "```"
      ],
      "id": "30746619"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "25dbc027"
      },
      "source": [
        "c) Compare your results from parts a and b with those of question 1."
      ],
      "id": "25dbc027"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "34ac850e"
      },
      "source": [
        "> Forward stepwise selection produced the same results as those of question 1, but backward stepwise selection chose the model with 6 parameters. This is due to how the predictors are created - remember that X2 is $x^2$, X3 is $x^3$ and so on. This means that the predictor variables have high collinearity. In a backward stepwise selection this means that the algorithm stops removing terms earlier, resulting in a model with more predictors compared to forward stepwise or best subset selection.\n",
        ">"
      ],
      "id": "34ac850e"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "630f4909"
      },
      "source": [
        "---\n",
        "## 3. Training and test error (3 points)"
      ],
      "id": "630f4909"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cdea3f49"
      },
      "source": [
        "This question will explore the relationship between training and test error and the number of features included in a model. We will again use a simulated dataset.\n",
        "\n",
        "a) Simulate dataset with features $p = 20$ and observations $n = 1,000$. Generate data for a response variable $Y$ according to the model\n",
        "\n",
        "$Y = X\\beta + \\epsilon$\n",
        "\n",
        "where $\\beta$ is random with some elements that are exactly zero.\n",
        "\n",
        "Split your simulated dataset into a training set containing $n=100$ observations and a test set containing $n=900$ observations."
      ],
      "id": "cdea3f49"
    },
    {
      "cell_type": "code",
      "execution_count": 97,
      "metadata": {
        "id": "3f113879",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 324
        },
        "outputId": "3e934070-ca43-4c03-9931-3708e4f57e91"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<style>\n",
              ".list-inline {list-style: none; margin:0; padding: 0}\n",
              ".list-inline>li {display: inline-block}\n",
              ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
              "</style>\n",
              "<ol class=list-inline><li>-1.1144403712824</li><li>0.52614872250706</li><li>0</li><li>1.03203241433948</li><li>-1.46307313907892</li><li>0.51514502428472</li><li>0</li><li>0</li><li>0</li><li>0</li><li>0</li><li>-0.527756507508457</li><li>0</li><li>-1.89949801936746</li><li>0</li><li>-0.712022013962269</li><li>0</li><li>0</li><li>-1.55381905566901</li><li>-0.679967186413705</li></ol>\n"
            ],
            "text/markdown": "1. -1.1144403712824\n2. 0.52614872250706\n3. 0\n4. 1.03203241433948\n5. -1.46307313907892\n6. 0.51514502428472\n7. 0\n8. 0\n9. 0\n10. 0\n11. 0\n12. -0.527756507508457\n13. 0\n14. -1.89949801936746\n15. 0\n16. -0.712022013962269\n17. 0\n18. 0\n19. -1.55381905566901\n20. -0.679967186413705\n\n\n",
            "text/latex": "\\begin{enumerate*}\n\\item -1.1144403712824\n\\item 0.52614872250706\n\\item 0\n\\item 1.03203241433948\n\\item -1.46307313907892\n\\item 0.51514502428472\n\\item 0\n\\item 0\n\\item 0\n\\item 0\n\\item 0\n\\item -0.527756507508457\n\\item 0\n\\item -1.89949801936746\n\\item 0\n\\item -0.712022013962269\n\\item 0\n\\item 0\n\\item -1.55381905566901\n\\item -0.679967186413705\n\\end{enumerate*}\n",
            "text/plain": [
              " [1] -1.1144404  0.5261487  0.0000000  1.0320324 -1.4630731  0.5151450\n",
              " [7]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 -0.5277565\n",
              "[13]  0.0000000 -1.8994980  0.0000000 -0.7120220  0.0000000  0.0000000\n",
              "[19] -1.5538191 -0.6799672"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<table class=\"dataframe\">\n",
              "<caption>A data.frame: 6 × 21</caption>\n",
              "<thead>\n",
              "\t<tr><th></th><th scope=col>x1</th><th scope=col>x2</th><th scope=col>x3</th><th scope=col>x4</th><th scope=col>x5</th><th scope=col>x6</th><th scope=col>x7</th><th scope=col>x8</th><th scope=col>x9</th><th scope=col>x10</th><th scope=col>⋯</th><th scope=col>x12</th><th scope=col>x13</th><th scope=col>x14</th><th scope=col>x15</th><th scope=col>x16</th><th scope=col>x17</th><th scope=col>x18</th><th scope=col>x19</th><th scope=col>x20</th><th scope=col>y</th></tr>\n",
              "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>⋯</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
              "</thead>\n",
              "<tbody>\n",
              "\t<tr><th scope=row>1</th><td> 1.6121462</td><td>-0.06232354</td><td>-0.8529892</td><td>-0.6289533</td><td>-1.25375539</td><td> 0.2784459</td><td>-0.85727115</td><td>-1.0380393</td><td>-0.07431788</td><td> 1.1696039</td><td>⋯</td><td>-0.284183231</td><td>-0.04692446</td><td> 0.75231624</td><td>-0.2388338</td><td> 0.5859516</td><td>-0.5788594</td><td> 0.3542662</td><td>-0.1278698</td><td> 0.12284860</td><td>-2.3572036</td></tr>\n",
              "\t<tr><th scope=row>2</th><td>-0.5912932</td><td>-0.85956918</td><td> 0.4191234</td><td>-0.4400123</td><td>-1.10531031</td><td>-0.6212755</td><td>-1.04795561</td><td> 0.5394848</td><td> 0.31067039</td><td> 0.3045565</td><td>⋯</td><td> 0.005294635</td><td>-0.79788145</td><td>-0.62964470</td><td> 0.8261991</td><td>-1.1866012</td><td> 2.7034960</td><td>-0.1840551</td><td> 0.4347243</td><td>-0.99645577</td><td> 3.2827437</td></tr>\n",
              "\t<tr><th scope=row>3</th><td>-0.7221665</td><td> 0.77325349</td><td> 0.6327792</td><td> 0.3279123</td><td> 0.58771375</td><td> 0.9466239</td><td> 1.35191122</td><td>-0.8663276</td><td> 0.91656613</td><td> 1.8614002</td><td>⋯</td><td> 1.585804018</td><td> 0.46803560</td><td>-1.43875859</td><td>-0.9038115</td><td> 2.3467646</td><td>-0.8134583</td><td>-0.9581906</td><td> 1.3986076</td><td>-0.08577181</td><td>-0.1969125</td></tr>\n",
              "\t<tr><th scope=row>4</th><td>-2.3762106</td><td>-0.73012714</td><td> 1.2426104</td><td>-0.9584096</td><td> 0.08107354</td><td> 0.7792539</td><td>-0.64568154</td><td>-0.4262454</td><td> 1.11165865</td><td>-0.8757133</td><td>⋯</td><td> 0.594952882</td><td>-0.34141192</td><td> 1.17840045</td><td> 0.6768907</td><td> 1.8649899</td><td>-0.7013750</td><td>-1.8685490</td><td>-0.7044939</td><td>-0.12502707</td><td>-1.0820421</td></tr>\n",
              "\t<tr><th scope=row>5</th><td>-0.6461708</td><td>-1.19135815</td><td>-0.7292524</td><td> 0.3907425</td><td> 0.45241382</td><td>-2.1570203</td><td>-0.46927108</td><td>-0.7344011</td><td>-1.19135498</td><td> 0.5722186</td><td>⋯</td><td> 1.240175029</td><td>-1.78038509</td><td> 0.70213007</td><td> 0.2696388</td><td> 0.1367579</td><td> 1.0203101</td><td> 1.1488594</td><td> 0.1038445</td><td> 0.51248933</td><td>-4.4028769</td></tr>\n",
              "\t<tr><th scope=row>6</th><td>-0.7254025</td><td>-0.01360758</td><td>-1.1735859</td><td> 0.9665259</td><td> 0.77204310</td><td>-1.6154565</td><td> 0.09176664</td><td> 1.3567897</td><td>-0.42705313</td><td> 0.0141786</td><td>⋯</td><td>-0.717254831</td><td> 0.21792386</td><td> 0.01006601</td><td> 1.1869282</td><td> 1.2060878</td><td>-0.7245896</td><td>-1.5859460</td><td>-0.2024745</td><td>-0.90069934</td><td> 1.0594188</td></tr>\n",
              "</tbody>\n",
              "</table>\n"
            ],
            "text/markdown": "\nA data.frame: 6 × 21\n\n| <!--/--> | x1 &lt;dbl&gt; | x2 &lt;dbl&gt; | x3 &lt;dbl&gt; | x4 &lt;dbl&gt; | x5 &lt;dbl&gt; | x6 &lt;dbl&gt; | x7 &lt;dbl&gt; | x8 &lt;dbl&gt; | x9 &lt;dbl&gt; | x10 &lt;dbl&gt; | ⋯ ⋯ | x12 &lt;dbl&gt; | x13 &lt;dbl&gt; | x14 &lt;dbl&gt; | x15 &lt;dbl&gt; | x16 &lt;dbl&gt; | x17 &lt;dbl&gt; | x18 &lt;dbl&gt; | x19 &lt;dbl&gt; | x20 &lt;dbl&gt; | y &lt;dbl&gt; |\n|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n| 1 |  1.6121462 | -0.06232354 | -0.8529892 | -0.6289533 | -1.25375539 |  0.2784459 | -0.85727115 | -1.0380393 | -0.07431788 |  1.1696039 | ⋯ | -0.284183231 | -0.04692446 |  0.75231624 | -0.2388338 |  0.5859516 | -0.5788594 |  0.3542662 | -0.1278698 |  0.12284860 | -2.3572036 |\n| 2 | -0.5912932 | -0.85956918 |  0.4191234 | -0.4400123 | -1.10531031 | -0.6212755 | -1.04795561 |  0.5394848 |  0.31067039 |  0.3045565 | ⋯ |  0.005294635 | -0.79788145 | -0.62964470 |  0.8261991 | -1.1866012 |  2.7034960 | -0.1840551 |  0.4347243 | -0.99645577 |  3.2827437 |\n| 3 | -0.7221665 |  0.77325349 |  0.6327792 |  0.3279123 |  0.58771375 |  0.9466239 |  1.35191122 | -0.8663276 |  0.91656613 |  1.8614002 | ⋯ |  1.585804018 |  0.46803560 | -1.43875859 | -0.9038115 |  2.3467646 | -0.8134583 | -0.9581906 |  1.3986076 | -0.08577181 | -0.1969125 |\n| 4 | -2.3762106 | -0.73012714 |  1.2426104 | -0.9584096 |  0.08107354 |  0.7792539 | -0.64568154 | -0.4262454 |  1.11165865 | -0.8757133 | ⋯ |  0.594952882 | -0.34141192 |  1.17840045 |  0.6768907 |  1.8649899 | -0.7013750 | -1.8685490 | -0.7044939 | -0.12502707 | -1.0820421 |\n| 5 | -0.6461708 | -1.19135815 | -0.7292524 |  0.3907425 |  0.45241382 | -2.1570203 | -0.46927108 | -0.7344011 | -1.19135498 |  0.5722186 | ⋯ |  1.240175029 | -1.78038509 |  0.70213007 |  0.2696388 |  0.1367579 |  1.0203101 |  1.1488594 |  0.1038445 |  0.51248933 | -4.4028769 |\n| 6 | -0.7254025 | -0.01360758 | -1.1735859 |  0.9665259 |  0.77204310 | -1.6154565 |  0.09176664 |  1.3567897 | -0.42705313 |  0.0141786 | ⋯ | -0.717254831 |  0.21792386 |  0.01006601 |  1.1869282 |  1.2060878 | -0.7245896 | -1.5859460 | -0.2024745 | -0.90069934 |  1.0594188 |\n\n",
            "text/latex": "A data.frame: 6 × 21\n\\begin{tabular}{r|lllllllllllllllllllll}\n  & x1 & x2 & x3 & x4 & x5 & x6 & x7 & x8 & x9 & x10 & ⋯ & x12 & x13 & x14 & x15 & x16 & x17 & x18 & x19 & x20 & y\\\\\n  & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & ⋯ & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n\\hline\n\t1 &  1.6121462 & -0.06232354 & -0.8529892 & -0.6289533 & -1.25375539 &  0.2784459 & -0.85727115 & -1.0380393 & -0.07431788 &  1.1696039 & ⋯ & -0.284183231 & -0.04692446 &  0.75231624 & -0.2388338 &  0.5859516 & -0.5788594 &  0.3542662 & -0.1278698 &  0.12284860 & -2.3572036\\\\\n\t2 & -0.5912932 & -0.85956918 &  0.4191234 & -0.4400123 & -1.10531031 & -0.6212755 & -1.04795561 &  0.5394848 &  0.31067039 &  0.3045565 & ⋯ &  0.005294635 & -0.79788145 & -0.62964470 &  0.8261991 & -1.1866012 &  2.7034960 & -0.1840551 &  0.4347243 & -0.99645577 &  3.2827437\\\\\n\t3 & -0.7221665 &  0.77325349 &  0.6327792 &  0.3279123 &  0.58771375 &  0.9466239 &  1.35191122 & -0.8663276 &  0.91656613 &  1.8614002 & ⋯ &  1.585804018 &  0.46803560 & -1.43875859 & -0.9038115 &  2.3467646 & -0.8134583 & -0.9581906 &  1.3986076 & -0.08577181 & -0.1969125\\\\\n\t4 & -2.3762106 & -0.73012714 &  1.2426104 & -0.9584096 &  0.08107354 &  0.7792539 & -0.64568154 & -0.4262454 &  1.11165865 & -0.8757133 & ⋯ &  0.594952882 & -0.34141192 &  1.17840045 &  0.6768907 &  1.8649899 & -0.7013750 & -1.8685490 & -0.7044939 & -0.12502707 & -1.0820421\\\\\n\t5 & -0.6461708 & -1.19135815 & -0.7292524 &  0.3907425 &  0.45241382 & -2.1570203 & -0.46927108 & -0.7344011 & -1.19135498 &  0.5722186 & ⋯ &  1.240175029 & -1.78038509 &  0.70213007 &  0.2696388 &  0.1367579 &  1.0203101 &  1.1488594 &  0.1038445 &  0.51248933 & -4.4028769\\\\\n\t6 & -0.7254025 & -0.01360758 & -1.1735859 &  0.9665259 &  0.77204310 & -1.6154565 &  0.09176664 &  1.3567897 & -0.42705313 &  0.0141786 & ⋯ & -0.717254831 &  0.21792386 &  0.01006601 &  1.1869282 &  1.2060878 & -0.7245896 & -1.5859460 & -0.2024745 & -0.90069934 &  1.0594188\\\\\n\\end{tabular}\n",
            "text/plain": [
              "  x1         x2          x3         x4         x5          x6        \n",
              "1  1.6121462 -0.06232354 -0.8529892 -0.6289533 -1.25375539  0.2784459\n",
              "2 -0.5912932 -0.85956918  0.4191234 -0.4400123 -1.10531031 -0.6212755\n",
              "3 -0.7221665  0.77325349  0.6327792  0.3279123  0.58771375  0.9466239\n",
              "4 -2.3762106 -0.73012714  1.2426104 -0.9584096  0.08107354  0.7792539\n",
              "5 -0.6461708 -1.19135815 -0.7292524  0.3907425  0.45241382 -2.1570203\n",
              "6 -0.7254025 -0.01360758 -1.1735859  0.9665259  0.77204310 -1.6154565\n",
              "  x7          x8         x9          x10        ⋯ x12          x13        \n",
              "1 -0.85727115 -1.0380393 -0.07431788  1.1696039 ⋯ -0.284183231 -0.04692446\n",
              "2 -1.04795561  0.5394848  0.31067039  0.3045565 ⋯  0.005294635 -0.79788145\n",
              "3  1.35191122 -0.8663276  0.91656613  1.8614002 ⋯  1.585804018  0.46803560\n",
              "4 -0.64568154 -0.4262454  1.11165865 -0.8757133 ⋯  0.594952882 -0.34141192\n",
              "5 -0.46927108 -0.7344011 -1.19135498  0.5722186 ⋯  1.240175029 -1.78038509\n",
              "6  0.09176664  1.3567897 -0.42705313  0.0141786 ⋯ -0.717254831  0.21792386\n",
              "  x14         x15        x16        x17        x18        x19       \n",
              "1  0.75231624 -0.2388338  0.5859516 -0.5788594  0.3542662 -0.1278698\n",
              "2 -0.62964470  0.8261991 -1.1866012  2.7034960 -0.1840551  0.4347243\n",
              "3 -1.43875859 -0.9038115  2.3467646 -0.8134583 -0.9581906  1.3986076\n",
              "4  1.17840045  0.6768907  1.8649899 -0.7013750 -1.8685490 -0.7044939\n",
              "5  0.70213007  0.2696388  0.1367579  1.0203101  1.1488594  0.1038445\n",
              "6  0.01006601  1.1869282  1.2060878 -0.7245896 -1.5859460 -0.2024745\n",
              "  x20         y         \n",
              "1  0.12284860 -2.3572036\n",
              "2 -0.99645577  3.2827437\n",
              "3 -0.08577181 -0.1969125\n",
              "4 -0.12502707 -1.0820421\n",
              "5  0.51248933 -4.4028769\n",
              "6 -0.90069934  1.0594188"
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "set.seed(1997)\n",
        "# parameters\n",
        "n <- 1000\n",
        "p <- 20\n",
        "\n",
        "beta <- runif(p, min = -2, max = 2) # vector of length 20.\n",
        "beta[sample(1:p, 10)] <- 0  # setting half of them to 0\n",
        "beta\n",
        "\n",
        "dat <- matrix(rnorm(n * p), nrow = n, ncol = p)\n",
        "colnames(dat) <- paste0(\"x\", 1:20)\n",
        "\n",
        "e <- rnorm(n, mean = 0, sd = 0.5)\n",
        "\n",
        "y <- colSums((t(dat) * beta)) + e\n",
        "\n",
        "dat <- data.frame(dat)\n",
        "dat$y <- y\n",
        "\n",
        "head(dat)\n",
        "\n",
        "train <- dat[1:100, ]\n",
        "test <- dat[101:1000, ]"
      ],
      "id": "3f113879"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5886639c"
      },
      "source": [
        "b) Perform best subset selection on the training set, and plot the associated training and test set MSE for the best model of each size."
      ],
      "id": "5886639c"
    },
    {
      "cell_type": "code",
      "execution_count": 95,
      "metadata": {
        "id": "a1dbdf04",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 857
        },
        "outputId": "6c0c69d2-52cd-4b3a-e6c1-306d9f69cc1b"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "plot without title"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3deWBU5dXH8TNZ2EHZBako7isq\nLiguVaxWS1CKuCHUBWgtFlutYq2vqfpWqUvFpVbcabVWrBsuVEWr4itVUXBFVAQEQZxsJJBA\nQnLemSSQQJaZ38y59z53+H3/gMnMM/dMMvOBSe6djChjLO0k6BvAWCZESIwZREiMGURIjBlE\nSIwZREiMGURIjBlESIwZZA+pfIP5JptWU1TmwxQt3ejDkMqich+maIkfQyqKKn2YUl3qwxAt\nK6oGVttDWrvefJNNq46u8WGKlvgBaUN0nQ9TtNCPIeuifvwzutGXfxTWRC0hLb9u9JlXforM\nJyQwQkILIaSa8XeuW//IKOS/UkICIyS0EEIqyVuoWpS3CNgiIYEREloIIenlt5VW/GPc5q9O\nVWXCStclXpN266MlPkypLF7vw5B10TIfplQW+jGkNOrLvV/sw5DK4mgS937SkAon5uWN/Wrz\nhyVRxtimipKFVHXJnSXrnhi9ef368oSVlCZek3brokU+TCkvWuvDkLLoGh+mlBf4MaQkWubD\nlLX+3PnRdQnXVCQL6YPh8aUXzEzw/1bj+D0SGL9HQgvh90jv58Xv5LGE5GGEhBZCSOvG3lm2\n4cmRK4EtEhIYIaGFEJIuvXb0WVd8hMwnJDBCQgsjJDxCAiMkNEIyi5DgCAmMkOwiJDBCsouQ\nwAgJjZDMIiQ4QgIjJLsICYyQ7CIkMEJCIySzCAmOkMAIyS5CAiMkuwgJjJDQCMksQoIjJDBC\nsouQwAjJLkICIyQ0QjKLkOAICYyQ7CIkMEKyi5DACAmNkMwiJDhCAiMkuwgJjJDsIiQwQkLb\nJiCVvjbXB0mEBEdIYMFCur2TyI7PGm+0aYQER0hggUJ6UuJ1+Nh2q00jJDhCAgsU0uG1kGSc\n7VabRkhwhAQWKKTedZCG2m61aYQER0hggUI6qA7SWNutNo2Q4AgJLFBI99VBett2q00jJDhC\nAgv2p3ZXtIk5Otl4o00jJDhCAgt4P9LSvz/aq9031lvdOkKCIySw4I9suFsuMN/qVhESHCGB\nBQ9p497ZXu9IIiQ4QgILHpI+KcPMN7tlhARHSGAOQNIh8qr5dreIkOAICcwFSHMjh9aYb7hx\nhARHSGAuQNLhMsN8w40jJDhCAnMC0uc5Azz9ihISHCGBOQFJx8td5ltuFCHBERKYG5BWduzp\n5UOdkOAICcwNSPp7ucZ80w0REhwhgTkCqbR3x5Xm294cIcEREpgjkPR2+bn5tjdHSHCEBOYK\npMrdsj813/imCAmOkMBcgaSPyQjzjW+KkOAICcwZSDWD5S3zrddHSHCEBOYMJH1dBnt1oBAh\nwRESmDuQ9GR52nzzdRESHCGBOQTpo+w9q8y3XxshwRESmEOQ9GcyzXz7tRESHCGBuQRpRYc+\na80HxCMkOEICcwmSXiHXmw+IR0hwhATmFKTi7p2+M5+ghJRChATmFCS9RX5lPkEJKYUICcwt\nSBsG5H5hPoKQUoiQwNyCpH+TM8xHEFIKERKYY5CqD4548JvACQmOkMAcg6QvyTHmMwgJj5DA\nXIOkP5IXzIcQEhwhgTkHaUHW3uYHChESHCGBOQdJz5EHrYcQEhwhgbkHaUnbHa0fKoQER0hg\n7kHSX8sU4yGEBEdIYA5CKuq2fYHtEEKCIyQwByHpjXKp7RBCgiMkMBchVfygzWLTIYQER0hg\nLkLSB+Rc0yGEBEdIYE5C2rhfZJ7lEEKCIyQwJyHpc3Kc5RBCgiMkMDch6fHysuEQQoIjJDBH\nIb0bGYjcrgQREhwhgTkKSUfJI3ZDCAmOkMBchbS4zc4tXJJChARHSGCuQtKJcqvZEEKCIyQw\nZyF936Wr2Z1JSHCEBOYsJL1WJlsNISQ4QgJzF9LaPu2WGQ0hJDhCAnMXkt4j5xsNISQ4QgJz\nGNLGfbLm2wwhJDhCAnMYkj4lp9gMISQ4QgJzGZIOkdkmQwgJjpDAnIY0N3KoybthEhIcIYE5\nDUlPk8cthhASHCGBuQ3p85xdLL7UhARHSGBuQ9IJcofBEEKCIyQwxyGt7tzDwAAhwRESmOOQ\n9Gq5Ov0hhARHSGCuQyrr3X552kMICY6QwFyHpHfIhLSHEBIcIYE5D6ly9+xP0x1CSHCEBOY8\nJH1cTkt3CCHBERKY+5BqjpA5aQ4hJDhCAnMfkr4hh6d5oBAhwRESWAgg6SnyVHpDCAmOkMDC\nAOmznD0q0xpCSHCEBBYGSHq+3JPWEEKCIySwUEBa0aFXaTpDCAmOkMA8hVS9MWFl5YnXbNx4\nheQns6ylKqMl6Vw92Yo3+DCkIrrWhykbC/0YUhat8GHKhmIfhmwsiVYmXNNADYS0tiRhhUWJ\n15SULOvWcVEy61qoOFqYxrWTrqDYhyFF0aS+YukW9WNIoS+fS3GBD0NKCqKJ7/2Gp1UBPbVT\nvVUmpjGET+3g+NQOLBTfI8UePrvmfpH6EEKCIySwkEDSR2RU6kMICY6QwMICqWZQ5O2UhxAS\nHCGBhQWSvibHpDyEkOAICSw0kPRH8lyqQwgJjpDAwgNpQdZeVSkOISQ4QgILDyQ9Vx5IcQgh\nwRESWIggLW3bN8XHECHBERJYiCDppXJjakMICY6QwMIEqajbdtGUhhASHCGBhQmSTpFfpzSE\nkOAICSxUkCp2avNVKkMICY6QwEIFSR+Uc1IZQkhwhAQWLkjVB0bmpTCEkOAICSxckPQF+WEK\nQwgJjpDAQgZJh8pL+BBCgiMksLBBejdyAHKD6yIkOEICCxskPUP+Bg8hJDhCAgsdpK/b9CtH\nhxASHCGBhQ6S/kpuQa9CSHCEBBY+SNEuXdF7mZDgCAksfJD0erkCvAYhwRESWAghlfdrtwy7\nBiHBERJYCCHpNPkZdgVCgiMksDBC2rhP1gfQFQgJjpDAwghJn5aTofWEBEdIYKGEpEfJbGQ5\nIcERElg4If03ciByswkJjpDAwglJR8hjwGpCgiMksJBCWpS7C3BFQoIjJLCQQtKfy+3JLyYk\nOEICCyuk1Z17JP/1ISQ4QgILKyS9Rn6f9FpCgiMksNBCKtuh/fJk1xISHCGBhRaS3iXjk11K\nSHCEBBZeSJW7Z3+S5FJCgiMksPBC0hkyPMmVhARHSGAhhlRzhLyW3EpCgiMksBBD0jfksJqk\nFhISHCGBhRmSDpMnk1pHSHCEBBZqSAtz9qhMZh0hwRESWKgh6QVydzLLCAmOkMDCDenbDr1K\nk1hGSHCEBBZuSHql/CGJVYQER0hgIYdU0qPTqsSrCAmOkMBCDklvk18mXkRIcIQEFnZIG3bN\n+SzhIkKCIySwsEPSR+X0hGsICY6QwEIPqWaQ/F+iNYQER0hgoYek/5GjEy0hJDhCAgs/JD1J\nZiZYQUhwhASWAZA+zNqrqvUVhARHSGAZAEnHyP2tLyAkOEICywRIy9v3XdvqAkKCIySwTICk\nl8kfW72ckOAICSwjIBV167y6tcsJCY6QwDICkv5JLmntYkKCIySwzIBUsVObr1q5mJDgCAks\nMyDpw3J2K5cSEhwhgWUIpOqDIu+1cikhoRESWIZA0hfl2JYvJCQ4QgLLFEg6VGa1eBkhwRES\nWMZAmp+1f4ufCSHBERJYxkDSs2R6SxcREhwhgWUOpK/b7ljewkWEBEdIYJkDSSfJzS1cQkhw\nhASWQZCi221f0PwlhARHSGAZBEn/V37b/AWEBEdIYJkEqbxfu6XNXkBIcIQElkmQ9F4Z2+z5\nhARHSGAZBWnjvlkfNHc+IcERElhGQdJn5cfNnU1IcIQEllmQ9Dh5pZlzCQmOkMAyDNI7kYHN\nfD6EBEdIYBkGSX8q/2h6JiHBERJYpkFalLtz0y0SEhwhgWUaJL1IpjY5j5DgCAks4yCt7tyj\nyReOkOAICSzjIGm+XLX1WYQER0hgmQepbIf232x1FiHBERJY5kHSv8iFW51DSHCEBJaBkKr2\nzv5ky3MICY6QwDIQkv5L8rY8g5DgCAksEyHpkfLaFh8TEhwhgWUkpDflsJrGHxMSHCGBZSQk\nzZMnGn9ISHCEBJaZkBbmDGh8FxESHCGBZSYkvVD+0ugjQoIjJLAMhfRtx56N7BASHCGBZSgk\nvUryGz4gJDhCAstUSKW9Oq3a/AEhwRESWKZC0qly0ebThARHSGAZC2nDrtmfbTpNSHCEBJax\nkPQf8tNNJwkJjpDAjCG9MG7Exe8i872DVHOIvFV/kpDgCAnMFtLsse+tfmY8ck97B0lfl8H1\nBwoREhwhgdlCGv8qOt9DSPpjebbuBCHBERKYKaSCvFd/dfplC5H5XkL6KGvPqtoThARHSGCm\nkBbl/W556b1nbb7hpUUJKyxMvCblzpKpdUOiBR5O2ZwvQwqjXn7FNufPFyyjPpfEaxpEJ4S0\nIPYPwNmzNystSFg08ZLU+7hd729qh0Q9HbMpf4b4MoWfCzwkiSnFyUKK5n0Z+3PiE62v2iIv\nn9qpXi7/G/+LT+3g+NQOzPSpXfXY52L39BlvAlv0FlJx987fKSGlECGB2f7U7onR86N3jK0A\ntugtJL1ZJikhpRAhgdlCqp4+ZsSVW/9WuVbzGFJF/9wvCSmFCAkscw8Rqmu6nElIKURIYJkO\nqfqgyFxCwiMksEyHpLPkWELCIySwjIekJ8iLhARHSGCZD2lB1v6VhIRGSGCZD0nPlgcJCY2Q\nwLYBSEvadDz4+LuqPJ9DSGiEZJcPkL5oI7GG1SRemWaEBEZIdvkA6QSp7VHPBxESGCHZ5T2k\n6tw6SOO8HkRIaIRkFyGBERIaIRk1lE/t8AgJbBuA9HnnuKNj+cMGJEIC2wYg6bKfH7h3/NhV\nryMkMEKyyw9IWh0tPjAy3/MxhARGSHb5BGnNszLM8zGEBEZIdvkFSY+UN7weQ0hghGSXb5De\nlKO9HkNIYIRkl2+Q9ET5t8djCAmMkOzyD9K8yCCPfwJOSGCEZJd/kHSkIL9wL4UICYyQ7PIR\n0uc5e3j7WgpCAiMku3yEpOfJg56OISQwQrLLT0hL2+7k6TRCAiMku/yEpL+SqV6OISQwQrLL\nV0jfd+5Z6uEYQgIjJLt8haRXyfUejiEkMEKyy19IJd22K/BuDCGBEZJd/kLSG+UK78YQEhgh\n2eUzpPId2y33bAwhgRGSXT5D0rvkIs/GEBIYIdnlN6TKXXO/8moMIYERkl1+Q9K/yblejSEk\nMEKyy3dI1QOzFng0hpDACMku3yHp03KqR2MICYyQ7PIfkg6Wt70ZQ0hghGRXAJBel2O8GUNI\nYIRkVwCQdKi84skYQgIjJLuCgPRu5BBPXnROSGCEZFcQkPQ0edKLMYQERkh2BQLpk+w9vXjR\nOSGBEZJdgUDSMfKwB2MICYyQ7AoG0pI2/T2YS0hghGRXMJD0l3Kn/RhCAiMkuwKCtLKjBy86\nJyQwQrIrIEh6pdxgPoaQwAjJrqAgFXfb3vyxQkhghGRXUJD0f+V31mMICYyQ7AoM0tod2q8w\nHkNIYIRkV2CQ9A652HgMIYERkl3BQdowIHex7RhCAiMku4KDpA/Jz2zHEBIYIdkVIKSN+2R/\nYjqGkMAIya4AIem/ZITpGEICIyS7goRUc7jMtRxDSGCEZFeQkPQlOcFyDCGBEZJdgULS4+VV\nwzGEBEZIdgUL6Z3IoYYvOickMEKyK1hImifP2I0hJDBCsitgSB9n7Yt8/q1HSGCEZFfAkHS0\n/N1sDCGBEZJdQUP6us3OZncmIYERkl1BQ9JfyN1WYwgJjJDsChzSyg59rB6ZhARGSHYFDkkv\nlylGYwgJjJDsCh5SUdfti2zGEBIYIdkVPCS9Tq62GUNIYIRklwOQynp3/M5kDCGBEZJdDkDS\n22SSyRhCAiMku1yAtGGXNl9bjCEkMEKyywVIer+cbzGGkMAIyS4nIG3cO/szgzGEBEZIdjkB\nSR+X0w3GEBIYIdnlBqSagyP/TX8MIYERkl1uQNJZclL6YwgJjJDscgSS/lBeS3sMIYERkl2u\nQHpLDkv7ReeEBEZIdrkCSX8iz6U7hpDACMkuZyB9lLV/ui86JyQwQrLLGUh6ljya5hhCAiMk\nu9yB9EXuLmnerYQERkh2uQNJx8s96Y0hJDBCssshSN+275veY5SQwAjJLocg6aVyc1pjCAmM\nkOxyCVK0S/ek1rUUIYERkl0uQdJ8uSadMYQERkh2OQWprFen1WmMISQwQrLLKUh6i/wmjTGE\nBEZIdrkFqeIH7ZalPoaQwAjJLrcg6b0yLvUxhARGSHY5BmnjXmm86JyQwAjJLscg6WNyZspj\nCAmMkOxyDVLNQZH3Ux1DSGCEZJdrkPR5OTnVMYQERkh2OQdJj5X/pDiGkMAIKcnWlSasqCTx\nmrRbEy1MfvFLMjjFMYVrUrwiUkm02IcppQV+DCmK+nLvA3d+6hVGE9/7ZalC2liVsNJ1idek\n3YZoCbD6JJmZ2pjiDaldD6o8WubDlKpCP4aURct9mLKh2IchVSXRxPd+wzOWbeGpnX6YdUBq\nLzrnUzswPrWzy0FIOkr+mdIYQgIjJLtchLQoZ/eqVMYQEhgh2eUiJL1Q7ktlDCGBEZJdTkJa\n0b5veQpjCAmMkOxyEpL+Wv6cwhhCAiMku9yEFO3cI4UXnRMSGCHZ5SYk/R/5Az6GkMAIyS5H\nIZV0T+FF54QERkh2OQpJb5LfwmMICYyQ7HIVUkW/dt+g1yEkMEKyy1VI+leZgF6FkMAIyS5n\nIVXulv05eBVCAiMku5yFpI/I2eA1CAmMkOxyF1L1gZH52DUICYyQ7HIXkj4rw7ArEBIYIdnl\nMCQ9Ul6H1hMSGCHZ5TKkN+UoaD0hgRGSXS5D0hPl38hyQgIjJLuchjQvMqgGWE5IYIRkl9OQ\ndKQ8AawmJDBCssttSJ/n7AG86JyQwAjJLrch6XnyQPKLCQmMkOxyHNLStjslfwMJCYyQ7HIc\nkv5Kpia9lpDACMku1yF937lnabJrCQmMkOxyHZJeJdcnu5SQwAjJLuchlXTbriDZpYSERUh2\nOQ9Jb5QrklxJSGCEZJf7kMp3bLc8uZWEBEZIdrkPSe+Si5JbSEhghGRXCCBV7pr7VVILCQmM\nkOwKAST9m5yb1DpCAiMku8IAqXpg1oJk1hESGCHZFQZI+rScmswyQgIjJLtCAUkHy9tJrCIk\nMEKyKxyQXpdjklhFSGCEZFc4IOlQeTnxIkICIyS7QgLpvcghiV90TkhghGRXSCDpafJkwjWE\nBEZIdoUF0ufZeyZ80TkhgRGSXWGBpGPk4URLCAmMkOwKDaQlbfsnuqmEBEZIdoUGkk6UOxOs\nICQwQrIrPJBWdkz0onNCAiMku8IDSa+UG1pfQEhghGRXiCAVd9u+9YcXIYERkl0hgqR/lN+1\nejkhgRGSXWGCtHaH9itau5yQwAjJrjBB0jvk4tYuJiQwQrIrVJA2DMhd3MrFhARGSHaFCpI+\nJGNbuZSQwAjJrnBB2rhP9ictX0pIYIRkV7gg6b9kRMsXEhIYIdkVMkg1h8vcFi8kJDBCsitk\nkPQlOaHFywgJjJDsChskPV5ebekiQgIjJLtCB+mdyKEtveickMAIya7QQdLh8kwLlxASGCHZ\nFT5IH2ft28KXjJDACEm1uOFLsOC2NOaHD5KOlr83fwEhgRFS7FT85aJlkxfF/rwznf+nQgjp\n6zY7N3//ExIYIdVDWiWzdNuDpL+Qu5s9n5DACGnbhrSyQ59mH8yEBEZI2zYkvVymNHc2IYER\n0jYOqajr9kXNnE1IYIS0jUPS6+TqZs4lJDBC2tYhlfXu+F3TcwkJjJBipy6bO3fuC3Jb7M/L\ntj1IeptManomIYERUuxU49KYH1JIG3Zp83WTMwkJjJBU8xuXxvyQQtL75fwm5xESGCHZFVZI\nG/fO/mzr8wgJjJDi1V6v+o1HPk5rflgh6Qw5feuzCAmMkGI90r88xmBw7DukS9KZH1pINYdF\n/rvVWYQERkiqL0T6rYjv4r9w+rAWX5+TTKGFpLPkpK3OISQwQlI9cbfi2E3sfnSNbtw/L435\n4YWkP5TXtjyDkMAISbVH/C1O5sr9sT+v6ZPG/BBDeksO2/JF54QERkiq2Y/F/viTLIn9eV9u\nGvNDDEl/Is9t8TEhgRGSatdHY3+c0jd+8s/t0pgfZkgfZe2/xRePkMAISfWAi2Nf7Lbnxk+e\nv1sa88MMSc+SRxt/SEhghKR6bccZX54qs2Onvug4Lo35oYb0Re4ujR8JhARGSKrRnUXk7NiJ\np3t0arKLHyjUkHS83NPoI0ICI6RYhX++/J/xH1r9fff/pDM/3JC+7dC30cOakMAIqVHrW/q9\no8kVbkh6mdzc8AEhgRGSXSGHFO3SvWHThARGSKpDGpfG/JBD0ny5ZvNpQgIjpNip3EOGDD68\nvjTmhx1SWa9OqzedJiQwQlK9sEvuiQ+Wpj8/7JD0VvnNppOEBEZIscofPTGr/dmz0n3shB5S\nxQ/aLas/SUhghFTXiil7S5/LPkxrfugh6b2yaYc0IYER0ubendhDBt6axvzwQ9q416YXnRMS\nGCE1VPHsIImkMT/8kPQxObPuBCGBEdKm3r1oexnwh6a/mCr5MgBSzUGR92tPEBIYIdX23S37\nSucL3timj2yo7Xk5ufZvQgIjJNXKp4fnZJ3wSNM7dXbeXGR+JkDSY+U/8b8ICYyQVHtmH3nD\nVxX1NVpRPGbktgdpjtQe3EFIYIQUO5WV0+yvLL7xgTHbHiT9sbyghARHSC3+yuK3x1Vsi5A+\nzDqgmpDgCKmlysbO10aQNlQkrKQs8Zq0K48WezxhpPy9oqJoncdT4pVF1/gwpaLAjyFron7c\n++uKfBhSURQtT7im4T+NRJCmTtXGkEqi20r/zRmwKujbwFyv4W0eE0CaP7Z0C0iV6xO2pizx\nmrSriBZ7PeI8+ev6onKvp8RaGy31Ycr6Aj+GrImu9WFKeZEPQ9YXRysSrml4IpsA0k0jzznn\nnOFn3JDg/63GZcb3SKor2vct5/dIYPweqYVK4/9/nfsy8qjNFEj6a/kzIYERUmttiz+1ixVt\nn5PT5xLv7zNCQgsrJKyMgfTP2j1qxyFfzZQiJDRCMssHSNU71O2bfsLrQYSERkhm+fHDhvqD\nPK7yehAhoRGSWT5AKo7UQfqj14MICY2QzPLje6Tjah21Te8ddZOIkNAIySw/IC3ZMQ6pj+eP\nckJCIySz/ICkZbePvXK4jPd6DCGhEZJZvkCKH/1dPlCmezyEkNAIySzfIOmizh0/9XYIIaER\nkln+QdLHZF9vH+iEhEZIZvkISS+UCzwdQkhohGSWn5AqDpKHvRxCSGiEZJafkPSLLu0XeDiE\nkNAIySxfIenjsofB23S0FCGhEZJZ/kLSn9e+S7VHERIaIZnlM6T1B8sDng0hJDRCMstnSPrl\ndu3mezWEkNAIySy/IekM2d2riYSERkhm+Q5JJ256qxfzCAmNkMzyH9L6QTLNmyGEhEZIZvkP\nSZd2a/e+J0MICY2QzAoAks6M7ObJHUhIaIRkVhCQdJKM8mIIIaERklmBQKo8Qv7iwRBCQiMk\nswKBpMu6t51nP4SQ0AjJrGAg6XORXe3vQ0JCIySzAoKkl0peem9V3UyEhEZIZgUFqWqI3GE9\nhJDQCMmsoCDpNz1y3zYeQkhohGRWYJD0hchOBbZDCAmNkMwKDpJeIcNsv00iJDRCMitASFVH\ny59NhxASGiGZFSAkXd4z9y3LIYSERkhmBQlJZ2f/IGo4hJDQCMmsQCHpVXKK4bdJhIRGSGYF\nC2njULnJbgghoRGSWcFC0u92yJljNoSQ0AjJrIAh6WvZ/cy+TSIkNEIyK2hI+j9ystXbnRMS\nGiGZFTik6h/JDUZDCAmNkMwKHJJ+1yd7ts0QQkIjJLOCh6T/ye69ymQIIaERklkOQNJr5bjW\nLk46QkIjJLNcgFR9olxnMYSQ0AjJLBcg6eods14xGEJIaIRklhOQ9I2c3ivTH0JIaIRklhuQ\n9I/yw/S/TSIkNEIyyxFI1T+W/LSHEBIaIZnlCCQt7J/1UrpDCAmNkMxyBZLOze31bZpDCAmN\nkMxyBpL+SY6tSm8IIaERklnuQKo5VX6f3hBCQiMks9yBpEU7Z81KawghoRGSWQ5B0nfadFua\nzhBCQiMks1yCpLfI4Mo0hhASGiGZ5RSkmhEyOY0hhIRGSGY5BUmLd4k8k/oQQkIjJLPcgqTv\ntum6JOUhhIRGSGY5BkmnymEpP4IICY2QzHINko6U36Y6hJDQCMks5yAVD4g8neIQQkIjJLOc\ng6QL2nf9OrUhhIRGSGa5B0nvkkNTexAREhohmeUgJB0tv05pCCGhEZJZLkIq2yvyZCpDCAmN\nkMxyEZJ+1H77xSkMISQ0QjLLSUh6jwyswIcQEhohmeUmJD1XLsaHEBIaIZnlKKSyveVReAgh\noRGSWY5C0o87dPocHUJIaIRklquQ9D45oBy8CiGhEZJZzkLSsfIL8BqEhEZIZrkLae0+8nfs\nGoSERkhmuQtJP+nQ6TPoCoSERkhmOQxJH5X9IBmEhEZIZrkMSc+T8chyQkIjJLOchlQxUKYD\nywkJjZDMchqSLurc8dPkVxMSGiGZ5TYkfUz2TR4HIaERklmOQ9Jxcn7SawkJjZDMch1SxUHy\nULJrCQmNkMxyHZJ+2aXdgiSXEhIaIZnlPCR9XPYoTW4lIaERklnuQ9JfyFnJLSQkNEIyKwSQ\n1h8s9ye1kJDQCMmsEEDSL7dr90Ey6wgJjZDMCgMknSG7J3MrCQmNkMwKBSSdKGcksYqQ0AjJ\nrHBAqhws9yReRUhohGRWOCDp0m5t30+4iJDQCMmskEDSmZHdEt7phIRGSGaFBZJeIqMSLSEk\nNEIyKzSQKo+UuxIsISQ0QjIrNJD0m+5t57W+gpDQCMms8EDS5yP9W38MExJaBkCqqkxY6brE\na9JufbTEhymVxesNNvIb+cmG1i5fFy0zmJKwAj+GlEZ9ufeLfRhSWRxN4t5PFVJ5WcKKShKv\nSbvSaKEPU8oKSw02UjxYbpcCYzYAABbqSURBVGrt8pJoscGUhBX4MaQ4usaHKaX+3PnRxPd+\nw3MJPrVrLYundqrLe+S+3crFfGqHlgFP7ZKIkJr0YtZOBS1fSkhohGRWuCDpZBlW0+KFhIRG\nSGaFDFLV0XJrixcSEhohmRUySLq8Z85bLV1GSGiEZFbYIOns7B9EW7iIkNAIyazQQdKr5JQW\nvk0iJDRCMit8kKqHyp+av4SQ0AjJrPBB0u/65LzZ7AWEhEZIZoUQkr6W3e/75s4nJDRCMiuM\nkPQaGdrc9ggJjZDMCiWk6h/JH5s5m5DQCMmsUELS1X2zZzc9l5DQCMmscELS17N7r2pyJiGh\nEZJZIYWk18lxTTZJSGiEZFZYIVWfJNdtfR4hoRGSWWGFpN/vmPXyVmcREhohmRVaSPpGTu9v\ntzyHkNAIyazwQtIb5IdbbpWQ0AjJrBBDqhku12xxBiGhEZJZIYakhf2zXmr8MSGhEZJZYYak\nc3N7rWj0ISGhEZJZoYakf5Jjqho+IiQ0QjIr3JBqTpWrGj4iJDRCMivckLRo56xZmz8gJDRC\nMivkkPSdNt2WbjpNSGiEZFbYIemtMnjT740mJDRCMiv0kGpGyBX1JwkJjZDMCj0kLd4l8kzd\nKUJCIySzwg9J323bdUntCUJCIySzMgCS3i6H1T7qCAmNkMzKBEg6Ui6L/0VIaIRkVkZAKt0z\n8pQSEh4hmZURkPTD9tt/TUh4hGRWZkDSv8ihGwgJjpDMyhBIOlouISQ4QjIrUyCV7SX/IiQ0\nQjIrUyDpR+07Dx906kyvxyghwRGSXd5D0okS7w+ezyEkNEKyy3tIlb1qIclCrwcREhoh2eU9\npPl1juQerwcREhoh2eU9pA8ICY+QzMoYSBt61EGa5vUgQkIjJLt8+GHDU7WOsiKTvH4EEhIY\nIdnlAyR95/R9TnrinQFy6JfeziEkMEKyyw9IdUc2rDlbujzi6RhCAiMku/yDpDq9o4wp83AM\nIYERkl1+QtKFB8qeH3g3hpDACMkuXyFpxaRIbj5y50EREhgh2eUvJNVnusuPmr7JrE2EBEZI\ndvkNSb87UXrNamVtGhESGCHZ5TskrZma69EuJUICIyS7/Iek6tUuJUICIyS7goDk1S4lQgIj\nJLsCgeTRLiVCAiMkuwKC5MkuJUICIyS7goLkxS4lQgIjJLsCg+TBLiVCAiMkuwKEZL5LiZDA\nCMmuICFZ71IiJDBCsitQSMa7lAgJjJDsChiS6S4lQgIjJLuChmS5S4mQwAjJruAh2e1SIiQw\nQrLLAUhmu5QICYyQ7HIBktUuJUICIyS73IBks0uJkMAIyS5HIJnsUiIkMEKyyxVIFruUCAmM\nkOxyB1L6u5QICYyQ7HIIUtq7lAgJjJDscgpSmruUCAmMkOxyC1J6u5QICYyQ7HIMUlq7lAgJ\njJDscg5SGruUCAmMkOxyD1Lqu5QICYyQ7HIQUsq7lAgJjJDschJSiruUCAmMkOxyE1Jqu5QI\nCYyQ7HIVUiq7lAgJjJDschZSCruUCAmMkOxyFxK+S4mQwAjJLpchobuUCAmMkOxyGhK4S4mQ\nwAjJLrchYbuUCAmMkOxyHRKyS4mQwAjJLuchAbuUCAmMkOwKAaSkdykREhgh2RUGSMnuUiIk\nMEKyKxSQktylREhghGRXSCAltUuJkMAIya6wQEpmlxIhgRGSXaGBlMQuJUICIyS7QgQp4S4l\nQgIjJLvCBCnRLiVCArOFVHjzuWdcuQiZT0hgZpBa36VESGC2kH4zefHKW0ZXAFskJDA7SK3u\nUiIkMFNIpTd8o/p93hfAFgkJzBBSa7uUCAnM/nukhacWAVskJDBTSC3vUiIkMHNIpb98aPPp\nkihzvO//mBsZ/23Qt2IbqeF/mMSQlk+4u2bzB2XFCSssSrwm7YqiBT5MKS7w53MptN3g7J3l\n4PebnBu1HdJ8hVFfvmL+3PlJfMUanhclhLTgnOcS/6fVOD61A7N9ahev2V1KfGoHZvvU7tOz\n54HzCQnMHlKzu5QICcwU0obxj8WfCvLH3x7mBaRmdikREpgppAV5tT0PbJGQwDyBFN+l1HZK\n40cCIYHxECG7QgypyS4lQgIjJLtCDWmrXUqEBEZIdoUb0pavUiIkMEKyK+SQtniVEiGBEZJd\noYfUaJcSIYERkl3hh9SwS4mQwAjJrkyAVLdL6fNLT7l0oadjaiMkuwgJzGtIta9SyhGRtk95\nO0cJyTJCAvMckupjEYnX1fMHICHZRUhgPkB6Sep60etBhGQXIYH5AGlmPaTHvR5ESHYREpgP\nkFbk1EHqnl/s7SBCsouQwHyApNfWOjp2e+k8CXnfWThCsouQwPyAVPPQId0OebB6zZTu0nHS\ncu/mEJJdhATmByTdtEO2bOqO0mYM9IsKkQjJLkIC8xVSbNz03SVr1KfeDCEkuwgJzGdIsS/e\njL0lMuy/XgwhJLsICcx3SLEv38xDRYbMtB9CSHYRElgAkGLNGRqnVNP82pQjJLsICSwYSDFK\nwyJywHTbz5CQ7CIksKAgqS4Yky27Tqs0HEJIdhESWHCQVL+ckCP9p5abDSEkuwgJLEhIqksm\ntZde+VaPTEKyi5DAgoWkujp/O+ky2eYFtIRkFyGBBQ1JtSC/m3SatMJgCCHZRUhgwUOKHznU\nV9qMafVd0ZOKkOwiJDAXIKmun7aT5I75LM0hhGQXIYG5AUm1cvpekjXs3bSGEJJdhATmCqT4\nkUOHiAyZncYQQrKLkMDcgRTrlSPSOnKIkOwiJDCnIMWPHBIZmOqRQ4RkFyGBOQZJ9YNREdl3\nelUqQwjJLkICcw6S6sdjcmTnVI4cIiS7CAnMQUiqX09qJ73z4S8yIdlFSGBOQlL9bnIH6ZFf\nhF2JkOwiJDBHIalG87tKp0nfIlchJLsICcxZSKqlU/tI2wnfJH8FQrKLkMAchhS7N6f2k9wx\nSb8hDCHZRUhgTkOK//quPSRr2LzkFhOSXYQE5jik+JFDB4uc8HYySwnJLkICcx6Sas3Mw5P7\n9V2EZBchgYUAktYdOXTQjEQH4RGSXYQEFg5Iqv83LCL7JThyiJDsIiSwsEBS/WhMtuwytaKV\nFYRkFyGBhQeS6uJJbWWHKS3fXkKyi5DAwgRJddmkDtKzxTf+IyS7CAksXJBUv89v+Y3/CMku\nQgILG6TYQ6zFN/4jJLsICSx8kOJHDjX/xn+EZBchgYURUv0b/w17f6tzCckuQgILJ6T4G//t\n0+SN/wjJLkICCyuk+EF4h2115BAh2UVIYOGFpE3e+I+Q7CIksFBD2uqN/wjJLkICCzmkujf+\nGzC19m4nJLsICSz0kFS/qnvjv+pph/Y9bpaHc+ojJLMICc5LSPVv/He4xHvY00HxCMksQoLz\nFlL8jf861TqSLp5/OoRkFiHBeQ1J9dY6SPKO14MIySxCgvMe0oP1kCZ85fEgQjKLkOC8h7Sk\nXa2jiMigqZ5OIySzCAnOe0h6e9xRu+emnxCRdqNmpvR2FklFSGYREpwPkPSt846/OP6Ozt9M\n2VWk76QFHo0hJLMICc4PSI12yM6b0Dn+FC/qxRRCMouQ4HyGpFoxY1i2tB02o9J8CiGZRUhw\nvkOKtWLKHiLdJmz9qqV0IySzCAkuCEix5k3qLrLPlNWWUwjJLEKCCwiS6vqZo3Ik+4QZdoez\nEpJZhAQXGKRYK6ceKNJ1whyjKYRkFiHBBQkp1ieTe4rsnb/MYgohmUVIcAFDin2eM0flStYJ\n09P/bAnJLEKCCxxSrKJpQ0S2G/NKoje1SBAhmUVIcC5A0vhTvN4ie+Z/nc4UQjKLkOAcgRRj\n8MqY9pI1ZFpZylMIySxCgnMGUqzi+FO8Lik/xSMkswgJziVIsRbm7yzyg8kpvXSJkMwiJDjH\nIMXuw1fGdBQZNK0UnkJIZhESnHOQYpXUv3QJ/DITklmEBOcipFjLpgwQ2XFy0ze2aCVCMouQ\n4ByFFLsv50zohL10iZDMIiQ4ZyHFKp8Re4rXdtiMJF+dTkhmERKcy5BiLZ+yu0ifSfOTWUtI\nZhESnOOQFHjpEiGZRUhw7kOqfXV6jrRJ+Op0QjKLkODCACnWt1MHJnzpEiGZRUhwIYGk8ad4\nPeNP8Va1uICQzCIkuPBAqn11eq5kt/jSJUIyi5DgwgQpVuG0g0W2nzCnueNaCcksQoILGSSN\nv3Spl8he+UuaXEBIZhESXPggxV+6NKpN/KVLa7c6m5CsIiS4MELS+lenb/XSJUIyi5DgQgop\n1mf5/UV2mrx48xmEZBYhwYUXUu1LlzrEn+LVv3SJkMwiJLgwQ9L6ly61j790ae1Vu3Y+0of3\nTickuwgJzDtIsT7P30Wk3+Rjat8X8Cnv5tTnKaSy4oQVFCVek3ZF0QIfphT7MqQoWujHmKgf\nQwqjnt77hbPO61j/RrW9Pf+iFSTxFWv45xyEVJO4tRVJLEq3jdE1PkypKanyYcj66FofptQU\n+jFkbXS9xxNK8uolLfN4UM2a6MbEi1KFlER8agfGp3ZQN9RDOvB6r95Zsz5+j2QXIYH5AGlB\nnaPtckX6T5jp4ThCsouQwHyAVPdfUo9FBTPGdBHpOGx6sUdzCMkuQgLzA5L+36TTbyiKn6h4\nZVI/kewhUz73Ygwh2UVIYL5A2mKH7CdThkREBkyagzzqk4qQ7CIkMP8hxVo6bVgbkZ5jZqT+\na/mbi5DsIiSwQCDFWjtzzPYi7YdNW2k3hZDsIiSwoCDFz5wzeU+RrEH584ymEJJdhAQWIKR4\ni6cOyYp/w/RKgl9DlFSEZBchgQUMKdb300d1Euk2anraDxBCsouQwIKHFKt85oQ+Iu1OmLo8\nrSmEZBchgTkBKVb1vPxBIrLP5GZ/eUpyEZJdhATmCqR4i6eekNZhRIRkFyGBuQQpVlqHERGS\nXYQE5hgkTecwIkKyi5DA3IMUL7XDiAjJLkICcxOSpnQYESHZRUhgzkJS/DAiQrKLkMBchqTg\nYUSEZBchgTkOKV7ShxERkl2EBBYCSJrsYUSEZBchgYUDkiZ1GBEh2UVIYKGBpIkPIyIkuwgJ\nLEyQ4rV2GBEh2UVIYGGDpK0cRkRIdhESWAghaUuHERGSXYQEFk5I8ZoeRkRIdhESWHgh6VaH\nEVXeduzAcUuSvzIhtRYhgYUakjY6jOjbU+K/z7XTwqSvSUitRUhgYYekmw4jitT9hvEfJn01\nQmotQgLLAEjxPry+Zx2k7KQfAITUWoQEliGQVH9eBymHkEwiJLCMgfTPOkhDk74CIbUWIYFl\nDKSa4XFHXRYlfQVCai1CAssYSFp199BDJgK/GI+QWouQwDIHEnfIWkZIYIRkFyGBERIaIZlF\nSHCEBEZIdhESGCHZRUhghIRGSGYREhwhgRGSXYQERkh2ERIYIaERklmEBEdIYIRkFyGBEZJd\nhARGSGiEZBYhwRESGCHZRUhghGQXIYEREhohmUVIcIQERkh2ERIYIdlFSGCEhEZIZhESHCGB\nEZJdhARGSHYREhghoRGSWYQER0hghGQXIYERkl2EBEZIaIRkFiHBERIYIdlFSGCEZBchgRES\nGiGZRUhwhARGSHYREhgh2UVIYISERkhmERIcIYEFDanKjwdfTYUfd5huqPFhSHVFlQ9T1I9/\n37SqAnnwpVqNL3d+ZQVy79tDYmwbjJAYM4iQGDOIkBgziJAYM4iQGDOIkBgziJAYMyickApv\nPveMKxcFfSvMmp03N+ibYNML40Zc/G7QN8Km5deNPvPKT5NfH05Iv5m8eOUtoyuCvhlGFY8Z\nmRmQZo99b/Uz43054MnrasbfuW79I6NKk75CKCGV3vCN6vd5XwR9O4y68YExmQFp/KtB3wKz\nSvIWqhblJf+sJ5SQalt4alHQN8Gmt8dVZAakgrxXf3X6ZQuDvhk2XX5bacU/xiV/UF9oIZX+\n8qGgb4JNZWPna2ZAWpT3u+Wl957ly6HZnlc4MS9v7FfJrw8rpOUT7vbjyGwfmjpVMwbSAtWN\nZ88O+nZYVHXJnSXrnhid/JOekEJacM5zQd8Eo+aPLc0USNG8L2N/Tnwi6Nth0QfD4z/KumBm\n0lcIJ6RPz54X9E2w6qaR55xzzvAzbgj6dhhUPTb2r9uGM94M+nZY9H5e/IePYzMc0obxj0Vj\nZcSPv0vjn8m5L/vygl+ve2L0/OgdYzPiblk39s6yDU+OXJn0FUIJaUFebc8HfTvMyoyndlo9\nfcyIK78J+lbYtPTa0Wdd8VHy60MJiTHXIiTGDCIkxgwiJMYMIiTGDCIkxgwiJMYMIiTGDCIk\n78uXwXUH2A4a2sqqIXsmubmqMR3aL6/fcLzOxzzZ4tozO6oenuyGWRoRkvfFHu/Tak/YQHpe\nRs9cV7/h391337Srd5KpLa2NQ7qt8XF883mHexO/rt6X3+7krt/HT9hAelg2HRaaL7WHFpX2\n79zS8W1xSFt0B+9wb+LX1fvy5Yt2P4ufiEMaODB+6tTuqkcf9eah7freVDm5b6ehi2OQ9nr/\nqA5dxxbHLn39hM7tD3ogdmLI0c/1O6J+Ky8e3andvrfW6ND407kl9RuuO0bvUnln09LNV625\ntl/b/Z7Y/NTu5WM69R71pZ4Uu+6ghm3VX2vluJ3a9v5phry0NaAIyfvyZf218oZuDWlov+Pe\nXz5CTrh2xRtdfhJ7TPfb86anL4/kqc7OPua5l38ht6gef8Bef6k/NPfpyI+fmX2pXK6L8uX+\n9zbUb7gO0tXyVv3Shqv+SUa/8vh+e9ZDejly4iMPDOiz6otT5b3PGrZVf63BO9z/2qP798qI\n31oSVITkfflSsX73fSqbQJIFqnPkyNiHo2OP9yHyr9ipc2SZHrRb/CE9PPZ8bag8tWkje+0U\nx3NaboE+JHM2b7gO0lE5JfVLN1+1pu9+sRMrc+shHbJLleo7bW7XC2WLbdVea41cGfvwqxu+\n9esLkokRkvfFIOlLMqUJpPj3L1/Jb2N//lZKdUjbGLWYkqdWyyUVse6Rd3Vom8r6bXwrv4j/\n9YA8vwWkF1atWvnuBXKR1i1tuOoymRRfcUQdpAL5Zf1V4pAabav2WpXd+8/24/3BMjpC8r44\nJB3VYenWkPrHTiyRG2N/TpZiHbJL/IJZMm2+1PeUDu27aRvvyvXxv16Ue7eAVFvOL9dr3dKG\nq75Tt35kHaSP5Q/1V4lDarStugFv7SLdRz7qy/sGZmyE5H21kFZ0Hq6HtApp1/gFL8p98+WC\nubVF65bU9p5cG//rBbl/C0i3zZr17znxn0/ULW246n/rsJxWB+lTuab+KnFIjbZVP2Djq7/d\nRw4p9/KLkOkRkvfVQtJb5dkjYpAOin/vooc3A6lj/OnVgzKzUH626ZoNkFbJhPhf98q/m/ke\nafPShqsulonxvw6sg1Qq58c/Wvp9LaRG22oYoHfLw2af8TYYIXlfHaSq/fsfGYN0fI+a2Pcy\n7ZuBJC/HTp2WtUoP2y7+f8z031c1fpzv1ze+kR93WNMapIarVvfYNeZyUaT+hw379yxVXRh7\ngjdOqhpvq/Za885crfHv1m72+OuQ0RGS99VB0rcikRikqXLjdx8ct29TSEf02+OvsyfL2aqv\n5x4w/aWrc89r/D+SvpB14rOzLoqvbg1Sw1X/R3765F/7D6qH9HzW4f+4d49eq/QaufZfjbZV\ne61VnQ944JV/HtkF+HWIbOsIyfvqIen5EoO04dId2w58bmLnJpAOPmLeUe27jiuLfTjnR51z\n97ipagtI+vJRHdse9KC2Dqnhqhuv3KHN/k9f3KZ+h+wLgzv0GvGF6vKDcvdstK26a304oldu\n3xEfePk1yPgIiTGDCIkxgwiJMYMIiTGDCIkxgwiJMYMIiTGDCIkxgwiJMYMIiTGDCIkxgwiJ\nMYP+H6NH3ZNXbf2AAAAAAElFTkSuQmCC"
          },
          "metadata": {
            "image/png": {
              "width": 420,
              "height": 420
            }
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "plot without title"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3deWAU9fnH8ScX9yEeqJTiURXv\nC2+sVrG2tYJaiq0g1CrS1qu1HthatbT+lKr1qPVAa61WWkWtJ17gbaUiCp6IiqgoqJuLBBLI\n9fx2kkACZLPz2Xlmdr7D5/1H2Ox+Z54luy9I2GFHlDEWOMn3HWAsCRESYwYREmMGERJjBhES\nYwYREmMGERJjBhESYwbZQ6pZZb7L9Wsqr45gilY1RDCkrrwmgilaGcWQ2vK6CKY0VkUwRKvL\nG4HV9pCWrzTf5fo1ppZFMEUro4C0KrUigilaFsWQFako/hhtiOQPhWUpQrKKkMAIyS5CAiMk\nNEIyi5DgCAmMkOwiJDBCsouQwAgJjZDMIiQ4QgIjJLsICYyQ7CIkMEJCIySzCAmOkMAIyS5C\nAiMkuwgJjJDQCMksQoIjJDBCsouQwAjJLkICIyQ0QjKLkOAICYyQ7CIkMEKyi5DACAmNkMwi\nJDhCAiMkuwgJjJDsIiQwQkIjJLMICY6QwAjJLkICIyS7CAmMkNAIySxCgiMkMEKyi5DACMku\nQgIjJDRCMouQ4AgJjJDsIiQwQrKLkMAICW2DgFT1zKwIJBESHCGB5RfSdb1EvvaQ8U7Xj5Dg\nCAksr5DuF68eb9nudf0ICY6QwPIKaf9mSDLedq/rR0hwhASWV0ibt0AaZrvX9SMkOEICyyuk\nvVogjbPd6/oREhwhgeUV0q0tkF623ev6ERIcIYHZQlr8hzE/uuAd//s7v0va0dahP/8ICY6Q\nwEwhNZ16/YqVd40CTsf+8T/vP1r+CtyDnCIkOEICM4VUOXy+avnwBcAel69cslGfz4ANcomQ\n4AgJzPZbu/Ouqar91/g1X526lVlbVr3yavlh9nWBqk1VhDyhufKaCIYsT1VFMGVlaRRDlqWW\nRzClpjyCISsrUrVZ17T9sZENUtnpw4eP+3DNp5UpP325j0z1tZAxpyv3C6n+l9dXrrh3zJr1\nK2uyVllVUzO7eFAq+8oArUiVh7r/1sqXRzCkOrUsgik1pVEMqUxVRzBleTQPfmpF1jW1fiG9\nPsJbevLDWf7eal/z0d9ny0RgEzz+jATHn5HATH9Gem249yCPgyGt2Lp4LrANHCHBERKYKaQV\n466vXnX/yCXAHlv+P9J02Q+5G2iEBEdIYLb/avfxpDE/Pv9NZH7rf+w7Vm5GtgIjJDhCAovJ\n/5D9tHffz833vSZCgiMksJhA0qvlBPN9r4mQ4AgJLC6QGvaWR813vjpCgiMksLhA0leLtl5u\nvvfWCAmOkMBiA0nPkt+a7701QoIjJLD4QKoaWDzPfPctERIcIYHFB5LeJwc3me+/OUKCIySw\nGEHSEXKr+f6bIyQ4QgKLE6RPem38pfkAL0KCIySwOEHSK+VE8wFehARHSGCxglS/l8wwn6CE\nlEOEBBYrSDq7aLvaTCsDREhwhAQWL0h6mlxsPoKQcoiQwGIGadnXurxrPoOQ8AgJLGaQ9B45\nxP7FJEKCIySwuEHSo+Xv5kMICY6QwGIH6eOeG39lPYSQ4AgJLHaQdLKcZD2EkOAICSx+kOr3\nLHjaeAghwRESWPwg6SuFOxi/mERIcIQEFkNI+jP5ve0QQoIjJLA4Qqoc0HW+6RBCgiMksDhC\n0n/JoaYvJhESHCGBxRKSfl/utBxCSHCEBBZPSIt6bpIyHEJIcIQEFk9IeqmcYjiEkOAICSym\nkFbtXPCM3RBCgiMksJhC0hcLBme4JYcICY6QwOIKScfLpWZDCAmOkMBiC6msf9f3rIYQEhwh\ngcUWkt4p37J6MYmQ4AgJLL6Q9AiZajSEkOAICSzGkN7vtnl5xhuhCAmOkMBiDEknyQSbIYQE\nR0hgcYa0aqfCl0yGEBIcIYHFGZI+X7BLncUQQoIjJLBYQ9KT5HKLIYQER0hg8YZUuln3hQZD\nCAmOkMDiDUlvl+8YDCEkOEICizmkpsPl7uBDCAmOkMBiDkkXdNuiIvAQQoIjJLC4Q9KL5bTA\nQwgJjpDAYg9p5Y6FLwcdQkhwhAQWe0j6XMFuQV9MIiQ4QgKLPyQdK1cEHEJIcIQE5gCk0k17\nfBRsCCHBERKYA5D0VvlesCGEBEdIYC5AajpM7gs0hJDgCAnMBUj6XtctA72YREhwhATmBCS9\nUM4MMoSQ4AgJzA1IKwcXzgowhJDgCAnMDUj6pOwe4MUkQoIjJDBHIOlouTr3IYQER0hgrkD6\nol+PRTkPISQ4QgJzBZLeLEflPISQ4AgJzBlIjUPlgVyHEBIcIYE5A0nfKtky168IIcEREpg7\nkHSi/CrHIYQER0hgDkGq2bZoTm5DCAmOkMAcgqRPyJDcnqqEBEdIYC5B0uPlupyGEBIcIYE5\nBWnpRr0X5zKEkOAICcwpSHqDjMxlCCHBERKYW5AaD5SHchhCSHCEBOYWJH2zZFA1PoSQ4AgJ\nzDFIeq6ciw8hJDhCAnMN0optil+HhxASHCGBuQZJH5N94ecrIcEREphzkHSk3IBuQkhwhATm\nHqQlfft8Bm5CSHCEBOYeJL1Ojge3ICQ4QgJzEFLjAfIIuAUhoRESmIOQ9I3irZZDGxASHCGB\nuQhJz5aJ0HpCgiMkMCchrdi6eC6ynpDgCAnMSUg6XfZD7jYhwRESmJuQ9Fi5GVhNSHCEBOYo\npE979/3c/2pCgiMkMEch6dVygv/FhARHSGCuQmrYWx71vZiQ4AgJzFVI+mrR1r5fTCIkOEIC\ncxaSnikX+l1KSHCEBOYupKqBxfN8LiUkOEICcxeS3icHN/lbSUhwhATmMCQdIbf6W0hIcIQE\n5jKkT3pt/KWvhYQER0hgLkPSK2Wsr3WEBEdIYE5Dqt9LZvhZR0hwhATmNCSdXbRdrY9lhARH\nSGBuQ9LT5GIfqwgJjpDAHIe07Gtd3s2+ipDgCAnMcUh6jxyS/cUkQoIjJLBQIdWtzNqy5dnX\ndNpRcmvWNbWpioBTfFVeE8GQ5amqCKasLI1iyLJU0EffTzXlEQxZWZGqzbqm7Y8NEFLt8qxV\nLMu+ptPe7bnxx9nWVKfKA07xVVl1BEOWpSojmLK8NIohFamgj76fqssiGLK8PJX90a/JFZKP\ngn5rp3q5nJRtCb+1g+O3dmCu/4ykWr9nwdNZlhASHCGBuQ9JXyncIcuLSYQER0hgCYCkP5NJ\nnS8gJDhCAksCpMoBXed3uoCQ4AgJLAmQdKoc2umLSYQER0hgiYCkR8mdnd1MSHCEBJYMSIt6\nbpLq5GZCgiMksGRA0kvllE5uJSQ4QgJLCKRVOxc8k/lWQoIjJLCEQNIXCwZn3hMhwRESWFIg\n6SlyacbbCAmOkMASA6msf9f3Mt1GSHCEBJYYSHqnfCvTi0mEBEdIYMmBpEfI1Ay3EBIcIYEl\nCNL73TYv7/gWQoIjJLAEQdJJ8rOObyAkOEICSxKkVTsVvtThDYQER0hgSYKkzxfsUtfR9YQE\nR0hgiYKkJ8nlHV1NSHCEBJYsSKWbdV/YwdWEBEdIYMmCpLfLdzq4lpDgCAksYZCaDpd71r+W\nkOAICSxhkHRBty0q1ruSkOAICSxpkPRiOW296wgJjpDAEgdp5Y6FL697HSHBERJY4iDpcwW7\nrftiEiHBERJY8iDpWLlynWsICY6QwBIIqXTTHh+tfQ0hwRESWAIh6a3yvbWvICQ4QgJLIqSm\nw+S+ta4gJDhCAksiJH2v65ZrvZhESHCEBJZISHqhnNn+U0KCIySwZEJaObhwVrtPCQmOkMCS\nCUmflN3r2z4jJDhCAksoJD1Brm77hJDgCAksqZC+6Ndj0ZpPCAmOkMCSCklvlqPWXCYkOEIC\nSyykxqHywJrLhIRGSGCJhaRvlXy9qvUiIcERElhyIelEObv1EiHBERJYgiHVbFv0WsslQoIj\nJLAEQ9InZJ+WpzYhwRESWJIh6fFyXfOvhARHSGCJhrR0o96LvV8JCY6QwBINSW+Qkd4vhARH\nSGDJhtR4oDykhJRDhASWbEj6ZsmgakLKIUICSzgkPVfOJaQcIiSwpENasU3x64SER0hgSYek\nj8m+DYQER0hgiYekI+UGQoIjJLDkQ1rSt8+nhIRGSGDJh6TXyXb7ffumCJ7khARGSHZFAGlB\nkaQb0RT6IEICIyS7IoD0bWnu36EPIiQwQrIrfEiNJS2Qxoc9iJDQCMkuQgIjJDRCMurwFkj/\nDH0QIYERkl0RQJrfy3N0KP+xAYmQwDYASLrolF23l5PCn0NIYIRkVxSQtDFVvmPR/NDHEBIY\nIdkVEaRlU+WE0McQEhgh2RUVpMY9CueFPYaQwAjJrqgg6QNybNhjCAmMkOyKDJLuL/8LeQwh\ngRGSXdFBekKODHkMIYERkl3RQdJvybPhjiEkMEKyK0JIL8rQcMcQEhgh2RUhJD1Sngh1DCGB\nEZJdUUKaUzAk1OOECAmMkOyKEpIe23busTAiJDBCsitSSG8V7oL8btEICYyQ7IoUkp4gU0Mc\nQ0hghGRXtJDeL96+PrwxhARGSHZFC0lPlr+FN4aQwAjJroghfdJ1UHgDCQmMkOyKGJKeLteH\nNoaQwAjJrqghLemxRWhPREICIyS7ooak58mfwhpDSGCEZFfkkEr7bBLWG4ETEhgh2RU5JL1Y\nfh/SGEICIyS7oodUuXHf0nDGEBIYIdkVPSS9XCaGM4aQwAjJrjxAWr5F989DGUNIYIRkVx4g\n6bVyVihjCAmMkOzKB6SVX+/yURhjCAmMkOzKByS9RU4NYwwhgRGSXXmB1DC46L0QxhASGCHZ\nlRdI+k8ZE8IYQgIjJLvyA6lx98I37McQEhgh2ZUfSHqf/MB+DCGBEZJdeYLUtF/BK+ZjCAmM\nkDI2ffxxZ8xG5ucJkk6X75mPISQwQsrUzHGvfvngqcgjnS9Ieog8Zz2GkMAIKVOnPo3Ozxuk\nF+Rg6zGEBEZIGSod/vSZPzyn7RyT9XVZq1qRfU3gVqYq179ymDxuPKZipfEOO2pFqjqCKXVl\nUQypSkXy6FdEMKSuIuXj0fcLacHw3yyuuuXHa/4EqEzFuKcK9vwq3/eBbVCV+4c0L/036Qkz\nV39euzxrFcuyrwlcdaq8g2uPkntsx5RV2+6vw5alKiOYsrw0iiEVqUge/bIIhiwvT2V/9Gv8\nQkoN/yD98fR7O1+1Vnn7GUn1zcLdbN/AmD8jgfFnpAw1jnsk/Ugf/wKwxzxC0uPlbtMxhARG\nSJm6d8zc1F/G1QJ7zCekBcZvYExIYISUqcY7xh53wafI/HxC0pPkdssxhARGSHblFdLHXbay\nHE9IYIRkV14h6S/kBsMxhARGSHblF9KSHlsaPi0JCYyQ7MovJP21XGU3hpDACMmuPENK9d60\nymwMIYERkl15hqQXyh/NxhASGCHZlW9IlRv3NXvOEBIYIdmVb0h6qfzWagwhgRGSXXmHtHzz\nnl8YjSEkMEKyK++Q9Go522gMIYERkl35h1Q7sBt0UFPmCAmMkOzKPyS9SX5mM4aQwAjJrhhA\nqvtGyYcmYwgJjJDsigEk/YeMMxlDSGCEZFccIDXsVPSOxRhCAiMku+IASe+RURZjCAmMkOyK\nBaSmvQpeNxhDSGCEZFcsIOkj8n2DMYQERkh2xQOSHigvBx9DSGCEZFdMIM2UQ4KPISQwQrIr\nJpD0cIHftny9CAmMkOyKC6RXCvZtCjqGkMAIya64QNKj5NGgYwgJjJDsig2kN4K/gTEhgRGS\nXbGBpD8U5D3LO4qQwAjJrvhAeq94h4BvYExIYIRkV3wg6Vi5I9gYQgIjJLtiBGlRl62DPa6E\nBEZIdsUIkk6QmwONISQwQrIrTpA+7z6gJvuqzBESGCHZFSdI+iu5JsgYQgIjJLtiBemrYG9g\nTEhghGRXrCDpb+T/AowhJDBCsitekCr6bVSefVWmCAmMkOyKFyT9g/wu9zGEBEZIdsUMUnX/\nXl/mPIaQwAjJrphB0ivlnJzHEBIYIdkVN0i1A7stznUMIYERkl1xg6R/ldNyHUNIYIRkV+wg\n1W1bsjDHMYQERkh2xQ6S/l1+muMYQgIjJLviB6lhx6J3cxtDSGCEZFf8IOm/5Ue5jSEkMEKy\nK4aQmvYsmJvTGEICIyS7YghJH5QROY0hJDBCsiuOkPQAmZXLGEICIyS7YgnpSfl2LmMICYyQ\n7IolJP2WPJPDGEICIyS74gnpJRmawxhCAiMku+IJSb8rj+FjCAmMkOyKKaQ5BUPw99QnJDBC\nsiumkPQ4uR8eQ0hghGRXXCG9XbgL/J76hARGSHbFFZKOlrvQTQgJjJDsii2kD0q2qwM3ISQw\nQrIrtpB0vNwKbkFIYIRkV3whfdL1a+AbGBMSGCHZFV9IeqZch21ASGCEZFeMIS3tsVk1tAEh\ngRGSXTGGpOfLZGg9IYERkl1xhlS+0UYVyHpCAiMku+IMSS+RS5DlhARGSHbFGlL1ZtAbGBMS\nGCHZFWtIOlnOB1YTEhgh2RVvSDVf6/6Z/9WEBEZIdsUbkv5FzvC/mJDACMmumENatU2Xj3wv\nJiQwQrIr5pD0VjnF91pCAiMku+IOqWFw0Xt+1xISGCHZFXdIepeM9ruUkMAIya7YQ2rco3Ce\nz6WEBEZIdsUekv5HjvO5kpDACMmu+EPS/eR//hYSEhgh2eUApMflO/4WEhIYIdnlACQ9VJ7z\ntY6QwAjJLhcgvSgH+1pHSGCEZJcLkPTb8qSfZYQERkh2OQHpVX9vYExIYIRklxOQdIQ86GMV\nIYERkl1uQHqrcFcfXyZCAiMku9yApD+Wf2VfREhghGSXI5DeL96+PusiQgIjJLscgaQ/lduy\nriEkMEKyyxVIH3cdlPWeEhIYIdnlCiQ9Tf6abQkhgRGSXc5AWtJjy2zPYEICIyS7nIGk58qV\nWVYQEhgh2eUOpFSfTbLshZDACMkudyDpRTKp8wWEBEZIdjkEqXLjvp0/vQgJjJDscgiSXiYX\ndHo7IYERkl0uQVq+ec8vOrudkMAIyWf1dVmrWpF9TeBWpipN9nOVnNnZzRUrTaZ03opUdQRT\n6kqjGFKViuTRr4hgSF1FysejnyukmuqslVdmXxO4qlSZyX5SA7vO7+TmsiqTKZ1XmaqIYEp1\naRRDKlLLIphSZfPgZ6kslf3Rb/teYsP+1k71ZpnQya381g6M39rZ5Rakum8Ud/IGxoQERkh2\nuQVJ75QTM99ISGCEZJdjkBp2Lnwj442EBEZIdjkGSe+VkRlvIyQwQrLLNUhN+xbMznQbIYER\nkl2uQdJH5ahMNxESGCHZ5RwkPUSez3ALIYERkl3uQXpBvpnhFkICIyS73IOkw2RmxzcQEhgh\n2eUgpNkF+3b8BsaEBEZIdjkISY+Whzu8npDACEm1ou1LMO+aAPNdhPRm4W4dftEICYyQ0peu\nT3+onrgg/fH6IH9PuQhJR8k9HV1NSGCE1AppqTyuGyKkBcU7dPQGxoQERkgbOCQdJ//o4FpC\nAiOkDR3Soi5bdfAUICQwQtrQIenP5cb1ryQkMELa4CF93r2DNzAmJDBC2uAh6dny5/WuIyQw\nQkpfOmfWrFnT5Zr0x3M2REip3ptWrXsdIYERUvpS+wLMdxWS/lYuXfcqQgIjJNVL2hdgvrOQ\nKtZ/A2NCAiMku5yFpH+UC9e5hpDACMmrebvG5+96K9B8dyFV91/3DYwJCYyQ0t21VU2awQHp\nn5B+GWS+u5D0z/Lrta8gJDBCUp1eMPAz1fPklDuOlgcDzHcYUu3AbovXuoKQwAhJ9cjtKtJ3\ncZNvNmnDbsMDzHcYkt4oP1/rc0ICIyTVTS9Lf5glf0t/vHjLAPNdhlS3bcnC9p8TEhghqRb9\nO/3hT7Io/fHWkgDzXYakt8tP2n9KSGCEpNpvavrDUQO8i1d3CzDfaUgNOxW92+5TQgIjJNXd\nz0h/sbs2v6P8T7cLMN9pSHq3HN/uM0ICIyTVST2nfXBM8/tSvd9zfID5bkNq2rPg9bbPCAmM\nkFRTW4vICekLD2za693MG2TNbUj6sBzd9gkhgRFSurKrz7vbe3u3f27/bJD5jkPSA+TlNZcJ\nCYyQ2rWy4zdL9JvrkGbIsDWXCQmMkOxyHZIeJs+svkhIYISkOrR9AeY7D+m/ctDqi4QERkjp\nSyX7DD1g/9YCzHcekn5PprdeIiQwQlI9pU/JkX9f7z9b47kP6bWC3Vu/hIQERkjpaqYeWdj9\nhMeDPnfch6Qj5b6WC4QERkgtfTZ5J9nynMxn+fZTAiC9XTS45Q2MCQmMkNY0+/RNZY/135jK\nfwmApCfKnc2/EhIYIbVV+9AQKQgwPwmQFnXZuvkJQUhghLS62b/YSLb9/UcB5icBkp4qU7xf\nCAmMkJr74qpdpPfJz2/YRzY093n3ATVKSHCEpFr3wIjiwiPuCvygJgKSniXXKiHBEZLqZkUH\nXfZhbWsB5icD0le9N6siJDhCSl8qLN6w37J47S6QywgJjpD4lsXrVNFvo3JCQiMkuxICSSfJ\nRYSERkh2JQVSdf9eXxISGCHZlRRIeoWc8074UwgJj5DMigJSTV+RwuGfhD6HkNAIyawoIF3T\n/O+Xuwd5KcBXhIRGSGZFAKl+o5ZXAv4R9iBCQos9pCGtb8N1304B5icF0uLWl9QuCHsQIaHF\nHpK82vxL/aQuAeYnBVJVUQukq8IeREhoMYfU7lzMeweYnxRI+sPmL0VRsPMX+oiQ0GIOad51\ncswpXuMvXpx5g6wlBlLpkLSjYhlq8D4WnUZIaDGHpPqd9w3mJwaSNj580ZRPfhC6JEJCiz0k\nXbFEteb2qxZmWOyr5EBqPmi14QQ5KNxZhIQWe0jz+1+u9fuI9H0903IfJQySNowOWRIhocUe\n0g92+1D/KTd+eNAPA8xPGiRP0oFhTiMktNhD6j9V9bhdVad+PcD8xEHShjEypDy8IYSEFntI\nXZ7Rhn7nqz7F15Faaj36u+FEGRLe85CQ0GIP6et/06e8UzHctqGe1XzdVv83irSkvUN7IhIS\nWuwhnbLFBVt9o0G/3J0/I7W05v8jNYwNTxIhocUe0pIDZNNZqj/qG+RNixMJSRvGyV6l4Qwh\nJLTYQ0pvXJf+8OoXQeYnE1KIkggJzQFItbP/k9L6QPMTCkkbfiJ7hiKJkNDiD+mq3iKz9Lcn\nBaGUVEjaGJIkQkKLPaRbZMTNaUh3FF8RYH5iIWnTz2XPlP0QQkKLPaTdf661aUj6mx0CzE8u\nJE/SHvaSCAkt9pC6zWiB9GRJgPkJhqRNvwhBEiGhxR5S/0daIE3rE2B+kiFp02my01LjIYSE\nFntIRxxa40Eq2/XIAPMTDSkMSYSEFntIzxZt90s5+Sd9Sl4KMD/ZkLTpdNlxiekQQkKLN6T9\n71WduZf3NgX7PRdkfsIhadMZxpIICS3ekOR67+OXc+cG/B8DSYfkSRpsKYmQ0FyAFLzEQ9Km\nM2Xw53ZDCAmNkMzKKyRtOstSEiGhxRzSyGvaCjB/A4CkTb+UHcwkERJazCG1L8D8DQGSNv1K\ndvjMaAghocUc0lnPthVg/gYByTvH7PZGkggJLeaQ+DPSenV2xr7fyPZB3pC2LUJCIySzYgBJ\nfytbL7IYQkhohGRWHCDphTaSCAkt3pB+9qzN/A0HUlrSVh8FH0JIaPGGZNUGBEl/ZyGJkNAI\nyayYQNKLZFBgSYSERkhmxQWSXiyDAp27QwkJj5DMig0kA0mEhOYopJnDZyHzNzBIeokM+jDQ\nEEJCcxNSxdiRhNRZk+XrgSQREpqbkC6/bSwhdVpa0gcBhhASmpOQXh5f2w5SQ33WqlZkXxO4\nVanKCKbUV6zytexyGfhe7kNqUtW5b+y/siiGVKdqIpiyqiKCIfWVqeyPftsftNkgVY+bq+0g\nVaZYB10sA2bn+z6wPNT2v8mzQbr2Wm0PqaY6axWV2dcEripVFsGU6rIqnwsvlc3n5DqkMlWR\n66ZIpVEMqUgti2BKVTQPfir7o9/2TXkWSHPHVa0FyUcb3s9IXlfKFu/kOIQ/I6E5+DPSFSNH\njx494vjLgD1umJD0Ktk8R0mEhOYgpCrvG8ETn0KetRsoJE/S2zkNISQ0ByE1x2/tfPXnHCUR\nEpqrkLA2WEh6tfR/K4chhIRGSGbFEpJeU5CLJEJCIySz4glJbyro/yY8hJDQCMmsmELSmws2\ngyUREhohmRVXSDoFl0RIaIRkVmwheZLewLYgJDRCMiu+kPSWwn5zoA0ICY2QzIoxJE/Sq8h6\nQkIjJLPiDElvLdwIkURIaIRkVqwh6d8KN5rtfzUhoRGSWfGGpLchkggJjZDMijkkT9IrftcS\nEhohmRV3SPr3wr5+JRESGiGZFXtInqT/+VtJSGiEZFb8Iem/i3xKIiQ0QjLLAUj67+K+vv4j\nFyGhEZJZLkDSu4v7+JFESGiEZJYTkNKSej6XfRUhoRGSWW5A0nv8SCIkNEIyyxFInqRns60h\nJDRCMssVSDotuyRCQiMks5yBpPeW9Hym8xWEhEZIZrkDKS2pR+eSCAmNkMxyCJLeV9Lj6c5u\nJyQ0QjLLJUh6f+eSCAmNkMxyCpI+2rXHzMy3EhIaIZnlFiRP0oyMNxISGiGZ5Rgknd61e0ZJ\nhIRGSGa5Bkkf69b10Qw3ERIaIZnlHCRP0iMd30JIaIRklnuQ9PFuXR/u8AZCQiMksxyElJbU\npUNJhIRGSGa5CEmf6NbloQ6uJiQ0QjLLSUj6ZIeSCAmNkMxyE5I+2b3Lg+tdSUhohGSWo5D0\nqQ4kERIaIZnlKiR9rmeXB9a5ipDQCMksZyHp871K1pFESGiEZJa7kDxJ/1nrCkJCIySzHIak\nL/Qqub/954SERkhmuQwpLanoX+0+JSQ0QjLLaUj6Yu+iqW2fERIaIZnlNiRP0l1rPiEkNEIy\ny3FI+lI7SYSERkhmuQ7Jk/TP1ouEhEZIZjkPSf/bZ7UkQkIjJLPch+RJurP5AiGhEZJZCYCk\nc/oV3eH9SkhohGRWEiB5kv6hhIRHSGYlApLO2bjwH4SER0hmJQOSvrZx4e2EBEdIZiUEUlpS\nwXnH7jjs7pDHeBESGCHZFUonNKMAABiySURBVDokfb23eF0Q9hxCgiMku8KHtLJvMyR5O+xB\nhIRGSHZF8DdSiyO5KexBhIRGSHaFD2luK6QpYQ8iJDRCsit8SHVbtEBaEPYgQkIjJLvCh6TT\nmx11mxf6IEICIyS7IoCkb4zb97gz5Wsfhz2HkMAIya4oIDW/IDtJtv8y5DGEBEZIdkUGSX8p\n+y0PdwwhgRGSXdFBahwl368PdQwhgRGSXdFB0lXfljFNYY4hJDBCsitCSLpsbzk/zDGEBEZI\ndkUJSb/aQa4KcQwhgRGSXZFC0g+3KLg9vDGEBEZIdkULSd/YqOSx0MYQEhgh2RUxJH2ma4//\nhjWGkMAIya6oIendhZvOD2kMIYERkl2RQ9IbZeCn4YwhJDBCsit6SPpb2bU8lDGEBEZIduUB\nUtN4OSCUN0MhJDBCsisPkLThBzI8jIOFCAmMkOzKBySt+aacGsIYQgIjJLvyAkkr95CL7ccQ\nEhgh2ZUfSPr51nKt+RhCAiMku/IEST/oXzjNegwhgRGSXfmCpLN7dXnKeAwhgRGSXXmDpDO6\n9HnNdgwhgRGSXfmDpFMLN7N9iy5CAiMku/IISa+Qbyy1HENIYIRkVz4h6Tmye4XhGEICIyS7\n8gqp6SQ5zPBrSUhghGRXXiFp3ffkR8hj2XmEBEZIduUXkq44SH5hNoaQwAjJrjxD0tId5f+s\nxhASGCHZlW9IuvjrBbcajSEkMEKyK++Q9O1+RffbjCEkMEKyK/+QdFaP7i+YjCEkMEKyKwaQ\n9OHiviYnTyIkMEKyKw6Q9I6CAYsMxhASGCHZFQtI+gfZzuDkSYQEFiqkhvqsVa3IviZwq1KV\nEUypr1gVwZCaVHWWFWfKPhWBx5QF3oOPqlM1EUxZFfyr4aPKVPZHv+0PWhDSiqqslVdmXxO4\nZamyCKZUlS2LYEhlqiLbih/IYaVBxwTegZ/KU5E8+tE8+Knsj351rpB8xG/twLJ+a5decqSM\nDnqwEL+1A+PPSHbFBZJW7S3nBRxDSGCEZFdsIBmcPImQwAjJrvhA0oVBT55ESGCEZFeMIOmb\nAU+eREhghGRXnCDps8FOnkRIYIRkV6wg6QNFmwQ4eRIhgRGSXfGCpDcFOXkSIYERkl0xg6S/\nk11yPnkSIYERkl1xgxTk5EmEBEZIdsUNkjaMzPnkSYQERkh2xQ6Sd/KkcU05jSEkMEKyK36Q\ntHJPuSinMYQERkh2xRBSzidPIiQwQrIrjpC8kyfdk8MYQgIjJLtiCck7edKT+BhCAiMku+IJ\nSWd26Y2fPImQwAjJrphCyunkSYQERkh2xRWSXi/boidPIiQwQrIrtpD0XNkNPHkSIYERkl3x\nhdT0U/kW9oUmJDBCsiu+kLyTJx0L3TtCAiMku2IMSVcMlZ8j6wkJjJDsijMkLd1J/ggsJyQw\nQrIr1pB08SDk5EmEBEZIdsUbkr69cdF9vhcTEhgh2RVzSPq/nt18nzyJkMAIya64Q9JHivvO\n9bmUkMAIya7YQ9I7fZ88iZDACMmu+EPSS/2ePImQwAjJLgcg6a9k3+rsqwgJjpDscgFS01gZ\n5ufZS0hghGSXC5C07ju+Tp5ESGCEZJcTkLRqiJyRfRUhgRGSXW5A0q8Gy5VZFxESGCHZ5Qgk\n7+RJf8+2hpDACMkuVyB5J0+anmUJIYERkl3OQNJnu3V/qfMVhARGSHa5A0kfLNrk3U4XEBIY\nIdnlECTv5EmfdHY7IYERkl0uQdKLZJfOsBASGCHZ5RQkPUP2X575VkICIyS73ILUMFKOznzy\nJEICIyS73IKktYfI2IwnTyIkMEKyyzFIumxPuTDTbYQERkh2uQbJO3nSNRluIiQwQrLLOUj6\nweYFd3Z8CyGBEZJd7kHSV3uVPNHhDYQERkh2OQhJn+7ae05H1xMSGCHZ5SIk/VfhZu91cDUh\ngRGSXU5C0r92ePIkQgIjJLvchKTnd3TyJEICIyS7HIXUdLJ8q3bdKwkJjJDschSSNhwrx6x7\n1wkJjJDschWS1qx/8iRCAiMku5yF5J086Q9rX0NIYIRkl7uQ9LNBcv1aVxASGCHZ5TAkfWed\nkycREhgh2eUyJO/kSc+3+5SQwAjJLqch6SPFfdqdPImQwAjJLrch6T/bnzyJkMAIyS7HIen/\nyXZfrL5MSGCEZJfrkPRs2Wf1yZMICYyQ7HIeUtM4Obz18SAkMEKyy3lI3smTTmh5PhASGCHZ\n5T4kXXGAnN58gZDACMmuBEDS1GD5k/crIYERkl1JgKQfbVlwmxISHCHZlQhI+ma/ogcICY6Q\n7EoGJH2uW/fnHr/i8bqQxyghWUZIYOFD0oeKikRkx/lhzyEkwwgJLAJI9QPFa7fMb7FvFCHZ\nRUhgEUD6r7T0QtiDCMkuQgKLANIDrZDuDXsQIdlFSGARQHqzFdLZ6723kHGEZBchgUUASY9p\ndlQiA6eE+zsiJLsICSwKSOWjC6Rg9AcTu8nO08KcQ0h2ERJYFJDSlJ7xXpH9ZEKhHPRieFMI\nyS5CAosG0uojG946WgpGvR/WEEKyi5DAooWkOmNvKZnQwbvsW0RIdhESWNSQtGnaN6TnxFC+\ngIRkFyGBRQ5JtW7K5rLp5BAeKEKyi5DA8gBJtXpyb9lhWpP1EEKyi5DA8gJJ9bMJxbLvM8ZD\nCMkuQgLLEyTV+aMK5Ih5pkMIyS5CAssbJNVZh0jhqI8MhxCSXYQElkdIqjN2ky5nlZsNISS7\nCAksr5C08Y4B0m+y1cGshGQXIYHlF1L6yT+5r9nBrIRkFyGB5RuSaqnZwayEZBchgeUfUuvB\nrC8FH0JIdhESWBwgWR3MSkh2ERJYPCDZHMxKSJkqu/LE4y9YgMwnJLC4QLI4mJWQMnX2xIVL\nrhqD/OsoIYHFBpJ3MGt/2XRyAAuElKGqyz5V/Wo48r0zIYHFCFLgg1kJqbPmH7Pmpe/GhqxV\n1WRfE7i6VGUEUxoqVkUwpDa1PIIpDaU+131yarHsMzPHIdWp2hy3RFpVEcGQhspUXdY1bdR8\nQKo67fY1lytTLPG9PKJADn0u3/fCidoOrsoOafGEG9v+pl9RlbXyyuxrArcsVRbBlKqyZREM\nqUxVRDClqhRY+/RQKTzurRyGlKciefSjefBT2R/96jUyskKaN/qR7H9ptY8/I4HF6mek1TUf\nzFoBD+HPSJl654Q54HxCAoslpOaDWTeGD2YlpAytOvXf3reC/OfvEIsnpNwOZiWkDM0b3tyj\nwB4JCSyukHI5mJWQ7CIksPhCwg9mJSS7CAkszpBaD2b9wO9qQrKLkMDiDQk7mJWQ7CIksLhD\nQg5mJSS7CAks9pCAg1kJyS5CAnMAku+DWQnJLkICcwJS6zuzPptlESHZRUhgjkDy9c6shGQX\nIYE5A6n1nVkXdbKAkOwiJDCHIGU9mJWQ7CIkMKcgZTmYlZDsIiQwtyB1fjArIdlFSGCuQfIO\nZu2a4WBWQrKLkMDcg5T5YFZCsouQwFyEpPrW9zs6mJWQ7CIkMDchdXwwKyHZRUhgrkLq6GBW\nQrKLkMCchdTBwayEZBchgTkMab2DWQnJLkICcxpSy8Gs+z3b+gkh2UVIYI5DWutgVkKyi5DA\nnIfU7mDWD54MeIIlXxGSWYQEFyak1oNZ3z5SpPDU8H83hGQWIcGFC6n5YNZi8Rof7iAlJMMI\nCS5kSOkHflyzIykI/bs7QjKLkOBCh6Q3tECS/4Y9iJDMIiS48CHd1wrpvrAHEZJZhAQXPqTq\nrVslZXljh8ARklmEBBc+JJ29TVrR/g9/UwpHLQxzDiGZRUhwEUDS2gevntmkOmN36TLhi/DG\nEJJZhAQXBaTVRzY0Tttaevl7k+NcIiSzCAkuSkjp31TzceEhPRUIySxCgosWkmrV5N6yFXjC\nP58RklmEBBc1JNWvJnaRXaAT/vmMkMwiJLjoIakuGlsgB71gPoWQzCIkuHxAUn11mMgRbxhP\nISSzCAkuP5BUZwzJ8n7heIRkFiHB5QuSNk3bXrpM+NJwCiGZRUhweYPkvUfKAOk9scpsCiGZ\nRUhweYTkvV/4RrKZn3Nn+oqQzCIkuLxCSo+f2F22noI8NTNHSGYRElyeIakunlAs+8ywmEJI\nZhESXN4htb7d0KvBpxCSWYQEFwNIqq8cJgWjFgSdQkhmERJcLCA1v/V+8YTPg00hJLMICS4m\nkLRp2nbSY2Lm09D6iJDMIiS4uEDyXlbaUjaeXJP7FEIyi5Dg4gMp/RxpPg1tfa5TCMksQoKL\nEyTvNLTdZMe2s1hgEZJZhAQXL0iqn04okv2ezmkKIZlFSHBxg6T6rvey0ms5TCEkswgJLn6Q\nVGcd2tEZnbNGSGYRElwcIanO2FNKJiwBpxCSWYQEF09I2jhtW+k5EZNBSGYRElxMIXnv3LW5\nbDK5FtiCkMwiJLjYQvLO6NxHvg68cxchmUVIcDGGpJqa2FV28v2yEiGZRUhwsYak+smEIjng\nOX9rCcksQoKLOSTVt0eJHDHXz0pCMouQ4GIPSfW/3glhPsy+jpDMIiQ4ByA1nxCmJPsJYQjJ\nLEKCcwKSNk7bRnpmOyEMIZlFSHBuQPJ1QhhCMouQ4FyB5L2s1FsGdfayEiGZRUhw7kDyTgjT\ntbMTwhCSWYQE5xIk1Y8nFMqBz2e4kZDMIiQ4tyCpvjUq4wlhCMksQoJzDZLqS0OlcNRHHdxA\nSGYREpx7kFQf3qHDE8IQklmEBOciJK2fMkB6rXdCGEIyi5DgnITknRCmn2y2zstKhGQWIcE5\nCkm1fGJ32WqtE8IQklmEBOcsJNXPJhTLru1eViIkswgJzmFIqu+NKpChL67+jJDMIiQ4pyGp\nzj5c5Ii3Wi4TklmEBOc4pPYnhCEkswgJznlIrSeEKdd7fzj09EXhjVkdIdlFSGChQmo5IUy/\noZKu+5ww5zRHSHYREljIkLwTwvSS5nYLd44SkmWEBBY6JNULWiDJ0rAHEZJdhAQWAaRLWiFt\n/7Opn4U6iJDsIiSwCCDNbHHUo0/6w5ajprwd2iBCsouQwCKApCc2Q3q2fs7kozdqxnTtnBzP\n/Nd5hGQXIYFFAanhhkN2GPV6y8U5147aOI1p86Mnz0Ge9b4iJLsICSwKSOu8INv49pRRm6Yx\n9Tli8ot1llMIyS5CAssDpOYWThk7KI2p1xGXzDC7A4RkFyGB5QuS18I7JmyTxtQzjQk52VLG\nCMkuQgLLJySvz6dN2Nk77mHoxBk1QacQkl2EBJZvSF5Lp501pECkeMjEhyuCTCEkuwgJLA6Q\nvL58eOKQQpGiIWdNy/n3TUh2ERJYXCB5Vc2YOLQkjWnnCdNSuUwhJLsICSxOkLyqZ1xyRNf0\nD03bTrjjU3RKqJAaG7JWXZN9TeDqUpURTGmoWBXBkNrU8gimNJRFMaQ6VRvBlFUVyOqq5y8b\n1s3DdOJNHyLbVabqsq5powZCWl6ZtbLy7GsCV5Eqi2BKZWlFBEPKU1F8xSpTUQwpi+T3UlGK\nblH63O+/2zuNaYtjr/2f721S2R/9tnfc47d2ncVv7cDi9q1d++rnXDuqH3BwHn9GsouQwOIM\nqXnLt6eM2iSNqb+Pg/MIyS5CAos7pOYWThk7MI2pd5aD8wjJLkICcwKSVxrTVlkOziMkuwgJ\nzBlIXgvvmLCtd3De0IkdHpxHSHYREphTkLy8g/MKREqGTHx43WcUIdlFSGDOQfJac3DeWdPa\nH5xHSHYREpiTkLy+WvfgvLprDt1j/CL/OyCkziIkMGcheXkH53URKfQOzms6yntjiF7zfW9M\nSJ1FSGBOQ/Ja3npwXv+W9yr6lu8NCamzCAnMeUhey2dcdEhRC6Qi308AQuosQgJLBCSv8S2Q\nignJJEICSwyku1sgDfO9ASF1FiGBJQZS0wjPUZ8FvjcgpM4iJLDEQNL6G4ftc/pi/+sJqbMI\nCSw5kPiCrGWEBEZIdhESGCGhEZJZhARHSGCEZBchgRGSXYQERkhohGQWIcEREhgh2UVIYIRk\nFyGBERIaIZlFSHCEBEZIdhESGCHZRUhghIRGSGYREhwhgRGSXYQERkh2ERIYIaERklmEBEdI\nYIRkFyGBEZJdhARGSGiEZBYhwRESGCHZRUhghGQXIYEREhohmUVIcIQERkh2ERIYIdlFSGCE\nhEZIZhESHCGBEZJdhARGSHYREhghoRGSWYQER0hghGQXIYERkl2EBEZIaIRkFiHBERJYviHV\nR/Hka6qN4gHTVU0RDGmsrY9gikbx55vW1yJPvlxriuTBr6tFHn17SIxtgBESYwYREmMGERJj\nBhESYwYREmMGERJjBhESYwa5CansyhOPv2BBvu+FWTOHz8r3XbBp+vjjzpid7zth0+I/jPnR\nBe/4X+8mpLMnLlxy1ZjafN8NoyrGjkwGpJnjXv3ywVMjOeAp7JpOvX7FyrtGVfnewElIVZd9\nqvrV8PfzfT+Muvy2scmAdOrT+b4HZlUOn69aPtz/dz1OQmpu/jHl+b4LNr08vjYZkEqHP33m\nD8+Zn++7YdN511TV/mu8/4P6nIVUddrt+b4LNlWPm6vJgLRg+G8WV93y40gOzQ69stOHDx/3\nof/1rkJaPOHGKI7MjqBrr9XEQJqn2nDCzHzfD4vqf3l95Yp7x/j/psdRSPNGP5Lvu2DU3HFV\nSYGUGv5B+uPp9+b7flj0+gjvn7JOftj3Bm5CeueEOfm+C1ZdMXL06NEjjr8s3/fDoMZx6T/d\nVh3/Qr7vh0WvDff+8XFcwiGtOvXfqXSJ+OfvKu93cuJTkfyH37C7d8zc1F/GJeJhWTHu+upV\n949c4nsDJyHNG97co/m+H2Yl41s7bbxj7HEXfJrve2HTx5PG/Pj8N/2vdxISY3GLkBgziJAY\nM4iQGDOIkBgziJAYM4iQGDOIkBgziJDC7xI5oOUA2yHDOlk1dLDP3dWP7dF9ceuOvXofcn/G\ntT/qqbq/3x2zABFS+KWf71OaL9hAelTGPLyidce/ufXWKb8bJNdmWutBuqb9cXxz+YCHE7+u\n4XdJt+/1+8q7YAPpH7L6sNBLpPnQoqqtemc6vs2DtFZ/4QMeTvy6ht8l8n63n3gXPEh77OFd\nOmYT1W8e/MK+3QZcUTdxQK9hC9OQdnzt4B79xlWkb33uiN7d97otfWHoNx8ZeGDrXh77Zq9u\nu/y5SYd5384tat1xyzF6v5ZXVi9ds2nTpIFdd713zbd2Tx3Sa/NRH+h30tsOadtX61ZLxg/q\nuvkPEvJfW/MUIYXfJbJykjyv60IaNvCw1xYfJ0dM+uz5Pt9PP6cHDr7igfMKhqvOLDrkkad+\nLlepHr77jje0Hpr7QMF3H5z5azlPF1wif3t1VeuOWyD9Tl5qXdq26Z9kzIx7dh3cCumpgiPv\num3bLZe+f4y8+m7bvlq3OmCLvz0zdbf+iXjXknxFSOF3idSu3H7nuvUgyTzVF+Wg9Kdj0s/3\noXJf+tJo+UT32s57So9If782TP6zeic7DvLwHFtSqrfLi2t23ALp4OLK1qVrNm0asGv6wpKS\nVkj7bFOv+kqX6/QUWWtfzVstkwvSn3542edRfUGSGCGFXxqSPimT14Pk/fzyoZyb/niuVOnQ\nrmlqaSX/+VJ+WZvuZpmtw7rUte7jc/m598tt8uhakKYvXbpk9snyC21Z2rbpJ3KWt+LAFkil\nclrrJh6kdvtq3qpuk61mRnF+sERHSOHnQdJRPT5eF9JW6QuL5PL0x4lSoUO38W54XKbMldb+\no8MGrN7HbPmj98tjcstakJorPm2ltixt2/SVlvUjWyC9Jb9v3cSD1G5fLQNe2kY2GTk1kvMG\nJjZCCr9mSJ/1HqH7dArpG94Nj8mtc+XkWc2lWpY096pM8n6ZLn9bC9I1jz/+xIvev0+0LG3b\n9H8tWI5tgfSOXNy6iQep3b5aBzQ8fe7Osk9NmF+EpEdI4dcMSf8sDx2YhrSX97OL7t8BpJ7e\nt1d/l4fL5Cert2yDtFQmeL/cIk908DPSmqVtmy6U071f9myBVCU/9T77+KtmSO321TZAb5R/\nmP2ON8AIKfxaINXvttVBaUiHb9qU/lmmeweQ5Kn0pWMLl+p+fb2/Y+64sL7983zXAd5Ovttj\nWWeQ2jZt3PQbaZcLClr/sWG3zapU56e/wRsv9e331bzVnB99qd5Pa1eG/HVIdIQUfi2Q9KWC\ngjSka+XyL14/bJf1IR04cIebZk6UE1SfK9n9jid/V3JS+7+RdHrhkQ89/gtvdWeQ2ja9SH5w\n/01bDWmF9Gjh/v+6ZYf+S/VimXRfu301b7W09+63zbj7oD7A2yGydSOk8GuFpD+VNKRVv/5a\n1z0eOb33epD2PnDOwd37ja9Of/rit3uX7HBF/VqQ9KmDe3bd6+/aOaS2TRsu2KLLbg+c0aX1\nBdnpB/Tof9z7qov3Khncbl8tW71xXP+SAce9HubXIPEREmMGERJjBhESYwYREmMGERJjBhES\nYwYREmMGERJjBhESYwYREmMGERJjBhESYwb9PyEhX48uLmh4AAAAAElFTkSuQmCC"
          },
          "metadata": {
            "image/png": {
              "width": 420,
              "height": 420
            }
          }
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "train_fit <- regsubsets(y ~ ., data = train)\n",
        "#summary(train_fit)\n",
        "train_fit_summary <- summary(train_fit)\n",
        "\n",
        "train_plot_data <- data.frame(num_pred <- 1:length(train_fit_summary$rss),\n",
        "                        mse <- train_fit_summary$rss / 100)\n",
        "\n",
        "ggplot(train_plot_data, aes(x = num_pred, y = mse)) +\n",
        "  geom_line() +\n",
        "    geom_point() +\n",
        "    xlab(\"Number of Predictors\") +\n",
        "    ylab(\"MSE\") +\n",
        "    theme_minimal()\n",
        "\n",
        "test_mat = model.matrix(y~., data = test)\n",
        "test_mse = rep(NA,length(train_fit_summary$rss))\n",
        "\n",
        "for(i in 1:length(train_fit_summary$rss)){\n",
        "\n",
        "    # extract the coefficients\n",
        "    coefi = coef(train_fit, i)\n",
        "\n",
        "    # make predictions using matrix multiplication of the test matirx and the coefficients vector\n",
        "    pred = test_mat[,names(coefi)]%*%coefi\n",
        "\n",
        "    # Calculate the MSE\n",
        "    test_mse[i] = mean((test$y-pred)^2)\n",
        "}\n",
        "\n",
        "test_plot_data <- data.frame(num_pred = 1:length(test_mse),\n",
        "                             mse = test_mse)\n",
        "\n",
        "# Plot the test MSE data\n",
        "ggplot(test_plot_data, aes(x = num_pred, y = mse)) +\n",
        "  geom_line() +\n",
        "    geom_point() +\n",
        "      xlab(\"Number of Predictors\") +\n",
        "        ylab(\"Test MSE\") +\n",
        "          theme_minimal()"
      ],
      "id": "a1dbdf04"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "24f9d2f1"
      },
      "source": [
        "c) Identify the best model according the test MSE. How does this model compare to the actual model used to simulate the dataset?\n",
        "\n",
        "*Note: If the model identified contains only an intercept or all of the features, then generate a new dataset (i.e. repeat part a with a different random seed) until the test set MSE is minimized for an intermediate model size.*"
      ],
      "id": "24f9d2f1"
    },
    {
      "cell_type": "code",
      "execution_count": 98,
      "metadata": {
        "id": "87db7d52",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 69
        },
        "outputId": "80ab06cd-918f-4283-980f-7085cc7aa03b"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "8"
            ],
            "text/markdown": "8",
            "text/latex": "8",
            "text/plain": [
              "[1] 8"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<style>\n",
              ".dl-inline {width: auto; margin:0; padding: 0}\n",
              ".dl-inline>dt, .dl-inline>dd {float: none; width: auto; display: inline-block}\n",
              ".dl-inline>dt::after {content: \":\\0020\"; padding-right: .5ex}\n",
              ".dl-inline>dt:not(:first-of-type) {padding-left: .5ex}\n",
              "</style><dl class=dl-inline><dt>(Intercept)</dt><dd>-0.13843287331594</dd><dt>x1</dt><dd>-1.11919382657127</dd><dt>x2</dt><dd>0.583095105056398</dd><dt>x4</dt><dd>0.961570765179235</dd><dt>x5</dt><dd>-1.50350325450737</dd><dt>x14</dt><dd>-1.83917534110568</dd><dt>x16</dt><dd>-0.631422802256244</dd><dt>x19</dt><dd>-1.62920411107757</dd><dt>x20</dt><dd>-0.733443640374377</dd></dl>\n"
            ],
            "text/markdown": "(Intercept)\n:   -0.13843287331594x1\n:   -1.11919382657127x2\n:   0.583095105056398x4\n:   0.961570765179235x5\n:   -1.50350325450737x14\n:   -1.83917534110568x16\n:   -0.631422802256244x19\n:   -1.62920411107757x20\n:   -0.733443640374377\n\n",
            "text/latex": "\\begin{description*}\n\\item[(Intercept)] -0.13843287331594\n\\item[x1] -1.11919382657127\n\\item[x2] 0.583095105056398\n\\item[x4] 0.961570765179235\n\\item[x5] -1.50350325450737\n\\item[x14] -1.83917534110568\n\\item[x16] -0.631422802256244\n\\item[x19] -1.62920411107757\n\\item[x20] -0.733443640374377\n\\end{description*}\n",
            "text/plain": [
              "(Intercept)          x1          x2          x4          x5         x14 \n",
              " -0.1384329  -1.1191938   0.5830951   0.9615708  -1.5035033  -1.8391753 \n",
              "        x16         x19         x20 \n",
              " -0.6314228  -1.6292041  -0.7334436 "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "which.min(test_mse)\n",
        "coef(train_fit, 8)\n",
        "# -1.1144403712824 0.52614872250706 0 1.03203241433948 -1.46307313907892 0.51514502428472 0 0 0 0 0 -0.527756507508457 0 -1.89949801936746 0 -0.712022013962269 0 0 -1.55381905566901 -0.679967186413705\n",
        "# x1, x2, x4, x5, x6, x12, x14, x16, x19, x20\n",
        "\n",
        "the best model accoridng to test MSE missed two features from the true model that was used to simulate the dataset (x6 which is meant to be -1.46307313907892, and -0.527756507508457)."
      ],
      "id": "87db7d52"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "894b4503"
      },
      "source": [
        "> *Write your response here*\n",
        ">"
      ],
      "id": "894b4503"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "35b974f0"
      },
      "source": [
        "**DUE:** 5pm EST, April 8, 2024"
      ],
      "id": "35b974f0"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "957a72c5"
      },
      "source": [
        "**IMPORTANT** Did you collaborate with anyone on this assignment? If so, list their names here.\n",
        "> *Someone's Name*"
      ],
      "id": "957a72c5"
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "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": "4.2.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}