08-control-statements-iteration-apply.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Control statements, iteration, and apply functions\n",
        "\n",
        "Joshua P. French  \n",
        "November 3, 2023\n",
        "\n",
        "To open this information in an interactive Colab notebook, click the Open in Colab graphic below.\n",
        "\n",
        "<a href=\"https://colab.research.google.com/github/jfrench/DataWrangleViz/blob/master/08-control-statements-iteration-apply.ipynb\"> <img src=\"https://colab.research.google.com/assets/colab-badge.svg\"> </a>\n",
        "\n",
        "------------------------------------------------------------------------\n",
        "\n",
        "In this module we will discuss control statements, iteration, and apply functions. These are fundamental tools used in R programming. A simple analysis will not need them, but as we analyze more complex data, these tools will be critical in completing our analyses.\n",
        "\n",
        "We start by making sure the necessary R packages are installed."
      ],
      "id": "79c14ba3-1628-423c-ba01-b02e37943871"
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [],
      "source": [
        "if(!require(microbenchmark, quietly = TRUE)) {\n",
        "  install.packages(\"microbenchmark\",\n",
        "                   repos = \"https://cran.rstudio.com/\")\n",
        "}"
      ],
      "id": "6c89d039"
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [],
      "source": [
        "if(!require(pbapply, quietly = TRUE)) {\n",
        "  install.packages(\"pbapply\",\n",
        "                   repos = \"https://cran.rstudio.com/\")\n",
        "}"
      ],
      "id": "845e2f17"
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {},
      "outputs": [],
      "source": [
        "if(!require(future.apply, quietly = TRUE)) {\n",
        "  install.packages(\"future.apply\",\n",
        "                   repos = \"https://cran.rstudio.com/\")\n",
        "}"
      ],
      "id": "332d7364"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Control statements\n",
        "\n",
        "A control statement is code that controls what is executed later in the program or script. The most well-known examples are “if-else” and “switch” statements.\n",
        "\n",
        "## `if`, `else if` and `else` ([YouTube](https://youtu.be/mHaYfe-bXME), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=f7efc57b-07f1-4f7c-8d4a-b0ad01017003))\n",
        "\n",
        "If-else statements in R rely on 3 functions: `if`, `elseif`, `else`.\n",
        "\n",
        "-   The `if` statement can be used by itself.\n",
        "-   The `else` statement must be used in conjunction with the `if` function.\n",
        "-   `else if` must be used in conjunction with the `if` statement and should be used with an `else` statement (otherwise there may be cases that are left out).\n",
        "\n",
        "An `if` or `else if` statement is used to perform a sequence of commands when its argument is `TRUE`.\n",
        "\n",
        "An `else` statement is used to perform a sequence of commands when all previous control statement arguments are `FALSE` in a group of control statements.\n",
        "\n",
        "A set of `if`, `else if`, and `else` statements will having the following structure:\n",
        "\n",
        "    if (logical statement) {\n",
        "      # code to run if logical statement is true\n",
        "    } else if (another logical statement) {\n",
        "      # code to run if another logical statement is true\n",
        "    } else {\n",
        "      # code to run if none of the other logical statements were true\n",
        "    }\n",
        "\n",
        "The curly braces, `{}`, are used to indicate the set of commands we want to execute if the logical statement is true.\n",
        "\n",
        "The `else if` and `else` statements should be one space after the `}` used to indicate the end if the previous `if` or `elseif` statement.\n",
        "\n",
        "It is common to indent the code within a control statement by one tab to make the code easier to read.\n",
        "\n",
        "### Categorization example\n",
        "\n",
        "According to the CareGivers of America at <https://caregiversofamerica.com/2022-generation-names-explained/>, here is a breakdown of birth years for seven recent generations.\n",
        "\n",
        "-   The Greatest Generation: born 1901-1927.\n",
        "-   The Silent Generation: born 1928-1945.\n",
        "-   The Baby Boomer Generation: born 1946-1964.\n",
        "-   Generation X: born 1965-1980.\n",
        "-   Millennials: born 1981-1996.\n",
        "-   Generation Z: born 1997-2012.\n",
        "-   Gen Alpha: born 2013 – 2025.\n",
        "\n",
        "Let’s say we have a vector with the birth year of an individual and we wish to identify the individual’s generation based on the dates above. We can do this using a sequence of control statements. For simplicity, we only use the last four generations (Millennials, Generation Z, Gen Alpha, and Other).\n",
        "\n",
        "We create a sequence of control statements to categorize the generation of a single number stored in the `year` object.\n",
        "\n",
        "We set `year` to `1950`."
      ],
      "id": "678b90ca-bb94-42b5-ba07-225595263268"
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [],
      "source": [
        "year <- 1950"
      ],
      "id": "9db61da9"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we create a sequence of control statements to categorize the `year` variable"
      ],
      "id": "8ffcbf99-ed76-4896-8a11-a42d26861f51"
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [],
      "source": [
        "if (year >= 1981 & year <= 1996) {\n",
        "  generation <- \"millennial\"\n",
        "} else if (year >= 1997 & year <= 2012) {\n",
        "  generation <- \"generation z\"  \n",
        "} else if (year >= 2013 & year <= 2025) {\n",
        "  generation <- \"gen alpha\"\n",
        "} else {\n",
        "  generation <- \"other\"\n",
        "}"
      ],
      "id": "8d963df8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now look at the result stored in `generation`."
      ],
      "id": "a4c5acf1-8780-49cf-b2c9-55effedfc7dd"
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [],
      "source": [
        "generation"
      ],
      "id": "603c46d9"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As expected, `generation` categorizes a person born in 1950 as being part of the `\"other\"` generation.\n",
        "\n",
        "It would be even better if we made a function to do this categorization so that we can do the categorization for many different values of `year`.\n",
        "\n",
        "We create a function, `determine_generation`, that takes the argument `year` and returns the generation corresponding to that year."
      ],
      "id": "1e8658eb-90d5-44c4-bc8f-c884a8cc0888"
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {},
      "outputs": [],
      "source": [
        "determine_generation <- function(year) {\n",
        "  if (year >= 1981 & year <= 1996) {\n",
        "    generation <- \"millennial\"\n",
        "  } else if (year >= 1997 & year <= 2012) {\n",
        "    generation <- \"generation z\"  \n",
        "  } else if (year >= 2013 & year <= 2025) {\n",
        "    generation <- \"gen alpha\"\n",
        "  } else {\n",
        "    generation <- \"other\"\n",
        "  }\n",
        "  return(generation)\n",
        "}"
      ],
      "id": "f046228f"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now use the function to determine the generation associated with the value 1985."
      ],
      "id": "28067ffe-2417-4d5d-82b0-460c8a7d49c0"
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {},
      "outputs": [],
      "source": [
        "determine_generation(1985)"
      ],
      "id": "5cdbba03"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The above example is trivial and unrealistic. E.g., our control statements won’t work with a vector of inputs. Additionally, it doesn’t handle cases like non-integer values very well (what generation is 1981.5?). However, it paints a picture of what can be done using control statements.\n",
        "\n",
        "## The `switch` function ([YouTube](https://youtu.be/ay5dwOMNlZ0), [Panopto]())\n",
        "\n",
        "The `switch` function is a compact way of evaluating a set of `if`, `elseif`, and `else` statements when we need to change our action depending on a single value argument that can only take a small number of values.\n",
        "\n",
        "The `switch` function takes the following arguments:\n",
        "\n",
        "-   `EXPR`: a single value that evaluates to a number or character string.\n",
        "-   … A list of what to do based on the value of `EXPR`.\n",
        "\n",
        "### Animal example\n",
        "\n",
        "Let’s create a function called `sound` that takes a single character value `animal` and returns the “sound” associated with that animal.\n",
        "\n",
        "The `switch` function takes the value in `animal` and then decides what to do."
      ],
      "id": "5fae1b5d-fd3f-4818-85be-9bc26dcf517c"
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {},
      "outputs": [],
      "source": [
        "sound <- function(animal) {\n",
        "  switch(animal, \n",
        "         \"dog\" = \"ruff\",\n",
        "         \"cat\" = \"meow\",\n",
        "         \"horse\" = \"neigh\",\n",
        "         \"frog\" = \"ribbit\",\n",
        "         \"cow\" = \"moo\",\n",
        "         \"duck\" = \"quack\",\n",
        "         \"rooster\" = \"keekarikee\",\n",
        "         \"coyote\" = \"yipp\",\n",
        "         \"donkey\" = \"hee-haw\",\n",
        "         \"we don't know what sound this animal makes\")\n",
        "}"
      ],
      "id": "302e4977"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "-   If `animal` is `\"dog\"`, then the `switch` function returns `\"ruff\"`.\n",
        "-   If `animal` is `\"cat\"`, then the `switch` function returns `\"meow\"`.\n",
        "-   The last input to the `switch` argument in the above example is what we want to return if the value of `animal` doesn’t match any of the previous choices.\n",
        "\n",
        "Let’s use this function."
      ],
      "id": "78c97d3f-4049-4d11-b517-7d6dbce96d70"
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {},
      "outputs": [],
      "source": [
        "sound(\"dog\")"
      ],
      "id": "4d0e0a1b"
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [],
      "source": [
        "sound(\"duck\")"
      ],
      "id": "c1aaec2e"
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {},
      "outputs": [],
      "source": [
        "sound(\"goose\")"
      ],
      "id": "c9ef46aa"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Iteration ([YouTube](https://youtu.be/xX6z0xqkBBc), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=97161e03-c736-4dde-a53a-b0ad0103fc2e))\n",
        "\n",
        "Iteration in programming is when we repeatedly (i.e., iteratively) run the same code block until some condition is satisfied.\n",
        "\n",
        "The two *main* iteration approaches in R are `for` and `while` loops.\n",
        "\n",
        "There is also a `repeat` loop that repeatedly runs the same set of code until the `break` function is used.\n",
        "\n",
        "-   The `break` function is used to “break out” of the current iteration of the inner-most loop.\n",
        "-   The `next` function is used stop the current iteration and advance the looping index.\n",
        "\n",
        "I have never once used the `repeat`, `break`, or `next` functions in over 20 years of R programming. While I’m sure they have their uses, they are not commonly used nor needed, so we will not be covering them.\n",
        "\n",
        "## `for` loops\n",
        "\n",
        "### Basic `for` loops ([YouTube](https://youtu.be/7RV5Aayi2xc), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=94a21d98-fd8d-44bc-861f-b0ad0104afba))\n",
        "\n",
        "`for` loops are useful if we know in advance the set of values that we want to iterate over. E.g., we may want to perform an action for integer values running from 1 to 100.\n",
        "\n",
        "The basic syntax of a `for` loop is shown in the psuedo code below.\n",
        "\n",
        "    for (i in vector_of_values) {\n",
        "      # some code to run\n",
        "    }\n",
        "\n",
        "We have the following components of the `for` loop:\n",
        "\n",
        "-   The loop starts with `for`, which indicates that we want to perform a `for` loop.\n",
        "-   Inside parentheses after `for`, we have the iterator name, which is the name of the variable that our loop is changing. In this example, we called the iterator `i`. Names like `i`, `j`, `k`, and `l` are commonly used, but we can use any name that won’t conflict with other object in memory.\n",
        "-   After the iterator name, we use the syntax “`in`”, which indicates that our iterator will take the values in `vector_of_values`.\n",
        "-   `vector_of_values` is a vector containing the values the iterator can take as we progress through the loop.\n",
        "-   After the parenthesis, we use the curly braces `{}` to indicate the code we want to run in our loop.\n",
        "-   The `#some code to run` is the code we wish to execute in our loop.\n",
        "\n",
        "Here’s what happens in a `for` loop:\n",
        "\n",
        "1.  The iterator takes the first value in `vector_of_values`.\n",
        "2.  Using that iterator value, the loop runs the code in the curly braces.\n",
        "3.  After that code runs, the `for` function will change the iterator to the next value in `vector_of_values`.\n",
        "4.  We repeat steps 2 and 3 until there are no more values in `vector_of_values`.\n",
        "5.  We exit the loop.\n",
        "\n",
        "Here’s an example of a simple `for` loop that prints the current iterator value for the values 1, 2, and 3."
      ],
      "id": "e4089728-271c-403c-a27e-7e03873600f2"
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {},
      "outputs": [],
      "source": [
        "for (i in 1:3) {\n",
        "  message(\"The current iterator value is \", i)\n",
        "}"
      ],
      "id": "0aa2b980"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `message` function prints a diagnostic message to the screen. Notice that the `message` function can combine text and numeric values into a single message. We have to manually include spaces in the text.\n",
        "\n",
        "Here is a loop that takes a character string of animal names and returns the sound made by the animal. We use the `sound` function we previously created."
      ],
      "id": "ddd77e93-7b3b-415b-bcca-cd42638d6839"
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {},
      "outputs": [],
      "source": [
        "for (animal in c(\"dog\", \"cat\", \"coyote\")) {\n",
        "  current_sound <- sound(animal)\n",
        "  message(\"The \", animal, \" makes the sound \", current_sound, \".\")\n",
        "}"
      ],
      "id": "091fccc7"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Q-q plot example ([YouTube](https://youtu.be/0sfB22OUMMM), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=a4fb684c-349d-4e11-8ff9-b0ad01054150))\n",
        "\n",
        "Let’s create a more realistic `for` loop example involving quantile-quantile plots (q-q plots).\n",
        "\n",
        "A q-q plot is a way to compare the quantiles observed for a data set with the theoretical quantiles we expect if the data comes from a certain distribution.\n",
        "\n",
        "-   We plot the pairs of quantiles in a scatter plot.\n",
        "-   If the points of the q-q plot fall close to a straight line, then it is plausible that the observed data are a sample from the theoretical distribution.\n",
        "\n",
        "Here’s an example of what a q-q plot might look like if the observed data are a random sample from the theoretical distribution being considered."
      ],
      "id": "e7b859f8-ecae-4ba1-95e4-5ee525f7292c"
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyMjIya\nmpqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAASdAHeZh94\nAAAgAElEQVR4nO3dh4KjOhZFUeFUzvz/35bJIhrMVd5r5nVX2WAE0mlACKxyALsp1wUAYkCQAAEE\nCRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJ\nEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQ\nQJAAAQQJEECQAAEECRBAkAABBAkQkHqQlJraAuddnzkx9/Rimknn3t28zK8fND1B/eqqbbE47eby\nxCOZFZ0xVdWPbM9WmZx7ukU1k+5ub6s/aHOQRmuzKkhiKxaOZFZ0xlRV76v+ybkXG/D+9rb6gzYH\nafTaqiCJrVg4klnRGXEF6ccJF2afe2tdkBKS3Ao3XudMHa5tld9On58Ol1fZCFT9avdiq3jnelCH\nR55fM3V85PNzvy8HpY73drbHUanzq/dR9aTDd9+XTGXtYu/nz/vnezvT66AuvYn6HzRcdG8tBi28\nvxGacv99iqJOt+EnV8vVpv1sh6wucu+fhHF5RuuQ34vVfa6trACkGqRHVdnHuqqPdeWrh9YMtBdb\n9TyfZnXp3pqc+5VVP17q2arps5f+UV17099t5qwW23z4qZmpyEhvolHD7S26txb9IA02Qr/cxUJ6\nn3xoXupth6qUX4I0WodqdVVESUo1SE1rqar6+mkW7zy/FFXdvqq/2GrnyrqmMTl3t4S7Pps6Dz9q\n/G4zZ1b8cmrfO3Uz3XoTjRquvuj+WqhekPoboX7zXH78+9P2r6Mi3vQg6aVUi0GaXIf+xgheokG6\nfRrAp5Xds6qqP//Yakc+dfX3X6ypsl1+Wqc6PMu/5ub+vJc9y/Z4qF7MHmVzzvLehzV/a+9WTf9d\nNej7573r57e/LpFFLvoTDU62eoueLNzkRugmKWZ4twXX1jvXg9TMOy7A4NeJdSjmPfdTHbiIVmWL\nU1WlZR1rL0+0tnGQHr2/5uaul/A+/NWnMeVvUzMM3/3MWTfZU9naruU0l+rf73rS/kSDhttf9Pyq\nDTeCavekzdnMqIh6VJp5xwUY/DqzDm+CFL62DtsfXrfLUQ1bm/Zif8ZBg5mYezYyC6/qR1H1UZOq\nA5O/2jebHYN2aDX6gN7KTq7axEao/vqrj7p6uenKMVy9/udOBmliHSY2RuAiWpUthm3odmjapfZm\n78X+jP2/puaWCJLSZxg1wokSTwVpctUmN0L916VJ6Guq4L8FaW4d8mhEtCpbDNrB52xBHc7XZ6+i\n+y/2Z+z9NTn370HKBstr/zXX9z3Zwgf1FjJZuMmN0P7+vlW9bMepgnfT6nuob0GaXAeCFL7m9OBW\nVeah/rXXKvov1qaCNDn3cXiOlI8+ayZITdmagg7PL3orMPVBvUVPFm5yI/TevJ9n9h7dtLdyuvYc\nqcjKYzJIc+tAkMJ3rTqOblmvuWj/bL+HL9amgjQ596jXTp9d+7D36N1b1Yd3K3cJEz1e5aS9iQYf\n1Fv0uHBzG6H7Z6E+C8smiqitedmp+Jn3Ly93kJdyiF33Sd1cc+tAkCLQv4RyLNtB0w9c/HUZvlib\nCtLk3N0Srvlc22kmHbzbzvnIteup1UWX9gN6Ew0+SF90v3Dj5XcboWv2x1fZ53CZKKIepEpWJObc\n/6TBXDPrQJAiUF/UP1WV+WgaRdkui1ZxHL5YmwrS5NzVP8+qG9mgz15rJh28e1fanG0rHNyY0Jto\n8EH6ovuF6y+/vxGaN5vOhuNUEbU1ry+zlkeHr3qB9bvDuabXgSDF4PWp7GM7zOz5+S07P1/VEf+p\nqvL+i5WpIE3PXQ6GU6dRL7JeiNP0SUM5VO7UnAPdz1l/nFo+nmjwQfqie4UbLL+/EZo3y/Oj43Xy\nk/Vpr9rwwGIpx5t+Cteba3IdCBKAPoIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAA\nggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACC\nBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgAALQVJAYH5o5fLBcbAI4Aez\nLZMgAavNN0yCBKy10C4JErDSUrMkSMA6i62SIAGrLDdKggSs8aVNEiRghW9NkiAB331tkQQJ+Op7\ngyRIwDcr2iNBAr5Y0xwJErBsVWskSMCidY2RIAFLVrZFggQsWNsUCRIwb3VLJEjArPUNkSABc4bt\ncOGOcoIEzBg0wzJFc1EiSMC00f5o8tWZiX/4fAMIEpybztFc4yRIwJRRGyRIwGbjJkiQgK2mWiDn\nSMA2kw2QXjtgk7n2x3UkYL1fmh9BAvp+an0ECej5rfERJED3Y9sjSIDm16ZHkIDOzy2PIAGt3xse\nQQIaO9odQQJqe5odQQIqu1odQQJK+xodQQIKO9scQQLy/U2OIAECLY4gAQINjiABAu2NICF5Es2N\nICF1Iq2NICFxMo2NICFV1RMYhNoaQUKa6mcCSTU1goQ0LT6l7tePMz2Lh4tA2pocsUcCdlD1HwQJ\n2EFp/0l9nvlZPFwEEqdEc0SQkKjFJ3n/8HFWZvFwEUicEuz7zgkSEiXdxAgSUiTewggSEiTfwAgS\n0mOgfREkJMdE8yJISI2R1kWQkBgzjYsgIS2G2hZBQlJMNS2ChJQYa1kECQkx17AIEtJhsF0RJCTD\nZLOyGaT3Wanjvf6QxU8hSJBntFVZDNI7U4VT9SEECXaZbVQWg3RR10+artmx/BCCBCuau44MtymL\nQcqqGV/Z4UWQYEd7H6zpJmUxSE123sfjVJCU7sdFAH3ST936siDTs5QO6t38dGSPBBvU4G/zSzI7\nS+mqzvVPL3UkSLCgbkgWDnFsdn9f2vW5fzl6I0gQYe3Izu4F2eep+el1JkiwoHp2nZ9jCPwsFTBB\nqfK5WxYWZGUWDxeBRFjqAyZIiJqtpkSQEDNrLYkgIWL2GhJBQrwstiOChGjZbEYECbGy2ooIEiJl\ntxERJMTJchsiSIiS7SZEkBAj6y2IICFC9hsQQUJ8HLQfgoTouGg+BAmxcdJ6CBIi46bxECTExVHb\nIUiIiqumQ5AQE2cthyAhIu4aDkFCPBy2G4KEaLhsNgQJsXDaaggSIuG20RAkxMFxmyFIiILrJkOQ\nEAPnLYYgIQLuGwxBQvg8aC8ECcHzobkQJITOi9ZCkBA4PxoLQULYPGkrBAlB86WpECSEzJuWQpAQ\nMH8aCkFCuDxqJwQJwfKpmRAkhMqrVkKQECi/GglBQpim24hSPNfO7iIQtqaJ9JJT/uwmSgQJIapb\nyCA5Sn/PSYEMz+LhIhAy1ftbTb5qF0FCeGaSQ5AcLALhUsMfCJK7RSBYavTT9JGeVQQJgVGjn7s9\nFL121heBQPXaxig5XEeyvQiEadg0nCVngCAhJGXL8CU8OoKEgLRHct5FiSAhHKr9w7smQpAQDC1H\n3rURgoRQTI8H8gRBQiA8GL6wgCAhDKrpq+McyatFICxdXx29dl4tAkHp7428ixFBQhCUp2dGHYIE\nr5V7H+VrF0OHIMFj7fkQQfqJt5sLdtVD6/Lc0766DkGCv5qhDMrXvroOQYK/lPafl311HYIEf6lu\np+Q7ggSPqVByRJDgM9/PjDoECf5Svp8ZdQgSvFQEKKRmQJDgoXAO6RoECR6q+rxDagYECf7pXT8K\nA0GCfwK6ftQgSPBPc1QXUDsgSPBQcDkiSPCQotdORkhbEOKqFIXVCAgSfBNk7RMkeCbMyidI8Eug\ndU+Q4JVQq54gwSfB1jxBgkfCrXiCBH8EXO8ECd4IudoJEnwRdK0TJHgi7EonSPBDXeehDQ1qECR4\noR7vHdxg1QZBgg+U/leI1U+Q4IFejoKsf4IE95r7JgiSsAA3JH6nhidHAdY/QYJr3b5IcY4kK8Qt\nid/09kX02skKcUviN9Vzt9reBq4jSQpzW+IHzcFcuAd1FYIEZ5RSvWO6kOudIMGR8svKwz85qhEk\nOKK6fu+QT45qBAkulLHp+hnCDlGBIME+/agu+H1RhSDBvmZIUPNzBAgSrGsO6VQeT1UTJFhXBag6\nTYqlqgkS7FJ6P0M8FU2QYJPSTo+iqmWCBJtCv6N8FkGCRVqPd2R1TJBgUXNIF18FEyRYpLT/4kKQ\nYFNcF480BAk2RdjNUCFIsCmSkXVjBAkWxVuxBAn2RFyvFoOk+kwsAl6LuVotBum6HKTVKUOgoq5V\nm4d2z+xoehHwV9yVavUc6akuphcBX0Vep3Y7G67qaXoR8JN+U2yM6LWDDc09fNFGiSDBgm5cUKxV\nS5Bgnj6eIdK6JUgwTvUfvxUlggTTmhHfBGn3LB4uArZ0CYr0DooSQYJZ7f2wcT1+a4ggwaheP0O0\nMSJIMEtpV4+irlWCBINU8yTIPOKjuhJBgjkqT+GorkSQYIz23K3oa5QgwRS9tzv6GiVIMESN/owZ\nQYIZ9X0TCfQzlAgSjOjum4i+n6FEkGBC/PdNDBAkGJBUP0OJIEFebyBDGpVJkCCu/3yGNCqTIEFa\n//kMidQlQYKoYk/UfEdsIj3fJYIEQVWKyj+rX5OpSYIEQar6nxr0NySAIEFKuSdqe75jvq98AkGC\njOZJQc2JkUqrGgkSZFR7oCZIzRFeMggSRHQ56rKUEoIEEfXBXNP5nViMCBJktGO8tetISSFI2KPO\njWqO5ZKtub1Buh7y/HVQh4dUgcaLgLdUdzSn0jyia+0M0r3YclmxEUWTlGx1BKbuYVDDcaoJ2hmk\no7rlT3XIb2rt18NuXgS81V13jfu53mvsDFLxb1D5zbCy/xilXCMBaYOkcoK0b5ZiE57UnSAlaRCk\npGtt96Hd866ynEO7NLXHdWl3NBT2dzYo9VdsxLtYkXKCFIq22zuRRwUt2N39nRVnSPnhJlSeiUXA\nY6lefx3hgix+1gz2BkHCj9J5huo6u4N0P5U9dy+h8kwtAt7pbj6iqip7g3Ssh4ZkokmidrzVje+u\nd0fUVWlnkK7q+C4251WdxYqUUzneanvoug5v6qq0M0iZept46hKV4yltcB3jgnoERjYQpCQ0+yGG\nM0zaGaRDvUcqBq4KonZ8o116zbUjPHrtGjLnSPdMXcWKlBMkr2jdC93eqLuhD6W9vXan+p4u0aF2\nBMkfg11RXueJEA2IXEdSJ9kRQgTJH210cu2YjlFBI4xswJI6Qc0uSRGjGQQJS1TvvyZWbsvkpR1B\nUn2OSwUjhsd0OeNUpxEkLNK7F9qxQRjh0A6L+r12OXUzgyDhi373AlUzbdehXe/wznGpYMSwYqmZ\nGQQJs8b37lExczi0w4z2O4662qBeZhEkzOh/i2WeUy1LBG6jKGWZRGmmFgE3ml1Rt0uiVhYIBenF\nOVJsRkGiUpbsCNK9dz2W+5Ei0x8YRJ18sWePdNBzxNe6xKYZD6SaHRMWSJ0jyaLWnGufW8f+aBV6\n7TChCRH9dWsRJEwY7Iaoj6/2BunvwMiG+AwO6KiO73YG6Y8hQnFpK1N7+CO1scLOIAk/PWhqEbCn\neTRQ77ssqYw16LVDo3nmVpsicrTeziCd1FusKDOLgBVKG6PaPJRYH2eHZTuD9MqOoldiJxYBC/SR\n3r3n1lETK+0+tKOzIQbdA4Ka/ziu24Ygobf/6Q0Moh5W44Js2gbd3W2vHfujjQhSysbd3frRBbWw\ngVSQHqe9Jfm6CIjTu7u1kd7tm1htb5AunCOFq+ul677wSH8T6+0MUpeju1iRcirRkonubv09bLB7\niNAtP6rX66i4sS887T1HbXd395abEoVLYIjQ32dv9JT9pjGq0YLuIuzw7IgK2E4gSPdi4CrnSMFR\no+5u/S1ss3us3S1/qUP+IEiBaS4VTXUUsfm32xmke1EJx6IuzmJFyqlJ4xafasLW/8He7u+/4rez\nUheh8kwsAvKa/rp8Yluz8X/ByIYU6QPrJt/DVgQpRfXp0aizLmfT/4ogpagZTDfuImLL/4jbKJKk\ntD/Hr2M7gpSk0TeI1S87KEokZA7tHkfRwd9UqHlT//Sx2X8ndI705jpS+NjqO0h1NnBoF4y5I3E2\n+h5CQboqvrEvDPpNsf03nBQnGmKdDX9iRcqpVIOU9v/B69hBKEgH2ScXU6umaAODuG9CEhdk09J/\n8Fb3sqPixIMgpaT3PRPaRmZ777YvSO+/4haKU/G8hpPkSRIVa4bKJ8+R2Nz77QrStT1Fel1Fn35C\nzRpRD2cY9tqxtQXsCdInR5fimSfPi8p4ZkMA6hF2g+tIbGwJO4L07p7BdeNxXCFQg7/Hv+BXO4J0\n6W6LPQnfIkvlmjEx6JtNLWNHkA7qVf/0Un/qIFUifRGQNR70zZYWsiNIXYU8T4y1C8RglB0bWopI\nkIa/7Eb92sF2FiNyaJeXz7YTRAXLmxj0zWaWI9PZ0PtZADUsbWrQN1tZkEz3913Jfrs5VSxtYkAD\nG1nSnguyN6Uuz7y8ICt7GYk6ljYx6JttLGrXEKFbVg8RymRzRCVLGweJTSxr56DV66kYtHoVPa7L\nqWVxoyCxhYVxG0USyn4GcmQQQUqAar6tfP4bKLATQUpAFZ7uOhKbV57VID3+TmVtni5fvnGWmpY0\nzA9b1wCLQXoftAccL9++RFXL6d1enudsXDMsBumistuz/Ol1z5ZHQlDXUupehrwLEtvWCItBytSz\n/fm5/EBJKltG/RVIVYrIkUkWg/RltLjS/bgI9NU93kq1J0psWUN2B+l+Kmro9JqZWMMeybbuuI7+\nOtP2BulY7T9U9j1Jn3OkezUV50gWqG5HxHVY83YG6aqO76Kqrmu+1uWoHbsdFkcVUeE/646Oq6uv\n3U9sVpN2BilT77lvf5vwuJTXkbLTH9eRzBiMYeC4zpqdQequUYjWElX+o+Y4rglS18vARjVrZ5AO\n9R7pya3mPujnSDUnSs17MEfmHOmeKdHvdaHSf6MHqT62I0d27O21O60a8rNrEVitDVK7T+K4zhKR\n60jqdBMqzuQisJLW162NaSjfcVquFHAbRSyay0Zdr123FdmexhGkOOjXjUajrNic5u0IkupzXKrE\nDcd499+DcQQpfFoFqIltx8a0gUO70PWeoToRJLalFQQpbHVHnTYeaDiBi1IlaHeQbkX395kHRDrR\nXSoaXjdqp3BSrgSJ3EZRXEmSKtB4EZij2qfVtR12wylcFCtJO4N0qZ5WzBAhB9qnmtS7pIn+Hraj\nNbtvo6juemXQqn2q2xVNd9exGS0SuI2i/4MIWsB3TXSmvvmonQK27D60a/ZIoidJNIHvuguwM9fx\n2Ig27e1s+CvPkR4Zo79t07oYZt+HNbsP7YyMbqARrFAf2M1sKzahXQQpWIu3+LMFLWNkQ8Dm/+li\nA9pGkGLE9rOOIEWIzWff3iBdMunzo9EiMInjOq/svo7E/UhOLPQ0sPFc2N1rJzrGbmoRmKK0P6fe\ngWVSQ4Rk0Ri+UIO/x+/Art2HdosPw/8VreGL2SCx5RzZfT/SccU3I+1bBMbmgsSGc2VvkO50Ntim\nbW5y5I2dQfqj184y7YtbuK/cI7tv7KPXzq7uPj7uK/cJvXZhaW/n47jOL7sP7ei1s2o+SGwzp3bf\n2Hf88i2WP6FRzJkNEpvMLbn7kcSKlNMqFijt//2X4RJBCsz0s07YYK5xG0VQlJp81gnbyzmCFJC5\nId9sLvekgvTgcVzmzQz5Zmt5YPeNfZwjWTMzMIiN5QOxG/tEv4+CtjHWDQri/MhDu4cI3fKjer2O\nSvRyEo1jbObLLdlUfhAYIvT32Rs9leijVmkdQ2rmq8TYUp4QCNK9GLjKOZJJ9XeIjZ9QzIbyxc4g\nnT6Hdi91yB8EyaS5Id9sJ2/sDNK9qNjyW/vOYkXKaSB9SvueZXLkqd2DVovfzkpdhMozsYjENV8N\nO/6KWLaSRxjZ4LtuT8RxnccIkue0mya4fuSxPUF6X8ofHweVCd9wTiMptUNGxod1bCLP7AlSVlZt\n9RwhvrFPXO+GCcZ7e25HkK7qWNxnnmXP/H1UN8elik/T4T0ezcAG8s+OIB1V8WzIh/or/2Rkg7Dm\npGjisI7t458dQaoq91KNsuOCrLQ2SOMHb7F5/LM7SAel/SKFlqIHieO6AOwI0qE4tHtVQxreKhMs\nFE2loLT/D16Hd3YE6VJ0NpyrG5GuDBESN/2YE7aNn3YE6Z21/d5XpZ6ChaKx5HOPOWHTeGrXBdlm\njJ1irJ2w2W+2TH7L+EpkiJA6CT9uNfnmorQ/x6/DP4y185Ea/D18Hd4hSJ6pz43qXwbv2S8OViJI\nXqlvJm+66uivCwZB8kcz1Lv+HzkKCUHyRfNkk2asN/dNBIUg+aLbGdWjvtkfhYQg+aF7wElzTEeO\ngkKQ3FPdkRznR6EiSE6p9qmPqj2g6x7y3U3mqnxYiyA51D2PIVftDbFqfANSIpsjaATJoV7/wtSd\nsN1k8BxBcqe9b2+mo06fDJ4jSO50N8BO3sCnTwXfESR3miHek/0L/YngO4LkULVLqs+Q2B8FjSA5\n1GZodn0T2RARIEiuTHxt2HgaW4XBXgTJjdl7yXsTWSkKJBAku7SBDPm39Yx3K0SIINnU7ofWJCnW\njRAngmRT203XfDf5worGug0iRZAsGl6BXVrRSDdBtAiSNUp1QZq8V6I3sa1SQQZBsqQ9M+rul1jo\ntYtwA0SOIFmij2Koflzo/I5w/WNHkOyo+xba/dDiJaT4Vj9+BMmOZo0YzxApgmSHGvy9YlKEhCBZ\nsmYow9pJ4B+CZMmqwXV5lKueBIJkzZczo3oi8+WACQTJK8muePAIkgWr9kXlhGbLAXMIknFrz44i\nW+3EECSzVt56tHYS+IogmbT21qN1E8BjBMmkaoT3uiDFss6JIkgGabdNfF2nSFY5WQTJoGZXtHzr\nkTYpgkWQDBoO+f42JcJFkExqk8T+KHYEyaSVl5BiWd2UESSz1gxqiGdtE0aQnEtqZaNFkFxLaV0j\nRpAcS2hVo0aQDFk54juCNUWBIBnB/bCpIUhGME41NQTJhJXPDAp9NdEhSCasC1LoawkNQTJhVZBC\nX0noCJIRK86Rgl9H6CwGSfWZWIQ3eL53aiwG6ZpQkL5eR4pgDaGzeWj3zI6mFxGK6FcwOVbPkZ7q\nYnoRYYh9/RJkt7Phqp6mFxGCyFcvSf702q0+gQpe3GuXKH+CZHkR7kS9cskiSOLor0uRiyBF/RCD\nb1eQAl41LCBIwr6MaQh4zbCEIMn6Msou3BXDMoIkazlI4a4XviBIshaDFO5q4RuCJGzhHCngtcI3\ndH8Lm++1C3il8BVBEjdzHSnodcI3BMmSCFcJGoIkZ2lIQ5hrhNUIkpTFIQ0hrhC2IEhSloY0hLg+\n2IQgCVm6gBTg6mAjgiRkIUgBrg22IkhC5oMU4MpgM4IkZe4cKcR1wWYEScpMr12Iq4LtCJKcqetI\nYa4JNiNIRkWzIviCIJkUy3rgK4JkUCSrgRUI0n5zY+zCWgvsQpD2mh1jF9JKYC+CtBfXj5ATpN3m\nRjQEtAoQQJB+V54bzQQpjDWAGIL0q/rcaDpIIawAJBGkXzXnRlPnSCGUH6II0lb1Tkj1fu332vlc\nfJhBkLZpY6Md0g2vI/lbehhDkLZpD+S4/wg6grSJvh/qvTCeBEkhSJv0Duhy7j9CgyBt0jug4/4j\ntAjSNnMHdL23kR6CtM3yN1v6W24YRpC2Wngwsc/FhlkESU6YpYYIgiQmyEJDCEGSEmKZIYYgrbP0\nlS3VBHbKAU8RpDWW++rKKWwVBX4iSGt8uXrkX4FhG0FaYX6A6rc3kAqCtMK3IHlWXDhAkFb4EiTP\nSgsXCNIai+dIvhUWLhCkNZZ67XwrK5wgSOvMXkfyr6hwgSDtE05JYRRB2iWYgsIwgrRHKOWEcQRp\nHuPrsBpBmsP4OmxAkOYwvg4bEKQvRWB8HdYgSF+KwPg6rEGQvhSB67BYgyDpi9X7Fhhfhw0IUrfQ\nfjcd4+uwAUEaLFTfJzG+DmsRpOEyvy+bHGGEIA2X+XXZ5AhjBGm4zG/LJkeYQJAGCyVH+AVB6hb6\ndXBdTo4wgyDpi/0WI3KEGQRpC1/LBecI0gaeFgseIEjr+VkqeIEgreZloeAJgrSWj2WCNwjSSh4W\nCR4hSOv4VyJ4hSCt4l2B4BmCtIZv5YF3CNIKnhUHHiJI3/lVGniJIH3lVWHgKYL0jU9lgbcI0hce\nFQUeI0jL/CkJvEaQFnlTEHiOIC3xpRzwHkFa4EkxEACCNM+PUiAIBGmWF4VAIAjSHB/KgGAQJH+L\ngIAQJF9LgKAQJD8LgMAQJB+Xj+AQJP8WjwAlHCS+Rwxykg3S/CPzyRG2SzdIcwsiR/hBqkFSg78t\nLhoxIki2l4woESTLC0acUg3S9DkSOcKP0g3SRK8dOcKvkg3SxHUkcoSfJRwkLxaKSBAkl8tENAiS\nu0UiIgTJ1RIRFYLkZoGIDEFysTxEhyDZXxwiRJDIEQQQJHIEAQSJHEFA8kEiR5CQepDIEUTYDNL7\nrNTxXn/I4qdYa97kCDIsBumdqcKp+hAvgkSOIMRikC7q+knTNTuWH2IrSLMP3ZJdDBJnMUhZNeMr\nO7ysBWn+oVuSSwEsBqlpz+/j0V6QFj6OHEGOxSAd1Lv56WgpSLMP3RJcBpBbDdJVneufXuroPEjk\nCJJsdn9f2vTcl7oA9ixi5nPGn0eOIMrqBdnnqfnpdR59itL9vIjhh2p/jl8HpEQ+smGm144cQVjk\nQZq+jkSOIM1FkL4fuZlt6eQI4hIMEjmCvPSCRI5gQHJBIkcwIbUgkSMYkViQyBHMiL7728rnInlJ\nBYkcwZSUgkSOYExCQSJHMCedIJEjGJRMkMgRTEolSOQIRiUSJHIEs9IIEjmCYUkEiRzBtBSCRI5g\nXAJBIkcwL/4gkSNYEFaQfni+EDmCDSEFafFB3ps+CRAWVJC2L54cwY6AgjT/2NRtn/mrRagAAAZf\nSURBVAPIizpI5Ai2xBwkcgRrAgrS1nMkcgR7ggrSpl47cgSLQgrSputI5Ag2hRUkm58AbBBpkMgR\n7IozSOQIlkUZJHIE22IMEjmCdREGiRzBvviCRI7gQHRBIkdwIbYgkSM4EVmQyBHciCtI5AiORBUk\ncgRXYgoSOYIzEQWJHMGdeIJEjuBQNEEiR3DJzya7fRHkCE5FEiRyBLc8DRIQmB9auXxwrHO5Dqku\nO90Vn+FjmbZKtU4Jkkd8LNNWqdYpQfKIj2XaKtU6JUge8bFMW6VapwTJIz6WaatU65QgecTHMm2V\nap0SJI/4WKatUq1TguQRH8u0Vap1SpA84mOZtkq1TgmSR3ws01ap1ilB8oiPZdoq1TolSB7xsUxA\ncAgSIIAgAQIIEiCAIAECCBIggCABAggSIIAgAQIIEiCAIAECCBIggCABAggSIIAgAQIIEiAgiiBd\nDyq7vJ0t3ck2vGQu19nRSpdLdlrZs2II0qX8AoHM0cZ9/vLdBbsdy3U+OFhyyc1KF9xW9rwIgvRU\n53fxT+TZzdIzF23qobJnseiH/UUX3Kx0uWSnlb0ggiCdqnVwU7VXdXSx4Iu6f/68qT/7i86drXTB\naWUv8a5AP3OzbdXFyYJP6pUX/zyf7C86d7bSegm8a7feFehXb3V0sdinm0pVTv9hdrTSHUeVvSSa\nIF3Lgx0X0guS0yUX3FX2rFiC9MrcHOXkBMk+h5U9K5IgvTN3+3qCZJnLyp4VcJD0758+Wr6ioi/b\nRZvKUg6S7cpeJYogvQ7Hl6tlu2lTVa/dy1GvXe4ySPYre5WAg9S6u+3DcdGm/sqz7bu62F90xVmQ\nHFf2rAiC9HK8aRMc2eAuSK4re1YEQTorpR9pWedkwYdyjdPqYSm4ruxZ3hVoO5VikN7l6G8HC665\n2tquK3uWdwUCQkSQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEE\nCRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJ\nEECQAAEECRBAkIxSml++5e5efcjsp2u/vC8HpQ6X97bP9e+r7wLFdjRqZ5AOqvqQ2U/vfr41i7lu\n+lyCJITtaF7TWLc32i9zaG9/cnR55fnrsipJ3YwESQjb0TwLQXpn6l79dFfq+9EdQRLHdjRPC9JF\nZX/lz9eDyupdx+fHw7V6/31QJ+3N+uu7q/kvmTq+ih/uJ1V/oXmXgqtqv+H8ov6ad8o/9clfp3L5\n2udWE3aluR+VOt4NbYiYESTzuiCdmpOY8gd1LF49tj+W71+0N/UglZNln53NX3UqdMn1IJ3Us/nx\nUc3ZLLg3eVb8+DcKUlea6/rzLPQQJPO6IB3fn6Z6KI6/Pj+9j8Xh2E1lz/yZqVv9fu/NLg+34sVz\nkQdVTHor35k+Qut2NNVP2uTN8vW3ewvMikDeikmwDUEyrwvSo/7tVJ7HvIvjuFN5cnOvdiPF+/qb\nXYM/Fe+9Vdb7zDVB6k3+GL09KI1SHNb9hiCZ1+9sqJpv0yXevNn7cfDmsJv6df87bgjSYPLx52oL\nvHwOL5/PHJsRJPOEg3Rs3p0+R3r2d2WjyReDlP8Vp1HZy9CWiBhBMm8qSMM3B5nK5+c4q8P1/hoE\nqe61e76KnUrv5Go0+VSQ9NLeLwfOkbYjSOaNg3TqTkWac6R2N6K/2c1xbM+RyteGQaqvI53Uqeop\nKN95tClZDtJpeGKkaBWbscnMGwep7Kr77EZOg167cjLtTaVe9RzXomftUvXaPfLn8BypuA5bjGz4\nq4/LDupa9MOpicmrF195FyRtgYeqj4890mYEybxxkOoTl7LN69eRqum6Nw+q3Qe115Eu9QnNo7/n\nuLdnOsVFoPJ60KnuPuhPXvzZfe6gNLd2YmxDkMybCFIxlkCdq3P6a9aObKhnaN98HNoglR1q5Qzn\nT+we+rFgpR79fT+WgyP+MnWu3h9OXvzZfe6wNOXIBnK0HUGKzv3PdQlSRJAAAQQJEECQAAEECRBA\nkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQ\nAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQ8A9gxSrZ7s4dhwAAAABJRU5ErkJggg==\n"
          }
        }
      ],
      "source": [],
      "id": "41d2f912"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here’s an example of what a q-q plot might look like if the observed data are a random sample from a distribution very different from the theoretical distribution being considered."
      ],
      "id": "3cdc6219-da35-449e-b423-d6baac95d92d"
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyMjIya\nmpqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAASdAHeZh94\nAAAgAElEQVR4nO3diZaiMABE0eA2aqv8/9+O7AFBRSobvHvOTGsrJJJUA2HR5AAWM6ErAKwBQQIE\nCBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIg\nQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAA\nAYIECBAkQIAgAQIECRAgSIDAyoNkzNgHPP443cR7Ps/uc/lTBR6/rs5XZX6c0fgb6t9+tTTfvnd2\nfZKxls8xYayh/rIZIfn8nm9m97n88QKbty7ub1/PaHaQXj7NV0GSfbBorOVzTBhrqFlrm8/vmd0X\nRid424GX97evZzQ7SC+/+2qxyD5YNNbyOSa4DNLOZN++9fO8PQXpxze+mXzqpe+CtB5r+zyN+zEz\nu3PbYJfD89HudC+b0NS/7X45OV1+PT4fHq/Vk+K31/3z+a14djCH3uy695x3ZveX5+fM7P/y6fIf\np50x+2s72V8x63tvVvVbh68+TpnJ2ooP63jfmVPvTf0ZDYvuLYdBD+8vjqbe/55VMYfLcM5VudZ7\nn8shq6vc+5PwWp93yzkNKw3SX9VU+7qh9nXTmT+rEa1fTk3XvudQPns+OFXPixb+V3eblyBVE91P\n3cxHy79n1cNTb9bZ3Z5V19/sV5sp/8brWGSk96aXjtsrurccep9luDj69S4K6c151/yqtxyqWpr3\nQXq3nNOw0iA1bV011PnZqI88PxUN1f7W/uXEdPmhfVa+p3uxHamy3tz7xbPPdxOOlt+VdR2fdb+/\n2a82U2aTdbz03vTSce2i+8uh91mGi6P6cSxn/3j2/fNLFS92kOxamrdB+ryco7fOIF2ezffsI9es\naqjnn0pru6VuvP4vR6e7Pn+cH89Nma63Fy8fzSA5dtGm7JfP3ml2t/LHVPnP17Jb2R939az/yu6c\n5a/z7r9adf1H1aFH6ljkov+mwc5Wr+jRyo0uju4txQSPtuLW587tIDXTvlZg8PSr5Ry5dGo6x6Fq\nkLKFrF+P9JX+08F0x+KPbuFU/W2s2/nxIUh/vR9T5ddlPXb/7lOztvub9epzyrrLHibr2H/ToOP2\ni55eOMPFaNo1abM381JFOyrNtK8VGDz9ajlHLp2aztG2QPvgfjntzbCvWL8cnc7UnTG/96f8EKTX\nH2PlT0bmzW/bGdSy0To2K4buTa8z6C2u0YUzshirH/+q2R57uenqMfx4o4vupT6fl3Pk0qnpHMMe\ncNk1vcp6sffL0em6lxYGaax8RZDM2zqOfOaxII0unNHFWP84NQm9j1X8tyB9Xs6RS6emcwxa8bmt\nb3bH863XTP1fjk5n/6V8/bs+nKb3tPdjtPzfg5QNih+tY/ZmRr1CRis3ujja549LNcq2H6t49157\nDfUpSJ+Xc+TSqekczcb9pWqKXf2016b9X45Od3jdds+Hk3wTpNHy98N9pJG5jQepqWVT5dE69t40\nmEGv6NHKjS6O3ovX48Tao3vvpXxfu49UZOVvNEhfLefIpVPTOc7VsM8l6zW29Uf3Mfzl6HQjo0nl\n216C9Og/Hf4YLf9l1O511u28+69eqjG8S7lKmKpj702DGfWKfq3c1OLo/izUe2HZSBWtT14OKj6n\n/ZeXK8hTeYpdN6duqq+Wc+TSqeks/QMg+7IVm1Hc4sdp+MvR6bpjldUBjdEGrmfXGAvSaPldWeep\nWbdvHbzaTvk3Xcf+mwYzsovuV+61fGtxtN1+fy/HHE4jVbSDVMmKxBz7cxpM9c1yjlw6NZ2lPiR/\nqJrir2nSslcVbbof/nJ0uq6F+yf99xq4nl1jLEij5Vd/nk13ZsPrrNu3Dl69GmvKiTr23zSYkV10\nv3L98geLo/7RDDbsx6poffL6MGu5dXivC6xfHU71xXKOXDo1nef+bKp9e5LY7fksO97u1fb6oWqw\n/i9Hpyt2BbL+OWD2z8qhdwB+LEjj5Zcnw5nDyyiy/TEO4zsN5alyh2YfaLSOgzcNZmQX3avcoPz+\n4mheLPeP9ufROdvvPVunBxal7C/2Llxvqi+Wc9zSqSkQMYIECBAkQIAgAQIECRAgSIAAQQIECBIg\nQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAA\nAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIEPATJAIn5oZfr\ngxOgCECJIAEC/oN03hlzuDotAvDNY5Cqzch9tUV5clIEEIjvIJ3M6ZHn95M5uygCCMR3kDLzKB4/\nzM5FEUAgvoPUDBO+DhcuHEsEQvIdpGMTpMxFEUAgXoN0+He+msvz4eP0frSBICExXoPUbrYZkz1c\nFAEE4vM40u12Ph8O5ZDD6W2OCBJSw5kNgABBAgQIEiBAkAABggQIECRAgCABAgQJ+Nabs0AJEvCd\n3lnXLy/+ML9FtYmlCGAmY/0/8eIP83OKICE6ZvBz/NUfZugQQUJ0CBIgQJAABfaRAAFG7QAJjiMB\nbhEkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAk\nQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAg\nAQIECVv15hthf5iZl0kiLAIb9/Y7yn+YnZdJIiwCG2es/2Wzcz1JhEVg28zgp2p+bieJsAhsG0EC\nBAgSoMA+EiDAqB0gwXEkIDYECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIg\nQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAA\nAYIECBAkQIAgAQIECRAgSIAAQQIECBIgQJAAAYIECBAkQIAgYd2M8dKbCBLWrEyRjygRJKyZsf73\nUJDrSSIsAqs02JAzg58OC/YySYRFYIVeNuQIEjDfy4YcQQJmG4kN+0jAXGNBYtQOmGl0Q27Nx5E+\nfjSChF9425CbKtr1JMM5ECS44G1DbqRoL5NU0/W5KAJb52lDbqRgL5OU/jKChLXyuWn3OJj9vZzD\n2Cy+ThkQH7/7SBdjLjn7SFgfz4MN9705PAgSVsf7qN0/k10JEtbG//D3bfd5H4ggITEhjiMdCRLW\nhlOEAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQ\nAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGC\nBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQ\nJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECA\nIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAEC\nBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEhImjFxdBaChISVKYoiSgQJCTPW/2ERJKTL\nDH4GRJCQLoIUvgisAEEKXwTWgH2k4EVgDdYzanfe5fl9Z3Z/qgq9FgFMW8lxpGvxMbLnhzHSJMWx\nbICvLQzS3lzym9nlF7OXVSknSEjOwiAVK6SbOam3UwkSEiMI0sFcCRI2bvGm3e1qspxNO2zc8sEG\nY/4VK6SrrEo5QUJyFg9/Z8UeUr67fD/9c5LdWV4rICSfB2RvB5Od83/FSuzDliBBQmI8BulWJuhk\njo/8fjBv10kECYlZHKTroRy5u3+e7lhsBJ6KoYk8f5iduFZASEuDtDflORom+5ykaoTcHKwn/Zct\n82sFhLQwSGezfxTd/myOn6crJ7xU23TViklYK2xHjH9pFwYpM4+vz8A9FntHlUe5mSetFbYinjO+\nbYIzG779YI+sfZN5v0IiSJgWzzVItoVB2tVrpNv7wYPaqYlP9nZ9FN9SQjwiuirWptlHumbvh7OX\nFAHY1hmk/FAPs0lPtYttISEiKw1SeRzJHGacITS/CMCyyn0kR2JbSojIKkftHIltKSEqKzuOZIyz\nkxHiW07AWwQJEGDTDhAgSIDAok273uZd4FoBIREkQIBNO0CAIAECgssoStn76yIWFAEkQBSkO/tI\n2LQFQbr2jsd+cz2Sy1oBIS1ZI+3sHPG1Ltgy1T6SFkFCYhi1AwQIEiCwNEj/dpzZACwN0j9OEQJ+\nGyzo3yBSevegsSKA6BlG7YDFzOJNu4N5TL5vAYKEhJh8cZDu2V56JHakCCBupv3vh+maxww2YON+\nvs8eQQJapvfjhyldIkhIxII7KBMkoGZeHvwwre3v8HtlviwCiIy1Q7M0SCf2kbBVZuLxD5N3Obou\nrdVUEUCczOST+dNn5pLvzf2+N1zYh20xb57NnkGxRffvuTa6ab9pjCAhdubt07lzKIJ0LU5cZR8J\nmzLsoovPtbvkd7PL/wgStuSlhy4M0rUI0L4YbDguqNXbIoDovHbQxVfIFs+Oxpx+rdLnIoDIjPRP\nzmwAZhrrngQJmGe0dxIkYJbxzsllFEhEJN9lPlEJgoQklB0sgihN1UCzafe3l578TZAw9PN1qFqT\nFRDtIz04jgSXFlw/pzRdvGqwgU07uBRHkN6ULgrS2fCNfXAoiiC9K1w22PBv/oy+KwLIo9hHelu4\nKEg77Z2LCRIGgo/afSibA7JIRNjjSJ/KJkjAZx975LIgPf4Vl1Acivs1HJQ7SQQJUfncIRcF6dzu\nIt3P0rufECTE5Iv+uCRIzxydinue3E4m454NWK1vuuOCID26e3BduB0XVuur3rggSKfustiD+BJZ\ngoRofNcZFwRpZ+71o7v5Z3bzZ/S5CCC0L/vigiB1w/q3A+faYZ2+7YqSIA2fLEaQEIeve6Jk0y4v\n720nRJAQhe87omawofdYgCAhBjP6oWb4+2q0325OkBCBOd1wyQHZizGnW14ekNUeRiJIiMCsXrjo\nFKFLVp8ilGlzRJAQ3rxOuPCk1fOhOGn1LN2uywkSwpvZB7mMAhgxtwsSJODV7B5IkIChH04vIEjA\nwC/djyABfT/1PoIE9PzW+QgSYPux7xEkwPJr11scpOuhGOI43Cfe/BuChDB+7nlLg7SvvhrJZNIk\nESQE8XvHWxiks9k/iiCd+VoXpG9Bv1sYpMw8XNyVmSAhgCXdbvFN9J3c3pwgwb9FvW5hkHb1GunG\npeZI3LJOp9lHumZG+r0uBAm+LexzS0ftDvWlfdI7FhMk+La0y0mOI5nDZWE13hYBuLa4x3FmA6Lm\n5+vFlpdBkBAxT194KShg0Z1We5bXZVGtsEbG+t91Kd7nQZDgiRn8dFuK75mwaQdPvARJM3OChHh5\nCJJqU2pxkC7F8PeRG0TCBef7SLJZSy6jKI4kqSr0WgS2y/WonW7OC4N0qu5WzClCcMTpcSThrBdf\nRnErf3LSKtKj7GaCyyj6DyQIEtwL3WX7m3bNGkm6k0SQ4Jy2ky0dbPhX7iP9ZZz9jbSI+9jiTTsn\nZzcQJDim7mIECVsk72Gc2YAN0ncwgoTtcdC/CBI2x0X3WhqkU6beP3opApBy0rsWH0fieiSkxU3n\nWjxqJz3HbqwIQMlR31KdIqRFkDbP1cmqrrrW4k27h6wqE0Vgg5xdPuGsZy2+Hmmv/WakkSKwPa4u\n6IvqyqbeJFcGGyDn6BLzyC5tsif5x6gd9NwEKbZbqPQv7Jsxavc4GrOv7+7wPngEaducBCm6e3r9\nOmr3yKy7OxAkjKq2bRzsIznuU4s37b4ftTsVa6/Hubp0iSBhRDNapx+1c92lFl/Yt//7drqsmvCe\n7e4ECaO6NVFqe92665E+T1e/5bHfjwXJ2f2PkQxHo3U+/jJ7DNKu3Qzc7VkjYYSrIMV5wcKvtTqb\nY/3obvYECa9cHT/Szk5Vxs/VOrXpuX5YgxGkjXJyRoOX3qQK0t83t+O6tW+6HwkSXrk4x85PZ1p8\nYR9nNkBJPtLkqS/JLuyTfh8FQYKIr660+BShS7439/vefH04aW4RwO+89STBKUL/nmujm5HeapUg\nQcJfRxIE6Vqc+sM+EuLjsR8tDNLhuWl3N7v8jyAhOj670cIgXYsAld/ad5x8/w8I0pbUJ3zLTwzz\n2osWn7RaPDsacxLVZ6QIrFWZHStFKW/W+DyzIaoiEFRz6LFZD8lPafB94jNBgn/W+qdJkDhJ3nvQ\nkiA9TuXDv53JxLeJJEjrVl0EWwWpjpM2SP470JIgZeVCqO4jxDf24WvNqsheL0mDFKD/LAjS2eyL\nC4yy7JY/9uYSuFZIhx2kXoLSzdGSIO1NcW/IP/Ov/J8zG/AV040v2ON1wlG7IL1nQZCqj32qzrJL\neeQS/tjDDM1D8XGkMJ1ncZB2xnqiQpDWyuotjm7OEajvLAjSrti0u1enNDxMJqwUQVqr7o+wq+M8\nobrOgiCdisGGY3Uh0plThPAFM/jprADvFgSpunNqOchwNuYmrBRBWivXQQrXcRYdkG3OsTOca4fv\nyE8FGpl7EJJThMxBen0sQVovZ98gVs7czWxdlc25dljA3X10g/YaggRv3N6LOmynIUjwxOlWXfA+\nQ5DgyWrHGX4tnyDhB25HvoP3GIIEH7rdIydtG77DECS4152a6qZtI+gvBAnudWd7u2jaKL6XjiDB\nueaUb+Ok08fRWQgSnLPO+XY388AIEpxzOc4QS1chSHBv1cMMFYIE95yd1BBPRyFI8MHNaXYR9ROC\nhGTF1E0IElIVVS8hSEhUXJ2EICFNkfURggR3jH2SnXjWDua5BEGCK92gtz5K0fUQggRXHJ6qGl8H\nIUhwpEmQgyRF2D8IEpwwxl2QYuweBAkOtIMMLoIUZe8gSHCg2z3S7yPF2TkIEvSGV/IpR+0i7RsE\nCXq9K/mkx5GiuKx8DEGCXBud9R8+ahEkiFnfZimftXqGOgQJYtVOkXGwGRZztyBI0Gq26jZwWpCN\nIEHLDH7KZxwnggQtV0GKvE8QJIgZ63/1XONFkCDm5MKJ6HsEQYJUMcigH2iIv0MQJAht5TK+VwQJ\nQtvcPyoQJOi4GbFLojcQJKg4+la+NDoDQYKGq2/lS6QvECRodNfybTFHBAkaw2v5pLNNAEGCRO9a\nPvlc40eQILHhcYYSQYKGg2GGlLoBQYLGVu5yMoEgQWWru0clgoQopdYFCBJilFwPIEiIUHodgCBB\nYdP7RwWChOXEI3YpNj9BwhL2BbGqVkuy9QkSflef8W2kJ6um2fgECb8rt+jy+nuQRM2WaNsTJPys\nXBflee/f4nkK5hECQcKP+l9/ZLadI4KEX1lro3LNpBi1S7fhCRJ+0nxDbJMnyeB3wu1OkPCDehXU\nfefylscZSgQJP2h3ipQ3Vk261QkS5msGu0UjDNZMk0WQMF8VJNUIgzXPdBEkzFe3j/JM1dSbnCDh\nB+obNCR1e4ZRBAmzGcPp3kMECTPVIWKzrocgYSb1Zt06WpsgYZZ2TcR2XQ9BwgzWkPemL+N7RZAw\ng3UQVtRGa2lqgoTvme6CCXLUR5DwLdNdWa4a+l5PQxMkfKu9AEk29L2idiZI+FJ7GR/jDCMIEr7T\nXnfE+mgMQcIXuh0j3RkN62pkgoQP2uGF7iYNktmK5hMJgoR3uqvJhbc4KWesmU00CBKmdRt05TO2\n66YRJExqLpdozvhmu24aQcIEa5iOce+PCBJG1WML1iqJ/aN3CBJGdLmxbl2nitEqW5cgYYRp94ys\nW0HKZr1GBAmv2rOBmrNUtbNeI4KEV80FR8oBb2vOa0SQ8KoKkvYGkN2MV4kgYUTVAOrV0ZrblSBh\nhPi+dc1cxfOLiccgmT4XRUBGvjZaeat6DNL5fZC+ThnStO5G9blpd8v2rouAhIs/ZStvU6/7SDdz\ncl0ElnOyg7T2JvU72HA2N9dFYJFuy5px71kYtUPH1GcDNaeq6mYsnFecCBI69ckM3TlCuvmuHUFC\no7lyojllVTdj2ZziRZBQGVw5QY7mIUgotasj6xIk0ZxVM4oaQUJuX7onP+d7I21JkNBec2R6o3ai\nWcvmFDeChPYuJ/W1sMoYbaYlCdLmdTfdkg4x5Pmm2pEgbVz/plvaNciWmpEgbVh3qn29h8T66GcE\naYuaBLUjCw5WRxtrRIK0LVaCuhTlDlZHW2tDgrQh/ZVQc8zIyepoc01IkDajOQeo25qzNurUS3xz\nLUiQNqO5l7exI+VmdbTBBiRIm9COz7XrJGvbjtucCBCkDbB2jezNO2u8QVyeeoYJIEgbYF+v16XI\nOowkL257CNL6dVdH9CPkrrgNIkjr17/MyPVtAzfadgRp/eoguRjlHi1skwjSutUbcvLbMEyW57yE\nSBGkNbO25tyMz70U6LqAaBGkNWuvMfJ0R/UNtxtBWjEz+OejvK0iSCvmOUibbjWCtGJeg7Sd2zOM\nIkgr1pxdl/vIkesCIkeQ1qp3EpDz1cXmW4wgrZLprpRwcnr3S3muC4geQVqfNkUM1/lDkNbGvmrP\nz9kMNFdOkFan3SfitCCvCNKq2ClSfzvLZJnOS0gBQVqR5uzU7vo9Tq/zhSClzzortf2msPoJ6yNf\nCFLqrMNF1m28/aSIlmoRpNQ1QwrdPVM9XcJXF44SQUqclaPuoa/lRzu1CFLiBkHyM8DQlY0aQUpc\nG6RmuNvfWdi0koUgpc1Y98P3ujaikfoIUspe7sngccHRRj0EKWXVcuqS5L1oNAhSwprsBFhctNAA\nQUpWd8We98W18cvKxxCkZHXXkPteXDTPK4KUqi5F5CgCBClF7Y2ITe53qK4s3G9xiSBI6Wnj43vE\nuyrdc3mJIEjpqTfo2iOx3gvHK4KUnHbniO26iBCk5FijDGzXRYMgJae+oDzEwDftMokgJaW9oNzT\njYgHpfstLikEKSHtOardTbe8Fu+1tMQQpIRYV0wE2EGiUd4hSOkwg3/+S8ckgpSOgKMMNMknBCkV\npjuhgRzFhyClot078j/KQIN8RpDS0ObH+6WwZem+C0wPQUpB94VHQS6pozk+I0jx6+0YBVg0tMYX\nCFLs2i8OCzPKwGXl3yFIsWsC5Pu2dV3x+AJBipx14CjEDhIt8SWCFDdjfRUsu0cRI0gRa787LNBm\nHe3wPYIUqeEFEyH2+WmG7xGkCJnudKBurRSgGgHKTBZBik73lXv1gaMwR2G33QizEaToNCfV5bk1\n0BCkGvgeQYrN8Koj/3dSbauB7xGk2LQBMnmw0bptt8BPCFJsmjPq7BGHMJXADAQpHoOxujzEBRN1\nTcIUmzKCFAfTBCh4iPItLv7lCFJQVmisS19Dhijf0NJXIkgB2SshU4/PhbroyKpVyMKTRZACqlc/\nzWGj3MQQpI0sezWCFE43PpdbI97tKwErhdkIUij1jlB/bWTIUaIIUhj9m9RZIw55HuYIbF2tYCWn\njiCF0ewevYx4Bx32Xv9yd4YgBVF9QHslFEEnDl+DhBEkvwZrneArIUsUlUgWQfLGDM9eyKP6oBFV\nJUUEyZNuGMEaUIjoc0ZUlSQRJE+qY0XNMJ19QlAU4qlJogiSB8Y6ZtSdBhRRjNa2wAMgSI402bF2\njAZBClm7gZjqkiiCJFd9wWveroh6+0Qmj2/3KKq6pIogidnnJph2x2hisCEOEVUlXQRJx7qFVi9I\nuRWfmA4c1WKqS7oIkkpvVWSvhbrBhSg/VpSVSg9BUjG9/60VU/dV5DGKtV6pIUjLdZt0/Y06+/84\n10apLemIEaTf9Qbl2rPn6pdya9Qu0o26PJUFnQKCNJfpjSQYOzm9X0dxH5MPYq5bYgjSd4bxyUcD\n1J39k4RU6pkCgjTNvKxmrCGF/qZcbzsvFSnVNXoEyS7WzkF/OLv+v4vPS5BSWhWVkqps9AhSV2h/\nlWJFaCw+r5t0aUkr9PEjSINCR9Y04+uh9ohRcpt0hdTqG70NB2mwJTZI0HiQrLSNjX0nI7kKR2+z\nQXpZkbwJ0mh8kgxQLdV6R2z1QZrq7S+7NmPbdIPjq29nmJDkP0CE0grS7E48uQMzjE3+Eq3+pCuI\nT2s9nyQiXoP09+9Qni9zOP39UsQPu/WTQ2pjQXrZ2FtReCyr/FDBeQzSY2c6+x+KmD/QPBKXd6+s\nNDl9G/iIIXgM0slkl1v56H7NzGl2EdOp+DifyVXS9nrV9j6xHx6DlJlb+/hmstlFiIOU5OGfxTb3\ngX3xGKRer33twsb2rtz523ZTu1zb61Xb+8S+JLRG+mVrbKPrnSksCGf87iNd7+Wj3/aRfkvFFtc7\nU1gS7vgc/t5b2267xy9FkIolWHYO+T2OdCqPI2WHfz8dR8IiLFSX0jqzAb9jmTpFkDaCReoWQdoG\nlqhjBGkTWKCuEaQNYKzTPYK0fixNDwjS6rEwfSBIa8ey9IIgrRyL0g+CtG4sSU8I0qqxIH0hSGvG\ncvSGIK0Yi9EfgrReLEWPCNJqsRB9IkhrxTL0iiCtFIvQr0iDBCTmh16uD453IT/DVsve7gefEGOd\n5tpqmxKkiMRYp7m22qYEKSIx1mmurbYpQYpIjHWaa6ttSpAiEmOd5tpqmxKkiMRYp7m22qYEKSIx\n1mmurbYpQYpIjHWaa6ttSpAiEmOd5tpqmxKkiMRYp7m22qYEKSIx1mmurbYpQYpIjHUCkkOQAAGC\nBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBFYRpPPOZKdHsNKD\nLMNTFvIzB/rQZclBG3vSGoJ0Kr9AIAu0cG+/fHfBYvvyM+8ClFwK86ELYRt72gqCdCj3FFYAAATS\nSURBVDPHR/En8him9CxEn/oz2a0o+s9/0YUwH7osOWhjv7GCIB2qzxCmac9mH6Lgk7k+/7+Yf/6L\nzoN96ELQxn4nugr9LMyyNacgBR/MPS/+PB/8F50H+9B2DaLrt9FV6FcPsw9R7C1Mo5qgf5gDfehO\noMZ+ZzVBOpcbOyFsL0hBSy6Ea+xJawnSPQuzlZMTJP8CNvaklQTpkYVb1xMkz0I29qSEg2R///Te\n8xEVu+wQfSrbcpB8N/ZXVhGk+25/D1V2mD5VjdrdA43a5SGD5L+xv5JwkFrXsGM4IfrUv3Jv+2pO\n/ouuBAtS4MaetIIg3QMv2g2e2RAuSKEbe9IKgnQ0xt7S8i5IwbvyE29rhKUQurEnRVeh+cwWg/Qo\nz/4OUHAt1NIO3diToqsQkCKCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQ\nIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiA\nAEECBAgSIECQAAGC5JSx/PItd9dqJpNzt548TjtjdqfHvPnG99V3iWI5OrUwSDtTzWRy7t3jS1PM\nedZ8CZIIy9G9prPO77QfprBefubodM/z++mrJHUTEiQRlqN7HoL0yMy1enQ15vPWHUGSYzm6ZwXp\nZLJ/5ePzzmT1quP5cHeuXn/szMF6sf767mr6U2b29+LB9WDqLzTvUnA27Tecn8y/5pXyf/vt90NZ\nvjXf6o1dba57Y/ZXRwtizQiSe12QDs1OTPnA7Ivf7tuH5esn60U7SOXbsufK5l+1K3TK7SAdzK15\n+FdN2RTce3tWPPz3EqSuNufv97PQQ5Dc64K0fzy76q7Y/no+euyLzbGLyW75LTOX+vXei10eLsUv\nj0UeTPHWS/nK+BZat6KpHllvb8q3X+4VmBWBvBRvwTwEyb0uSH/1s0O5H/MotuMO5c7NtVqNFK/b\nL3Yd/lC89jBZb57fBKn39r+Xlwe1MYbNut8QJPf6gw1V922GxJsXew8HLw6Hqe/Xf/sZQRq8/XW+\nVoGn5+bl7ZZjNoLknjhI++bV8X2kW39V9vL2t0HK/xW7Udnd0ZJYMYLk3liQhi8OMpVPT3E0u/P1\nPghSPWp3uxcrld7O1cvbx4Jk1/Z62rGPNB9Bcu81SIduV6TZR2pXI/aL3RT7dh+p/N0wSPVxpIM5\nVCMF5St/bUreB+kw3DEy9IrZWGTuvQapHKp7rkYOg1G78m3Wi8bc6ynOxcjaqRq1+8tvw32k4jhs\ncWbDv3q7bGfOxTicGXl79ct73gXJKnBXjfGxRpqNILn3GqR6x6Xs8/ZxpOp93Ys7066D2uNIp3qH\n5q+/5ri2ezrFQaDyeNChHj7ov734v5vvoDaX9s2YhyC5NxKk4lwCc6z26c9Ze2ZDPUH74t+uDVI5\noFZOcHzG7s/eFqzUZ39f9+XJEf8yc6xeH769+L+b77A25ZkN5Gg+grQ613+ha7BFBAkQIEiAAEEC\nBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgS\nIECQAAGCBAgQJECAIAECBAkQIEiAAEECBAgSIECQAAGCBAj8B1Y9H5FV7f9EAAAAAElFTkSuQmCC\n"
          }
        }
      ],
      "source": [],
      "id": "492d6f2b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Q-q plots are often used to evaluate whether data seems to have come from a normally-distributed population. However, each random sample from the population will vary. How much variability can we expect in our q-q plot if the data are a random sample from a normally-distributed population?\n",
        "\n",
        "In the code below, we:\n",
        "\n",
        "1.  Use the `set.seed` function to make our example reproducible.\n",
        "2.  Use the `mfrow` argument of the `par` function to create a plotting grid of size 2x2.\n",
        "3.  Use a `for` loop to iterate over the values 1, 2, 3, 4. For each value of the iterator:\n",
        "    1.  Use the `rnorm` function to draw a random sample of 20 values from a normal population with mean 0 and standard deviation 1.\n",
        "    2.  Use the `qqnorm` function to create a q-q plot for the random sample based on a normal population.\n",
        "    3.  Use the `abline` function add a 45 degree reference line for convenience.\n",
        "4.  Use the `par` function to reset the plotting grid to size 1x1."
      ],
      "id": "3336cb95-d7f2-436f-bf26-aa57ac5595a7"
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {},
      "outputs": [],
      "source": [
        "# make example reproducible\n",
        "set.seed(33)\n",
        "# create 2x2 plotting grid\n",
        "par(mfrow = c(2, 2))\n",
        "# draw normal samples, create q-q plot\n",
        "for (i in 1:4) {\n",
        "  # draw 20 values from a N(0, 1) population\n",
        "  s <- rnorm(20)\n",
        "  # draw q-q plot for sample using normal distribution\n",
        "  qqnorm(s)\n",
        "  # draw reference line\n",
        "  abline(0, 1)\n",
        "}\n",
        "# reset plotting grid to 1x1\n",
        "par(mfrow = c(1, 1))"
      ],
      "id": "3ec4d784"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### KNN example ([YouTube](https://youtu.be/oDn1NBKB_Ss), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=1b26f627-78b5-47ae-9bbd-b0ad0106357a))\n",
        "\n",
        "We now use a `for` loop to determine the K nearest neighbors of a set of locations. For each location, we determine the K locations that are nearest to the current location. These K locations are the nearest neighbors of the current location.\n",
        "\n",
        "We start be generating a random set of 25 coordinates where the x- and y-values can take any value between 0 and 1. In the code below, we:\n",
        "\n",
        "-   Use `set.seed` to make sure the example is reproducible.\n",
        "-   Use the `runif` function to draw a random sample of 50 values between 0 and 1 based on a uniform distribution.\n",
        "-   Store these values in a `matrix` with 25 rows and 2 columns that has the name `coords`."
      ],
      "id": "0480361b-7cd2-4841-a69c-93cce13b2502"
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {},
      "outputs": [],
      "source": [
        "set.seed(1989)\n",
        "coords <- matrix(runif(50), ncol = 2)"
      ],
      "id": "9d0118b2"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We plot the coordinates using the `plot` function."
      ],
      "id": "5df915a8-0c5d-46e9-a1d6-e1c31ddfaef1"
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {},
      "outputs": [],
      "source": [
        "plot(coords, xlab = \"x\", ylab = \"y\")"
      ],
      "id": "d7297695"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will determine the 3 nearest neighbors of each observation using a `for` loop. We will store our results in a `list` with 25 elements. Each element of the `list` will be a `numeric` `vector` of length 3 indicating the 3 nearest neighbors of the current location. We could store this in a `matrix` for efficiency, but we use a `list` to get broader exposure to the R programming language.\n",
        "\n",
        "We start by creating `nnlist`, which is the `list` that will store the nearest neighbors for each location. We can create an empty vector using the `vector` function, specifying `mode = \"list\"`, then indicating the length of the `list`."
      ],
      "id": "87317c06-8509-42b6-b04f-cb73cbece129"
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {},
      "outputs": [],
      "source": [
        "nnlist <- vector(mode = \"list\", length = nrow(coords))"
      ],
      "id": "8023eba4"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Why did we specify `length = nrow(coords)` instead of `length = 25`? This allows us to adapt our code to a new set of locations that may not be the same size.\n",
        "\n",
        "*It is wise not to “hard code” values (like the size of a `matrix` or `list` or even the values for a loop) into our code to make it easier to use the code for a new setting.*\n",
        "\n",
        "To find the nearest neighbors of each location, we need to find the distances between locations.\n",
        "\n",
        "The `dist` function can be used to determine the distances between the coordinates contained in a `matrix` object. The function can compute many kinds of distance, but the default is `\"euclidean\"`, which is what we need. By default, the `dist` function doesn’t return the diagonal elements of the distance matrix nor the upper diagonal elements. A trick to get the complete distance matrix is to use the `as.matrix` function on the object returned by the `dist` function. The value contained in row `i` and column `j` of the distance matrix is the distance between the $i$th location and the $j$th location.\n",
        "\n",
        "With all that being said, we compute the distance matrix for the locations in `coords`."
      ],
      "id": "ef098371-d31b-4e43-8344-580e331953a3"
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {},
      "outputs": [],
      "source": [
        "d <- as.matrix(dist(coords))"
      ],
      "id": "80da4cbe"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This matrix is too big to look at easily. However, from the code below, we can see that the distance between the 3 and 17th locations is 0.474391."
      ],
      "id": "d06ea623-b342-4c27-aedd-f860c5a69abc"
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {},
      "outputs": [],
      "source": [
        "d[3, 17]"
      ],
      "id": "38cdae46"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "How do we find the locations nearest the current location? We need to find the locations that have the smallest distance from the current location. The easiest way to do this in R is using the `order` function.\n",
        "\n",
        "The `order` function takes a vector of values and indicates which positions of the vector contain the smallest values, in ascending order. The `decreasing` argument of `order` can be changed to `TRUE` to order the values in descending order.\n",
        "\n",
        "Here’s a quick example that uses the `order` functions to order the values of a `numeric` `vector` from smallest to largest."
      ],
      "id": "2592f6ee-4a44-436b-a0ee-25ed867863f6"
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {},
      "outputs": [],
      "source": [
        "order(c(5, 1, 2, 4))"
      ],
      "id": "1a5aa017"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The second element of the vector we gave to `order` has the smallest value. The third element has the second smallest value. The fourth element has the next smallest value. The first element has the largest value.\n",
        "\n",
        "We can use the `order` function to order the distances associated with each location. Let’s do this for the first location in `coords`. We pass the vector of distances between the first location and all other locations (including itself) to the `order` function. This vector of distances is simply the first row of `d`."
      ],
      "id": "f835ed3b-96ff-44ff-8710-9411fb0b639b"
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {},
      "outputs": [],
      "source": [
        "order(d[1,])"
      ],
      "id": "927f181f"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Unsurprisingly, the first location is “closest” to itself. The next closest location is location 23. The next location after that is location 21, then 21, then 15, etc. A K nearest neighbors algorithm doesn’t usually allow a location to be its own neighbor. So, in our case, we want to select the second through fourth elements of the ordered distances. We can do that using the code below."
      ],
      "id": "1486fa3f-af57-4ed6-9eec-8c4a4451e14c"
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {},
      "outputs": [],
      "source": [
        "od1 <- order(d[1,])\n",
        "od1[2:4]"
      ],
      "id": "64a34aea"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we will use a loop to repeat this process for all locations in `coords`. For each location we will order the distances, determine the indices of the 3 nearest neighbors, and then place the nearest neighbor indices in the appropriate elements of `nnlist`."
      ],
      "id": "5555768a-f1cc-46c4-a38e-3b1b06a5265d"
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {},
      "outputs": [],
      "source": [
        "for (i in 1:length(nnlist)) {\n",
        "  # order the distances associated with the ith location\n",
        "  odi <- order(d[i, ])\n",
        "  # extract the 3 nearest neighbors for observation i\n",
        "  # (not including the location itself)\n",
        "  nni <- odi[2:4]\n",
        "  # put the nn of observation i in the ith position of nnlist\n",
        "  nnlist[[i]] <- nni\n",
        "}"
      ],
      "id": "194fa898"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We look at the first 3 elements of `nnlist`."
      ],
      "id": "c7fde65a-1f24-4960-98b7-369dbae29548"
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {},
      "outputs": [],
      "source": [
        "nnlist[1:3]"
      ],
      "id": "c3a6c60c"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The nearest neighbors of the first location are locations 23, 21, and 15.\n",
        "\n",
        "The nearest neigbhors of the second location are locations 11, 25, and 18.\n",
        "\n",
        "The nearest neighbors of the third location are locations 21, 12, and 15.\n",
        "\n",
        "### Nested `for` loops ([YouTube](https://youtu.be/AvO74v6_0tc), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=6eeb144f-d690-402d-bea5-b0ad0106c0f9))\n",
        "\n",
        "A nested `for` loop is one that uses one or more `for` loops inside of another `for` loop. Here is some pseudo code to illustrate a nested `for` loop.\n",
        "\n",
        "    for (i in i_vector) {\n",
        "      # some code\n",
        "      for (j in j_vector) {\n",
        "        # some code\n",
        "        for (k in k_vector) {\n",
        "          # some code\n",
        "        }\n",
        "        # some code\n",
        "      }\n",
        "      # some code\n",
        "    }\n",
        "\n",
        "The psuedo code above has 3 loops. The outermost loop is iterating over `i`, the innermost loop is iterating over `k`, while the second loop is iterating over `j`.\n",
        "\n",
        "The innermost loop iterates the fastest. This is because for each value if `i` and `j`, we iterate `k` through all the values in `k_vector`. The second loop iterates the second fastest because it iterates through all the values in the `j_vector` for a single value of `i`, but not as fast as innermost loop, which will have iterated `k` through the all the values in `k_vector` `length(j_vector)` times!\n",
        "\n",
        "It is helpful to use proper formatting to make the code easier to process. Make sure to indent the code inside a `for` loop to indicate that it is run instide a certain `for` loop.\n",
        "\n",
        "We will use a double loop to compute the distance matrix of a set of coordinates.\n",
        "\n",
        "Recall that if $(x_i, y_i)$ and $(x_j, y_j)$ are two points, then the Euclidean distance between them is computed as\n",
        "\n",
        "<span id=\"eq-dij\">$$\n",
        "d_{i,j} = \\sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}.\n",
        " \\qquad(1)$$</span>\n",
        "\n",
        "Let’s create a small set of 5 coordinates."
      ],
      "id": "a290552b-b0f0-414b-9b65-ed518fcc1070"
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {},
      "outputs": [],
      "source": [
        "coords <- cbind(c(1, 4, 2, 5, 3), c(9, 4, 3, 2, 1))\n",
        "coords"
      ],
      "id": "07ff9740"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Each row of `coords` represents a location. For each row of `coords`, we need to iterate over all the locations of `coords`. For each combination of locations, we need to use [Equation 1](#eq-dij) to compute the distance between the locations.\n",
        "\n",
        "We create a matrix, `d`, to store the distances we compute in our loop."
      ],
      "id": "04d23ade-6ab6-4fdd-b0fd-d3a295972479"
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {},
      "outputs": [],
      "source": [
        "d <- matrix(nrow = nrow(coords), ncol = nrow(coords))"
      ],
      "id": "2831c56e"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now compute the distances for all pairs of locations."
      ],
      "id": "32c6abfa-547d-441b-9b72-3f8abbae6fc2"
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {},
      "outputs": [],
      "source": [
        "# compute distances\n",
        "for (i in 1:nrow(coords)) {\n",
        "  for (j in 1:nrow(coords)) {\n",
        "    d[i, j] <- sqrt((coords[i, 1] - coords[j, 1])^2 + \n",
        "                      (coords[i, 2] - coords[j, 2])^2)\n",
        "  }\n",
        "}\n",
        "d # print d"
      ],
      "id": "ded7214a"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can compare this to the output of the `dist` function to confirm that we have computed the distances correctly."
      ],
      "id": "e9361477-afad-48d3-ba12-67500bc986cb"
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {},
      "outputs": [],
      "source": [
        "as.matrix(dist(coords))"
      ],
      "id": "8f376768"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## `while` loops ([YouTube](https://youtu.be/U-IEqpNgIrU), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=0a4fa32f-d2a8-424e-8145-b0ad0107ad54))\n",
        "\n",
        "A `while` loop is a loop that repeats a set of actions while a condition is `TRUE`. The loop stops running when the provided condition is `FALSE`.\n",
        "\n",
        "A `while` loop can be used to do the same thing as a `for` loop but is better suited for situations for which we don’t know how many times we need to iterate.\n",
        "\n",
        "A common use of a `while` loop is to iteratively improve an approximation one iteration at a time until the change in the approximation is smaller than some standard of precision.\n",
        "\n",
        "The basic syntax of a `while` loop is shown in the psuedo code below.\n",
        "\n",
        "    while (condition) {\n",
        "      # some code to run\n",
        "    }\n",
        "\n",
        "We have the following components of the `while` loop:\n",
        "\n",
        "-   The loop starts with `while`, which indicates that we want to implement a `while` loop.\n",
        "-   Inside the parentheses after `while`, we have a `condition`.\n",
        "-   If `condition` is `TRUE`, we run the code in between the curly braces.\n",
        "-   We then return to the `condition`. If `condition` is still `TRUE`, then we run the code in between the curly praces again.\n",
        "-   We keep repeating the previous step until `condition` is `FALSE`.\n",
        "\n",
        "We can read a `while` loop in the following way: While our condition is `TRUE`, run the following code.\n",
        "\n",
        "Here’s what happens in a `while` loop:\n",
        "\n",
        "1.  We check whether our condition is `TRUE`.\n",
        "2.  If the condition is `TRUE`, then we run the code in the curly braces. This code will perform some actions that should potentially have an impact on the condition.\n",
        "3.  We repeat steps 1 and 2 until our condition is `FALSE`.\n",
        "4.  We exit the loop.\n",
        "\n",
        "We provide a trivial example of a `while` loop that prints the current value of an iterator, `i`.\n",
        "\n",
        "-   We initially set `i` to `1`.\n",
        "-   The `while` loop runs as long as `i` is less than 4.\n",
        "-   We print the current value of `i` and then increment `i` (meaning we increase it by 1). - We repeat printing the current value of `i` until `i` is 4."
      ],
      "id": "ae05e7d0-5af9-4956-a8d2-47478f7df9e1"
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {},
      "outputs": [],
      "source": [
        "# set iterator\n",
        "i <- 1\n",
        "while (i < 4) {\n",
        "  message(\"The current iterator value is \", i)\n",
        "  i <- i + 1\n",
        "}"
      ],
      "id": "f57ac522"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`while` loops tend to be more dangerous than `for` loops because it is easy create an “infinite loop” \\`that never stops running.\n",
        "\n",
        "An infinite loop occurs when our `while` condition is always satisfied, so the loop keeps running indefinitely or until an error occurs (like running out of memory).\n",
        "\n",
        "If we want to crash our computer, we can run the following code:\n",
        "\n",
        "    i <- 2\n",
        "    while (i > 0) {\n",
        "      i = i^2\n",
        "    }\n",
        "\n",
        "In the code above, we start by setting `i` equal to `2`. We then check whether `i` is more than 0. If it is more than 0 then we update `i` to be the current value of `i` squared. Our condition will always be `TRUE`, so we keep squaring the current value of `i` indefinitely.\n",
        "\n",
        "If we are stuck in an infinite loop, then:\n",
        "\n",
        "-   Pressing `Ctrl + c` (PC) or `Cmd + c` (Mac) may stop a loop.\n",
        "-   Sometimes, a “STOP” button is available in the Console. Clicking this button will attempt to stop a loop.\n",
        "-   In R Studio, we can click “Session” in the Menu and then click “Interrupt R” to stop the code being executed.\n",
        "\n",
        "### Bisection example ([YouTube](https://youtu.be/my39zmDzUe4), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=37981baa-776b-4e1b-bb28-b0ad01087114))\n",
        "\n",
        "The bisection method is a method for finding the “zeros” or “roots” of a continuous function, i.e., where the function equals 0.\n",
        "\n",
        "The bisection method requires that we start the algorithm at two values such that the function has opposite signs for those values.\n",
        "\n",
        "Suppose we have a continuous function, $f(x)$.\n",
        "\n",
        "The bisection method is implemented in the following way:\n",
        "\n",
        "1.  Find two points, $a$ and $b$, such that $a < b$ and $f(a)f(b) < 0$.\n",
        "2.  Compute the average of $a$ and $b$, $m = (a + b)/2$\n",
        "3.  Determine whether $| f(m) | > \\epsilon$, where $\\epsilon$ is some very small positive number like 0.0001. If our condition is satisfied, perform the following steps:\n",
        "    1.  If $f(m)f(a) < 0$, then a root exists between $a$ and $m$. Set $b = m$ and return to step 2.\n",
        "    2.  Otherwise, a root exists between $m$ and $b$. Set $a = m$ and return to step 2.\n",
        "\n",
        "When our condition is NOT satisfied, we have found an (approximate) root of our function.\n",
        "\n",
        "We now create a function, `bisect`, that implements the bisection method. The function will take arguments:\n",
        "\n",
        "-   `lower`: The smaller starting value used to implement the bisection method.\n",
        "-   `upper`: The upper starting value used to implement the bisection method.\n",
        "-   `fun`: A continuous function that can take both positive and negative values.\n",
        "-   `eps`: A tolerance value indicating how close we must be to a true “zero”.\n",
        "\n",
        "We implement the bisection method below. We:\n",
        "\n",
        "1.  Start by checking whether `fun(lower) * fun(upper) > 0`. If this is `TRUE`, then then the bisection algorithm may fail and we might end up in an infinite loop!. We use the `stop` function to exit the function with an error if we have bad starting values.\n",
        "2.  Compute `m` using the `lower` and `upper` values.\n",
        "3.  Check if `m` is a root. If it’s not a root:\n",
        "    1.  Use a control statement to decide whether to update the `lower` or `upper` value to `m`.\n",
        "    2.  Compute `m` using the `lower` and `upper` values.\n",
        "    3.  Return so step 3.\n",
        "4.  Return the final value of `m`, which is an (approximate) root of `fun`."
      ],
      "id": "a61da19b-2525-41a6-9fa8-a16780fc9c1a"
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {},
      "outputs": [],
      "source": [
        "bisect <- function(lower, upper, fun, eps = 0.0001) {\n",
        "  # make sure our initial lower and upper values produce\n",
        "  # opposite signs when plugged into fun\n",
        "  if (fun(lower) * fun(upper) > 0) {\n",
        "    stop(\"\\nBig yikes!\",\n",
        "         \"\\nThose starting values are sus.\",\n",
        "         \"\\nYou need to glow up those values so  that fun(lower) * fun(upper) < 0.\")\n",
        "  }\n",
        "  # compute m\n",
        "  m <- (lower + upper)/2\n",
        "  # determine if we're at (or really close to a root)\n",
        "  while (abs(fun(m)) > eps) { # if we're not close to a root\n",
        "    # determine whether a root is below m or above m\n",
        "    if (fun(m) * fun(lower) < 0) { # if the root is below m\n",
        "      upper <- m\n",
        "    } else { # if the root is above m\n",
        "      lower <- m\n",
        "    }\n",
        "    # determine a new candidate for m\n",
        "    m <- (lower + upper)/2\n",
        "  }\n",
        "  return(m)\n",
        "}"
      ],
      "id": "78269381"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Consider the continuous polynomial function $$\n",
        "f(x) = (x - 3)(x + 1).\n",
        "$$\n",
        "\n",
        "From looking at this equation, we can easily verify that $f(x) = 0$ when $x = 3$ or $x = -1$.\n",
        "\n",
        "We create a function for $f(x)$ below."
      ],
      "id": "078a9788-6239-475a-93b3-9dd8334d292b"
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "metadata": {},
      "outputs": [],
      "source": [
        "f <- function(x) {\n",
        "  (x - 3) * (x + 1)\n",
        "}"
      ],
      "id": "be08a5b3"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now plot $f(x)$ for for a sequence of $x$ values between -2 and 4."
      ],
      "id": "54b682e5-250a-4745-a8e6-2b946d53ea08"
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "metadata": {},
      "outputs": [],
      "source": [
        "xs <- seq(-2, 4, len = 1000)\n",
        "plot(xs, f(xs), type = \"l\", xlab = \"x\", ylab = \"f(x)\")"
      ],
      "id": "a3c0c1c4"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will try to find a root of `f` using `bisect` using starting values of `2` and `5`."
      ],
      "id": "5aeca70f-c1c4-422d-9510-f2c4d7dc37b8"
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "metadata": {},
      "outputs": [],
      "source": [
        "bisect(lower = 2, upper = 5, fun = f)"
      ],
      "id": "ecdbd3f3"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If we choose different starting values then we can find the other root of $f$."
      ],
      "id": "17aa23fa-ac30-4756-9fe1-296cdce10f7e"
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "metadata": {},
      "outputs": [],
      "source": [
        "bisect(lower = -2, upper = 1, fun = f)"
      ],
      "id": "70039122"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What happens if we choose bad starting values?"
      ],
      "id": "a7abd105-152c-4549-a91b-bcf6076fcbf7"
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "metadata": {},
      "outputs": [],
      "source": [
        "bisect(lower = -4, upper = 5, fun = f)"
      ],
      "id": "b78ed902"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Loops versus vectorization ([YouTube](https://youtu.be/0AlZW5ZGvvg), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=4bf55aac-da5d-4824-92f2-b0af01266aec))\n",
        "\n",
        "Loops in R are notoriously slow.\n",
        "\n",
        "When possible, we should prefer *vectorized operations* over loops.\n",
        "\n",
        "A **vectorized** operation uses a fast, pre-compiled function to apply an operation to all elements of a vector, matrix, list, etc., without using a loop.\n",
        "\n",
        "Here is a trivial example demonstrating the the speed and power of vectorization over loops.\n",
        "\n",
        "Let’s say we want to compute the cumulative sum of a set of numeric values. The cumulative sum of a vector of numeric values is the sum of each value and all the previous values in the vector. The cumulative sum of the vector $[1, 2, 4]$ is $[1, 3, 7]$.\n",
        "\n",
        "We create a function, `cumsum_loop`, to compute the cumulative sum of a numeric vector. We describe our approach:\n",
        "\n",
        "-   The function takes, `x`, which we assume is a numeric vector.\n",
        "-   We create a vector, `results`, to store our cumulative sums.\n",
        "-   We create a `for` loop to iterate `i` over every element of `x`. For each value of `i`, we:\n",
        "    -   Check if we are considering the first element of `x`. If we are, then `results[i]` is simply the value in the first position of `x`.\n",
        "    -   If we are not considering the first element of `x`, then the cumulative sum of the first `i` elements of `x` is the value in the $i$th position of `x`, `x[i]`, plus the previous value of the cumulative sum, `results[i - 1]`.\n",
        "-   After we iterating over all positions of `x`, we return the `results` vector."
      ],
      "id": "904970c3-8c7c-43ff-bba4-b6b248670fc4"
    },
    {
      "cell_type": "code",
      "execution_count": 40,
      "metadata": {},
      "outputs": [],
      "source": [
        "cumsum_loop <- function(x) {\n",
        "  # create vector of results\n",
        "  results <- numeric(length(x))\n",
        "  for (i in 1:length(x)) {\n",
        "    if (i == 1) {\n",
        "      results[i] <- x[i]\n",
        "    } else {\n",
        "      results[i] <- results[i - 1] + x[i]\n",
        "    }\n",
        "  }\n",
        "  return(results)\n",
        "}"
      ],
      "id": "50ba41e5"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let’s try this function on the vector $[1, 2, 4]$."
      ],
      "id": "da785a9e-f7bd-4724-9885-860445f224df"
    },
    {
      "cell_type": "code",
      "execution_count": 41,
      "metadata": {},
      "outputs": [],
      "source": [
        "cumsum_loop(c(1, 2, 4))"
      ],
      "id": "36844aa5"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What is wrong with this function? Nothing, really. Except that it’s slow compared to the equivalent vectorized version of this calculation.\n",
        "\n",
        "The `cumsum` function is a vectorized version of a cumulative sum that relies on fast, pre-compiled code. We can see that it gives the same results as the function we created."
      ],
      "id": "af636b71-204b-4afe-9192-a3d21a1fb72a"
    },
    {
      "cell_type": "code",
      "execution_count": 42,
      "metadata": {},
      "outputs": [],
      "source": [
        "cumsum(c(1, 2, 4))"
      ],
      "id": "57a9fa6b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use the `microbenchmark` function in the **microbenchmark** package (Mersmann 2023) to compare the timings of `cumsum_loop` and `cumsum`. The `microbenchmark` function will run the provided functions 100 times, by default.\n",
        "\n",
        "-   We sample 1000 values from a standard normal distribution.,\n",
        "-   We use `microbenchmark` to time how long it takes to compute the cumulative sum of the sampled values using both cumulative sum functions. We repeat this timing 100 times for each function.\n",
        "-   We plot the results for comparison."
      ],
      "id": "354e8666-4bd2-4da5-af00-c73be9ece2d1"
    },
    {
      "cell_type": "code",
      "execution_count": 43,
      "metadata": {},
      "outputs": [],
      "source": [
        "# attach microbenchmark\n",
        "library(microbenchmark) \n",
        "# draw 1000 values from a N(0, 1) distribution\n",
        "x <- rnorm(1000)\n",
        "# compute the cumulative sum 100 times using both approaches\n",
        "cumsum_timings <- microbenchmark(cumsum_loop(x), cumsum(x))\n",
        "# plot timing results\n",
        "plot(cumsum_timings)"
      ],
      "id": "40db466f"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "While the exact timing results will depend on the computer the timings are performed on, the plot above will clearly show that the vectorized `cumsum` function is substantially faster than the `cumsum_loop` function.\n",
        "\n",
        "These results may seem trivial, but for large data sets where these operations are performed thousands or millions of times, vectorized operations can be many orders of magnitude faster than loops.\n",
        "\n",
        "# The vectorized \\*apply functions ([YouTube](https://youtu.be/nl-ktoJriH8), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=92d3b13f-35c1-4829-bd79-b0ad010b20ed))\n",
        "\n",
        "The `*apply` functions are a set of vectorized functions that can be applied different object types. The “`*`” in `*apply` is a “wildcard” character, meaning that it can be any set of characters.\n",
        "\n",
        "`*apply` functions will typically perform operations much faster than equivalent operations implemented using loops. `*apply` functions also make code simpler and more elegant.\n",
        "\n",
        "[Table 1](#tbl-apply) summarizes the various `*apply` functions available by default in R.\n",
        "\n",
        "| name     | description                                                                                |\n",
        "|:-----|:-----------------------------------------------------------------|\n",
        "| `apply`  | Applies a function to the margins of a data frame, matrix, or array.                       |\n",
        "| `lapply` | Applies a function to each element of a list.                                              |\n",
        "| `sapply` | A version of `lapply` that attempts to simplify the results to a vector, matrix, or array. |\n",
        "| `vapply` | A safer version of `sapply` that returns a pre-specified type of return value.             |\n",
        "| `tapply` | Applies a function to a vector or data frame based on the levels of a categorical variable |\n",
        "| `mapply` | A multivariate version of `sapply`.                                                        |\n",
        "\n",
        "Table 1: A summary of the various `*apply` functions available in R.\n",
        "\n",
        "### The `apply` function ([YouTube](https://youtu.be/iWq_EsrXhjQ), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=98faa071-9645-4c1c-afe8-b0ad010c0e64))\n",
        "\n",
        "The `apply` function applies a function to a margin (dimension) of a data frame, matrix, or array.\n",
        "\n",
        "The primary arguments of the `apply` function are:\n",
        "\n",
        "-   `X`: A data frame, `matrix`, or `array`.\n",
        "-   `MARGIN`: a vector indicating the dimension to apply the function to. `1` indicates rows, `2` indicates columns, `3` indicates the third dimension of an array, `1:2` indicates row and columns, etc.\n",
        "-   `FUN`: the function to be applied to the margins of `X`.\n",
        "-   `...`: additional arguments that need to be passed to `FUN`.\n",
        "-   `simplify`: a logical value indicating whether results should be simplified, if possible. The default is `TRUE`.\n",
        "\n",
        "We provide several examples using the `apply` function below.\n",
        "\n",
        "Consider the matrix `mat` below."
      ],
      "id": "c33496b9-3c8a-45b3-a7b1-df7a33fe4320"
    },
    {
      "cell_type": "code",
      "execution_count": 45,
      "metadata": {},
      "outputs": [],
      "source": [
        "mat <- cbind(1:5, 6:10)\n",
        "mat"
      ],
      "id": "b0b333f6"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can compute the mean value for each row of `mat` by running the command below."
      ],
      "id": "25048e70-b7ec-4536-a36f-08702ffe61c1"
    },
    {
      "cell_type": "code",
      "execution_count": 46,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(X = mat, MARGIN = 1, FUN = mean)"
      ],
      "id": "7c74ca40"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can compute the mean value for each column of `mat` by running the command below."
      ],
      "id": "33fdd8a5-552d-46f3-858f-29c5d2a02422"
    },
    {
      "cell_type": "code",
      "execution_count": 47,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(mat, 2, mean)"
      ],
      "id": "2a862a36"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that the `rowMeans` and `colMeans` functions would have produced the same results, with simpler code. However, the examples above are simple examples to illustrate usage of the `apply` function.\n",
        "\n",
        "What if we want to apply a function that has arguments we need to specify? E.g., let’s say we want to compute the 0.8 quantile of a set of values in a numeric vector `x`. We could run the command `quantile(x, prob = 0.8)`. How can we use a function like this with the `apply` function?\n",
        "\n",
        "The code below uses `apply` to compute the 0.8 quantile for each column of `mat`."
      ],
      "id": "7ca779b9-cd1b-4796-a180-21964ec54432"
    },
    {
      "cell_type": "code",
      "execution_count": 48,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(mat, 2, quantile, prob = 0.8)"
      ],
      "id": "fda86a62"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lastly, we consider a basic application of the `apply` function a 3-dimensional array.\n",
        "\n",
        "We create a 3-dimensional array with dimension $3\\times 2 \\times 4$ below."
      ],
      "id": "e548e839-80d4-4ef7-b3dd-68f62d139fbc"
    },
    {
      "cell_type": "code",
      "execution_count": 49,
      "metadata": {},
      "outputs": [],
      "source": [
        "A <- array(1:24, dim = c(3, 2, 4))\n",
        "print(A)"
      ],
      "id": "b5b06a51"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let’s compute means for each row of `A`."
      ],
      "id": "bf7a5971-a60b-4010-92ce-c1b12357708a"
    },
    {
      "cell_type": "code",
      "execution_count": 50,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(A, 1, mean)"
      ],
      "id": "29d86758"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`11.5` comes from averaging all the values in the 1st row of `A` across all columns and the third dimension: `(1 + 4 + 7 + 10 + 13 + 16 + 19 + 22)/8 = 11.5`. More simply, it is equivalent to running the following command."
      ],
      "id": "902ba1d1-34bd-4626-954e-64b08493ae03"
    },
    {
      "cell_type": "code",
      "execution_count": 51,
      "metadata": {},
      "outputs": [],
      "source": [
        "mean(A[1,,])"
      ],
      "id": "6b391191"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`12.5` and `13.5` come from averaging all the values in the second and third rows of `A`, respectively.\n",
        "\n",
        "Let’s compute the sums across the third dimension of `A`."
      ],
      "id": "80faebea-32ce-4eb2-9601-b0ed03d5e7c6"
    },
    {
      "cell_type": "code",
      "execution_count": 52,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(A, 3, sum)"
      ],
      "id": "f3023fe0"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`21` comes from adding all the values in the first element of the third dimension of `A`. More specifically, it is equivalent to running the following command."
      ],
      "id": "c98d6ad3-7c11-4f38-a13f-e5069a004855"
    },
    {
      "cell_type": "code",
      "execution_count": 53,
      "metadata": {},
      "outputs": [],
      "source": [
        "sum(A[,,1])"
      ],
      "id": "d1078b3d"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let’s compute the means across the row and columns of `A`."
      ],
      "id": "b9fc200e-baa7-4485-b760-846a77565715"
    },
    {
      "cell_type": "code",
      "execution_count": 54,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(A, 1:2, mean)"
      ],
      "id": "a5372ec1"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The command above returns the values obtained by averaging all of the values in a particular row and column of `A` across the third dimension of `A`.\n",
        "\n",
        "`10` comes from averaging the values in `A[1, 1, 1]`, `A[1, 1, 2]`, `A[1, 1, 3]`, `A[1, 1, 4]`, i.e., `(1 + 7 + 13 + 19)/4 = 10`. Alternatively, it is equivalent to running the following command."
      ],
      "id": "e3abb159-22fa-47a4-8d53-432c1247a4eb"
    },
    {
      "cell_type": "code",
      "execution_count": 55,
      "metadata": {},
      "outputs": [],
      "source": [
        "mean(A[1, 1, ])"
      ],
      "id": "0409b13b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The `lapply`, `sapply`, and `vapply` functions ([YouTube](https://youtu.be/4RjlOShPeU0), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=53c2f52e-8419-4846-83f7-b0ad010c9c25))\n",
        "\n",
        "The `lapply` function applies a function to each element of an atomic vector or list. It is very convenient for performing analyses in simulation studies. The `lapply` function always returns a list.\n",
        "\n",
        "The arguments of the `lapply` function are:\n",
        "\n",
        "-   `X`: An atomic vector or list.\n",
        "-   `FUN`: the function to be applied to each element of `X`.\n",
        "-   `...`: additional arguments that need to be passed to `FUN`.\n",
        "\n",
        "We start with a simple example that computes the sample mean of a `list` containing numeric vectors.\n",
        "\n",
        "We first define the our `list`."
      ],
      "id": "d431458b-f423-458e-a385-e54de961b99a"
    },
    {
      "cell_type": "code",
      "execution_count": 56,
      "metadata": {},
      "outputs": [],
      "source": [
        "simple_list <- list(1:10, 11:20)\n",
        "simple_list"
      ],
      "id": "3e173152"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We then use `lapply` to compute the sample mean of each element of `sample_list`."
      ],
      "id": "66613610-c0fc-43fa-971f-748873fdb7ed"
    },
    {
      "cell_type": "code",
      "execution_count": 57,
      "metadata": {},
      "outputs": [],
      "source": [
        "lapply(simple_list, mean)"
      ],
      "id": "7d911f27"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `sapply` function attempts to coerce the results of `lapply` to a `vector`, `matrix`, or `array` (whichever is simplest) using the `simplify2array` function. Technically, the `sapply` function has two additional arguments compared to `lapply`:\n",
        "\n",
        "-   `simplify`: a logical value indicating whether the `sapply` should attempt to simplify the result using `simplify2array`. The default is `TRUE`.\n",
        "-   `USE.NAMES`: a logical value indicating whether names should be associated with the results. The default is `TRUE`.\n",
        "\n",
        "We repeat the previous example using `sapply`."
      ],
      "id": "f72e19cf-fbc7-4a5e-badc-3b26c2d537b9"
    },
    {
      "cell_type": "code",
      "execution_count": 58,
      "metadata": {},
      "outputs": [],
      "source": [
        "sapply(simple_list, mean)"
      ],
      "id": "da448a34"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Instead of returning a `list` object with elements `5.5` and `15.5`, we get a vector with the same elements.\n",
        "\n",
        "The `sapply` function *seems* like it is the more convenient version of `lapply`, but this is a trap. The `sapply` function can produce unexpected results because it is not required to return an object of a specific type. The result may change depending on the input.\n",
        "\n",
        "Let’s create a second list, `simple_list2`, with two elements: the first element is a vector of 3 numeric values and the second element is a vector of 4 numeric values."
      ],
      "id": "eeff4e79-091c-4944-9938-ec661d68d0f1"
    },
    {
      "cell_type": "code",
      "execution_count": 59,
      "metadata": {},
      "outputs": [],
      "source": [
        "simple_list2 <- list(1:3, 4:7)\n",
        "simple_list2"
      ],
      "id": "c1f1b1d6"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create a function, `add2`, that takes a `numeric` vector `x` and adds 2 to all the elements in `x`."
      ],
      "id": "8a9e4b8a-2b74-4bab-9a44-7d2ee673e922"
    },
    {
      "cell_type": "code",
      "execution_count": 60,
      "metadata": {},
      "outputs": [],
      "source": [
        "add2 <- function(x) {\n",
        "  x + 2\n",
        "}"
      ],
      "id": "c8b39f63"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We confirm that `add2` works as intended."
      ],
      "id": "53a20200-5335-4005-9483-d8e853831497"
    },
    {
      "cell_type": "code",
      "execution_count": 61,
      "metadata": {},
      "outputs": [],
      "source": [
        "add2(2:3)"
      ],
      "id": "dc68463e"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If we use `sapply` to apply `add2` to each element of `simple_list`, the result is returned as a `matrix` because the results can be simplified by `simplify2array`."
      ],
      "id": "370cf267-6f05-4268-bb38-66aad22930b8"
    },
    {
      "cell_type": "code",
      "execution_count": 62,
      "metadata": {},
      "outputs": [],
      "source": [
        "sapply(simple_list, add2)"
      ],
      "id": "4cdd870b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If we apply `add2` to each element of `simple_list2`, the result is returned as a `list` because the results cannot be simplified by `simplify2array`."
      ],
      "id": "68652e9d-e42c-4568-9a1b-a34fc6a149f9"
    },
    {
      "cell_type": "code",
      "execution_count": 63,
      "metadata": {},
      "outputs": [],
      "source": [
        "sapply(simple_list2, add2)"
      ],
      "id": "8deadf43"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We don’t always know that data type the `sapply` function will return, which is dangerous if we use those results in any future calculations.\n",
        "\n",
        "The `vapply` function attempts to be a safer version of `sapply` by requiring the user to indicate the data type expected from each application of `FUN`. However, `vapply` also attempts to simplify the output using `simplify2array`, so it can also create issues in later code.\n",
        "\n",
        "The `vapply` function essentially adds one additional argument to `sapply`.\n",
        "\n",
        "-   `FUN.VALUE`: a template for the value returned by `FUN`.\n",
        "\n",
        "If `FUN` doesn’t have the same structure as `FUN.VALUE`, then `vapply` returns an error.\n",
        "\n",
        "Let’s create a trivial example. We create two lists, `l1` and `l2`.\n",
        "\n",
        "-   `l1` has two elements containing `integer` vectors.\n",
        "-   `l2` has two elements containing `double` vectors."
      ],
      "id": "c7256ebc-4a96-4962-8479-c8529e071b84"
    },
    {
      "cell_type": "code",
      "execution_count": 64,
      "metadata": {},
      "outputs": [],
      "source": [
        "l1 <- list(seq.int(1, 4), seq.int(4, 7))\n",
        "l1"
      ],
      "id": "31c170df"
    },
    {
      "cell_type": "code",
      "execution_count": 65,
      "metadata": {},
      "outputs": [],
      "source": [
        "l2 <- list(seq(1.5, 4.5, by = 1), seq(4.5, 7.5, by = 1))\n",
        "l2"
      ],
      "id": "3970cf84"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We apply the `max` function to `l1` using `vapply`. We specify `FUN.VALUE = integer(1)` to indicate that each time we apply the `max` function to an element of `l1`, we expect to get an `integer` vector of length 1."
      ],
      "id": "ba888d92-79f7-4d3a-b914-fc7a32f90faf"
    },
    {
      "cell_type": "code",
      "execution_count": 66,
      "metadata": {},
      "outputs": [],
      "source": [
        "vapply(l1, max, FUN.VALUE = integer(1))"
      ],
      "id": "4d31a93e"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We get an error if we do the same thing for `l2` because the `max` function will return an atomic vector of type `double` of length 1. The error produced by `vapply` indicates that it was expecting to get an `integer` type by got a `double` type instead."
      ],
      "id": "0236ba7d-2101-48ea-af6b-632c404b50d3"
    },
    {
      "cell_type": "code",
      "execution_count": 67,
      "metadata": {},
      "outputs": [],
      "source": [
        "vapply(l2, max, FUN.VALUE = integer(1))"
      ],
      "id": "62a73323"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can correct the problem by changing `FUN.VALUE` to `double(1)`."
      ],
      "id": "ee990015-f431-42a5-9497-96a9525d6ed3"
    },
    {
      "cell_type": "code",
      "execution_count": 68,
      "metadata": {},
      "outputs": [],
      "source": [
        "vapply(l2, max, FUN.VALUE = double(1))"
      ],
      "id": "17457238"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We should not use `sapply` or `vapply` inside any function to avoid problems with the output change data types. Instead, use the `lapply` function to keep the output consistent.\n",
        "\n",
        "### Bootstrap confidence interval example ([YouTube](https://youtu.be/IBuJSRu9Quw), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=88db6193-d823-4c75-a00b-b0ad010d3fb3))\n",
        "\n",
        "A confidence interval procedure produces plausible values for a target parameter with an associated confidence level.\n",
        "\n",
        "-   If we repeatedly obtain samples from a population, a confidence interval procedure will produce intervals containing the target parameter with probability equal to the confidence level.\n",
        "\n",
        "Confidence intervals are typically constructed on the assumption that we know the population distribution from which the sample came.\n",
        "\n",
        "A bootstrap confidence interval procedure can be used to construct confidence intervals when the population distribution is unknown.\n",
        "\n",
        "More generally, the bootstrap method allows us to approximate the sampling distribution of a statistic by using the observed data to produce simulated data sets.\n",
        "\n",
        "Suppose we have a sample of $n$ data values, $x_1, x_2, \\ldots, x_n$. Suppose we want to estimate the mean, $\\mu$, of the population from which the sample came. A natural estimator of $\\mu$ is the sample mean, $\\bar{x}$.\n",
        "\n",
        "The procedure for producing a bootstrap confidence interval for $\\mu$ with confidence level $1-\\alpha$ is:\n",
        "\n",
        "1.  Sample with replacement a data set of size $n$ from the observed data values. This is a bootstrap sample.\n",
        "2.  Compute $\\bar{x}$ for the random sample.\n",
        "3.  Perform steps 1 and 2 $B$ times, where $B$ is a number like 1000 or 10000.\n",
        "4.  Compute the $\\alpha/2$ and $1-\\alpha/2$ quantiles of the sample means obtained from the bootstrap samples.\n",
        "\n",
        "Suppose we have observed the following sample:\n",
        "\n",
        "3.318, 4.940, 5.516, 2.958, 6.098, 6.858, -0.539, 2.544, 4.091, 1.884\n",
        "\n",
        "Let’s construct a 95% bootstrap confidence interval for the population mean, $\\mu$, using the sample data. We will build a function to construct our bootstrap confidence interval that is based on `lapply`. Our function will take three arguments:\n",
        "\n",
        "-   `x`: A vector of numeric values.\n",
        "-   `conf.level`: The confidence level desired for the interval. The default value will be `0.95`.\n",
        "-   `B`: The number of bootstrap samples from which to build the confidence interval. The default value will be `1000`.\n",
        "\n",
        "The function will perform the following actions:\n",
        "\n",
        "-   Determine `alpha` from `confidence.level`.\n",
        "-   The `lapply` function is applied to each element of the atomic vector `1:B`.\n",
        "-   We supply an *anonymous* function to the `FUN` argument of `lapply`. An anonymous function is a temporary function we create to perform an action but for which we don’t want to go to the trouble of storing in memory for long-term use. In our case, the anonymous function takes the argument `i`, which comes from the elements of `1:B`.\n",
        "-   Inside our anonymous function, we sample `n` values with replacement from `x` using the `sample` function, making sure to set the `replace` argument of `sample` to `TRUE`. This produces a bootstrap sample.\n",
        "-   We compute the sample mean of our bootstrap sample and return that value from our anonymous function.\n",
        "-   We store the $B$ sample means produced by the `lapply` function in the `bstats` object.\n",
        "-   We use the `unlist` function to convert the list of sample means in `bstats` to a `numeric` vector.\n",
        "-   We use the `quantile` function to return the `alpha/2` and `1-alpha/2` quantiles from our vector of sample means."
      ],
      "id": "d34bd139-543e-485b-b854-be157ce70c07"
    },
    {
      "cell_type": "code",
      "execution_count": 69,
      "metadata": {},
      "outputs": [],
      "source": [
        "bootstrap_ci <- function(x, conf.level = 0.95, B = 1000) {\n",
        "  # determine alpha\n",
        "  alpha <- 1 - conf.level\n",
        "  # compute sample mean of B bootstrap samples\n",
        "  bstats <- lapply(1:B, FUN = function(i) {\n",
        "    # draw bootstrap sample\n",
        "    boot_sample <- sample(x, size = length(x), replace = TRUE)\n",
        "    # compute sample mean of bootstrap sample\n",
        "    mean(boot_sample)\n",
        "  })\n",
        "  # convert list to vector\n",
        "  bstats <- unlist(bstats)\n",
        "  # compute desired quantiles for bootstrap ci\n",
        "  quantile(bstats, prob = c(alpha/2, 1-alpha/2))\n",
        "}"
      ],
      "id": "b10edc14"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now apply our function to the observed sample of data."
      ],
      "id": "f5048633-f3c2-48e1-86fd-df3b4eefd8cc"
    },
    {
      "cell_type": "code",
      "execution_count": 70,
      "metadata": {},
      "outputs": [],
      "source": [
        "set.seed(43)\n",
        "x <- c(3.318, 4.940, 5.516, 2.958, 6.098, 6.858, -0.539, 2.544, 4.091, 1.884)\n",
        "bootstrap_ci(x)"
      ],
      "id": "0ccd49ae"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Our 95% bootstrap confidence interval for the population mean is \\[2.44, 4.98\\]. The sample was a random sample from a normal distribution with mean 3 and standard deviation 2, so our interval (in this case) did contain the true value.\n",
        "\n",
        "## The `tapply` function ([YouTube](https://youtu.be/CLLK9a8d_jA), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=1058278d-dbb8-4edd-b2ef-b0ad010de1b8))\n",
        "\n",
        "The `tapply` function is used to apply a function to all the values associated with a certain level of a `factor` variable. It is useful for summarizing information for different groups of data.\n",
        "\n",
        "The `tapply` function takes the following arguments:\n",
        "\n",
        "-   `X`: an object for which a `split` method exists. Usually, this is a `vector` or `data.frame` object.\n",
        "-   `INDEX`: a `factor` variable with the same length as `X`. This will be used to “split” `X` into different groups.\n",
        "-   `FUN`: The function to be applied to each group of `X` after the `split` function is used.\n",
        "-   `...`: additional arguments that will be passed to `FUN`.\n",
        "-   `default`: an argument we can probably ignore. Run `?tapply` in the Console for more details about this argument.\n",
        "-   `simplify`: a logical value indicating whether the result should be simplified to a scalar array. The default is `TRUE`. Otherwise, a `list` is always returned.\n",
        "\n",
        "We will use the `tapply` function to summarize chicken growth across different feed types.\n",
        "\n",
        "The `chickwts` data set in the **datasets** package is a data frame with 71 observations of 2 variables. The data provide the results of an experiment to measure and compare the effectiveness of various feed supplements on the growth rate of chickens. The following variables are included in the data frame:\n",
        "\n",
        "-   `weight`: the chick’s weight (g).\n",
        "-   `feed`: a `factor` indicating the feed type.\n",
        "\n",
        "We use `tapply` to compute the average chicken weight of the chicken’s receiving each feed type."
      ],
      "id": "9cfcf2f0-1436-44d1-b499-cd82bde76518"
    },
    {
      "cell_type": "code",
      "execution_count": 71,
      "metadata": {},
      "outputs": [],
      "source": [
        "tapply(X = chickwts$weight, INDEX = chickwts$feed, FUN = mean)"
      ],
      "id": "f0d1fe59"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The average weight of the chickens eating the `casein` feed type was 323.6 g. The average weight of the chickens eating the `soybean` feed type was 246.4 g.\n",
        "\n",
        "We can use the `tapply` function to compute a more complex summary of each group’s weights using the `summary` function."
      ],
      "id": "d571ad44-9a8b-432f-bd5d-5d7642d09411"
    },
    {
      "cell_type": "code",
      "execution_count": 72,
      "metadata": {},
      "outputs": [],
      "source": [
        "tapply(X = chickwts$weight, INDEX = chickwts$feed, FUN = summary)"
      ],
      "id": "0b55ca76"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Frequency table example ([YouTube](https://youtu.be/X6kNWXR1hTY), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=c778d1da-1509-481f-8f3d-b0ad010e6cb5))\n",
        "\n",
        "We can use the `tapply` function to construct a frequency table for a `numeric` variable based on a set of bins that break the variable into different intervals. Plotting a frequency table produces a histogram.\n",
        "\n",
        "The `faithful` data set in the **datasets** package summarizes the duration and waiting time between eruptions for the Old Faithful geyser in Yellowstone National Park. The data are a data frame with 272 observations of the variables:\n",
        "\n",
        "-   `eruptions`: eruption duration (minutes).\n",
        "-   `waiting`: waiting time until next eruption (minutes).\n",
        "\n",
        "We will create a frequency table for the `eruptions` variable. For simplicity, we start by extracting the `eruptions` variable from the `faithful` data frame."
      ],
      "id": "538c1002-79a2-4bfb-82ef-48d7289bc643"
    },
    {
      "cell_type": "code",
      "execution_count": 73,
      "metadata": {},
      "outputs": [],
      "source": [
        "eruptions <- faithful$eruptions"
      ],
      "id": "2272368b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we use the `pretty` function to automatically create breakpoints for the `eruptions` variable.\n",
        "\n",
        "The `pretty` function has two primary arguments:\n",
        "\n",
        "-   `x`: a `numeric` object. The range of the object is used to create the breakpoints.\n",
        "-   `n`: the desired number of intervals to be created. This is often ignored by the function.\n",
        "\n",
        "The `pretty` function attempts to produce a set of “pretty breakpoints” based on `x` and `n`. Run `?pretty` in the Console for more details about this function.\n",
        "\n",
        "We pass `eruptions` to the `pretty` function and ask it to produce approximately 10 intervals."
      ],
      "id": "80628f73-f089-4a0e-a0db-a711b0ec1824"
    },
    {
      "cell_type": "code",
      "execution_count": 74,
      "metadata": {},
      "outputs": [],
      "source": [
        "breaks <- pretty(eruptions, n = 10)\n",
        "breaks"
      ],
      "id": "a47ae127"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We end up with 9 breakpoints, which will break our `eruptions` into 8 intervals.\n",
        "\n",
        "To use the `tapply` function to count the number of values in each interval, we need to associate each eruption value with a level based on our breakpoints. Instead of using a complex loop with numerous `if` statements to categorize each eruption value, we can use the `cut` function.\n",
        "\n",
        "The `cut` function takes two primary arguments:\n",
        "\n",
        "-   `x`: a numeric vector that is assigned a `factor` label by “cutting”.\n",
        "-   `breaks`: the breakpoints that define the intervals that “cut” `x` into different groups. Each group is assigned a different label.\n",
        "\n",
        "We pass the `eruptions` variable and the `breaks` variable returned by the `pretty` function to `cut` to associate each eruption value with an interval. We assign these labels the name `interval_labels`. We use the `head` function to view the first values in `interval_labels`."
      ],
      "id": "7e62e2be-8448-4f5c-83d2-c095bccbed59"
    },
    {
      "cell_type": "code",
      "execution_count": 75,
      "metadata": {},
      "outputs": [],
      "source": [
        "interval_labels <- cut(eruptions, breaks)\n",
        "head(interval_labels)"
      ],
      "id": "0cab7037"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The first eruption value was in the interval (3.5, 4), the second in the interval (1.5, 2\\], the third in the interval (3, 3.5\\], etc.\n",
        "\n",
        "Now, we use `tapply` to break the `eruptions` variable into 8 groups based on `interval_labels`, then use the `length` function to determine the number of eruption durations in each interval (which is the number of values observed in that interval)."
      ],
      "id": "3843c68d-7e7b-47be-8386-4cecd8bb7ac3"
    },
    {
      "cell_type": "code",
      "execution_count": 76,
      "metadata": {},
      "outputs": [],
      "source": [
        "count <- tapply(eruptions, interval_labels, length)\n",
        "count"
      ],
      "id": "577b8d30"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "There were 55 eruptions that lasted between 1.5 and 2 minutes. There were 37 eruptions that lasted between 2 and 2.5 minutes, etc.\n",
        "\n",
        "## The `mapply` function ([YouTube](https://youtu.be/kAMAVfNvi8s), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=d0288e67-c616-4b43-9397-b0ad010f8d49))\n",
        "\n",
        "The `mapply` function is intended to be a multivariate version of `sapply`. More specifically, it will vectorize over multiple arguments instead of a single one. The main arguments of the `mapply` function are:\n",
        "\n",
        "-   `FUN`: the function that will be applied to the arguments being vectorized over.\n",
        "-   `...`: the arguments to be vectorized over.\n",
        "-   `MoreArgs`: a `list` with named elements that supply additional arguments need for `FUN`.\n",
        "\n",
        "There are also `SIMPLIFY` and `USE.NAMES` arguments that we do not discuss. Run `?mapply` in the Console for more details.\n",
        "\n",
        "The `mapply` function is more complicated than the other `*apply` functions. I have rarely used it because there are often simpler ways to solve the problem that do not require the complexity of `mapply`.\n",
        "\n",
        "We will walk through a trivial example that uses `mapply`.\n",
        "\n",
        "Suppose we have two vectors of the same size. The first vector contains the initials of several individuals. The second vector contains the initials of their romantic partner."
      ],
      "id": "7ac6ce81-cc83-4aee-8bba-d6ff23e7f10e"
    },
    {
      "cell_type": "code",
      "execution_count": 77,
      "metadata": {},
      "outputs": [],
      "source": [
        "partner1 <- c(\"PM\", \"JI\", \"RX\", \"HP\", \"XE\")"
      ],
      "id": "a2133a6c"
    },
    {
      "cell_type": "code",
      "execution_count": 78,
      "metadata": {},
      "outputs": [],
      "source": [
        "partner2 <- c(\"QN\", \"AA\", \"OJ\", \"DD\", \"GD\")"
      ],
      "id": "743cedf9"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will use `mapply` to `paste` the names together. We set `USE.NAMES` to `FALSE` for better formatting."
      ],
      "id": "d1e3e8f9-b705-4c48-b5be-5732be9cf7d1"
    },
    {
      "cell_type": "code",
      "execution_count": 79,
      "metadata": {},
      "outputs": [],
      "source": [
        "mapply(FUN = paste, partner1, partner2, USE.NAMES = FALSE)"
      ],
      "id": "c65c81c9"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This worked as expected, but the formatting is unpleasant.\n",
        "\n",
        "We can improve the formatting of our result by changing the `sep` argument of `paste` to `\" + \"` or `\" and \"`. We can accomplish this change by passing this information in a named list through the `MoreArgs` argument."
      ],
      "id": "c0b969d4-8b27-4cb8-b003-873d4cf7aceb"
    },
    {
      "cell_type": "code",
      "execution_count": 80,
      "metadata": {},
      "outputs": [],
      "source": [
        "mapply(FUN = paste, partner1, partner2, MoreArgs = list(sep = \" + \"), USE.NAMES = FALSE)"
      ],
      "id": "6cdadded"
    },
    {
      "cell_type": "code",
      "execution_count": 81,
      "metadata": {},
      "outputs": [],
      "source": [
        "mapply(FUN = paste, partner1, partner2, MoreArgs = list(sep = \" and \"), USE.NAMES = FALSE)"
      ],
      "id": "a72470cd"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that we must name any arguments provided after `...` in a function. The reason for this is that, otherwise, the arguments after `...` are assumed to be something intended to be passed through `...`.\n",
        "\n",
        "This is why the code below doesn’t work as intended."
      ],
      "id": "7f66df8b-75da-44f6-a470-09a0147ab20b"
    },
    {
      "cell_type": "code",
      "execution_count": 82,
      "metadata": {},
      "outputs": [],
      "source": [
        "mapply(FUN = paste, partner1, partner2, list(sep = \" and \"), FALSE)"
      ],
      "id": "2beab2c9"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What is actually happening is that `\" and \"` and `\"FALSE\"` are being passed into `paste` with the values from `partner1` and `partner2`, which then pastes the four individual character strings into a single string.\n",
        "\n",
        "## Parallel computing, progress bars, `*apply` functions ([YouTube](https://youtu.be/2LXW8d-lOEM), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=92a4020b-f79c-456d-a91c-b0ad0110622e))\n",
        "\n",
        "We can sometimes speed up the computations performed by an `*apply` function using parallel processing (processing multiple tasks *in parallel* instead of one-at-a-time).\n",
        "\n",
        "Parallel processing is particularly helpful with using `lapply`, `sapply`, or `vapply` to perform complex analysis on a sequence of simulated data sets. The analysis can be performed on each simulated data set separately, which means the operation can be parallelized.\n",
        "\n",
        "We provide small examples to illustrate the advantages of parallel computing when using `*apply` functions.\n",
        "\n",
        "We start by creating a simple function, `sleep_iterator`.\n",
        "\n",
        "-   The function takes the argument `i`, which is not used in the function, so that it is compatible with the `lapply` function.\n",
        "-   Inside the function, we use `Sys.sleep(1)` to pause for one second.\n",
        "-   We have the function return `NULL`.\n",
        "\n",
        "We define `sleep_iterator` below."
      ],
      "id": "ae8894a5-5cbf-41a2-9710-64e0bcedb273"
    },
    {
      "cell_type": "code",
      "execution_count": 83,
      "metadata": {},
      "outputs": [],
      "source": [
        "sleep_iterator <- function(i) {\n",
        "  Sys.sleep(1) # sleep for 1 second\n",
        "  return(NULL) # return NULL\n",
        "}"
      ],
      "id": "7813ba14"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we use the `lapply` function to run `sleep_iterator` 6 times.\n",
        "\n",
        "-   `lapply` will run `sleep_iterator` for `i = 1`, `2`, `...`, `6`.\n",
        "-   We time the operations performed by `lapply` using `system.time`.\n",
        "\n",
        "Unsurprisingly, applying `sleep_iterator` six times results in an elapsed time of about 6 seconds since we pause for 1 second for each of the six values of `i`."
      ],
      "id": "786ace3e-d774-40e0-8fff-36a804ffe107"
    },
    {
      "cell_type": "code",
      "execution_count": 84,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(lapply(1:6, sleep_iterator))"
      ],
      "id": "b65a0ceb"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The **parallel** package (R Core Team 2022) provides parallelized equivalents of several `*apply` functions.\n",
        "\n",
        "-   `parApply` is a parallelized version of `apply`.\n",
        "-   `parLapply` is a parallelized version of `lapply`.\n",
        "-   `parSapply` is a parallelized version of `sapply`.\n",
        "\n",
        "Run `?parallel::clusterApply` to see additional details about these and related parallelized functions.\n",
        "\n",
        "We load the **parallel** package."
      ],
      "id": "fac21c58-35af-4c65-b438-7a7e27ae785c"
    },
    {
      "cell_type": "code",
      "execution_count": 85,
      "metadata": {},
      "outputs": [],
      "source": [
        "library(parallel)"
      ],
      "id": "2fe3dd1f"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The available options and behavior for performing parallel processing depends on our computer’s operating system. In general, it is easier to perform parallel computing on Macs and Linux computers. Windows computers cannot run forked processes through R.\n",
        "\n",
        "We want to use `parLapply` to run `sleep_iterator` 6 times using parallel processing.\n",
        "\n",
        "The primary arguments of `parLapply` are:\n",
        "\n",
        "-   `cl`: a cluster object returned by the `makeCluster` function that defines how the processes will be run in parallel.\n",
        "-   `X`: the atomic vector or list to be iterated over.\n",
        "-   `fun`: the function to be applied to each element of `X`.\n",
        "-   `...`: additional arguments to be passed to `fun`.\n",
        "\n",
        "We use `makeCluster` to setup a parallel processing approach that runs 2 processes in parallel. The code is supposed to be general enough that it will work on Windows, Mac, or Linux computers."
      ],
      "id": "5d0a5b83-290b-4a17-9699-6e11f76256f2"
    },
    {
      "cell_type": "code",
      "execution_count": 86,
      "metadata": {},
      "outputs": [],
      "source": [
        "cl <- makeCluster(getOption(\"cl.cores\", 2))"
      ],
      "id": "cda82af5"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We once again use `system.time` to time how long it takes to run `sleep_iterator` 6 times."
      ],
      "id": "1cefbaef-4f5b-46e9-adc0-9f164d324058"
    },
    {
      "cell_type": "code",
      "execution_count": 87,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(parLapply(cl, 1:6, sleep_iterator))"
      ],
      "id": "055a40f7"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since we are running `sleep_iterator` in parallel using 2 processes, it only takes about 3 seconds to complete our computations instead of the original 6 seconds!\n",
        "\n",
        "The **future.apply** package (Bengtsson 2023b) provides parallelized implementations of several `*apply` functions. The advantage of the **future.apply** package is that code can be created that generalizes from being “evaluated on the local machine, in parallel a set of local machines, or distributed on a mix of local and remote machines” Bengtsson (2021).\n",
        "\n",
        "We start by loading the **future.apply** package."
      ],
      "id": "7e3b2be9-df55-4c68-8c0f-c97faccc8c52"
    },
    {
      "cell_type": "code",
      "execution_count": 88,
      "metadata": {},
      "outputs": [],
      "source": [
        "library(future.apply)"
      ],
      "id": "fa1badd6"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To perform a parallelized operation using **future.apply**, we must use the `plan` function to describe how we want the parallel processing to be performed.\n",
        "\n",
        "The code below sets up a plan to run 2 processes in parallel using a `\"multisession\"` approach. Run `?future::plan` in the Console for more details."
      ],
      "id": "bc678b10-9a25-4b57-a8d2-411dd2b9da2b"
    },
    {
      "cell_type": "code",
      "execution_count": 89,
      "metadata": {},
      "outputs": [],
      "source": [
        "plan(multisession, workers = 2)"
      ],
      "id": "65ed1a6f"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we use the `future_lapply` function to run the `sleep_iterator` function 6 times. The function works identically to the `lapply` function, but uses the parallelization plan to speed up the computations. Similar to before, it only takes about 3 seconds to run the `sleep_iterator` function 6 times because we are running two processes in parallel."
      ],
      "id": "1e4a5e8f-512b-43fa-ba1a-61fb5e41ba36"
    },
    {
      "cell_type": "code",
      "execution_count": 90,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(future_lapply(1:6, FUN = sleep_iterator))"
      ],
      "id": "a14771b7"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "When an `*apply` function takes a long time to return its results, it can be nice to have a progress bar to assess how long the computations will take.\n",
        "\n",
        "The **pbapply** package (Solymos and Zawadzki 2023) creates equivalents of many `*apply` functions that add progress bars showing the time remaining to complete the computation.\n",
        "\n",
        "-   Many of these functions are equipped to use parallel processing through the **parallel** package. We must set the `cl` argument with a `makeCluster` object to use this functionality.\n",
        "-   Run `?pbapply::pbapply` in the Console for a complete list of available functions.\n",
        "\n",
        "We start by loading the **pbapply** package."
      ],
      "id": "ecdcf8a2-589c-419d-9072-c6a4cac7500d"
    },
    {
      "cell_type": "code",
      "execution_count": 91,
      "metadata": {},
      "outputs": [],
      "source": [
        "library(pbapply)"
      ],
      "id": "dc57f30d"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `pblapply` function is an equivalent of `lapply` that includes a progress bar. We can provide a cluster object from the `makeCluster` function to perform the computations in parallel.\n",
        "\n",
        "We use `pblapply` to run `sleep_iterator` 6 times. We use the previously created cluster object, `cl`, which runs 2 processes in parallel. Similar to before, our code runs in about 3 seconds while also showing a progress bar. (We need to run this code interactively to see this progress bar.)"
      ],
      "id": "9d3f3753-0cac-4f29-bdf0-29f30c7d8a0a"
    },
    {
      "cell_type": "code",
      "execution_count": 92,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(pblapply(1:6, FUN = sleep_iterator, cl = cl))"
      ],
      "id": "af05f4de"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# References\n",
        "\n",
        "Bengtsson, Henrik. 2021. “A Unifying Framework for Parallel and Distributed Processing in r Using Futures.” *The R Journal* 13 (2): 208–27. <https://doi.org/10.32614/RJ-2021-048>.\n",
        "\n",
        "———. 2023a. *Future: Unified Parallel and Distributed Processing in r for Everyone*. <https://CRAN.R-project.org/package=future>.\n",
        "\n",
        "———. 2023b. *Future.apply: Apply Function to Elements in Parallel Using Futures*. <https://CRAN.R-project.org/package=future.apply>.\n",
        "\n",
        "Mersmann, Olaf. 2023. *Microbenchmark: Accurate Timing Functions*. <https://github.com/joshuaulrich/microbenchmark/>.\n",
        "\n",
        "R Core Team. 2022. *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing. <https://www.R-project.org/>.\n",
        "\n",
        "Solymos, Peter, and Zygmunt Zawadzki. 2023. *Pbapply: Adding Progress Bar to \\*Apply Functions*. <https://github.com/psolymos/pbapply>."
      ],
      "id": "d1bfb53e-d090-4a4a-83c3-9602354c8f6d"
    }
  ],
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernelspec": {
      "name": "ir",
      "display_name": "R",
      "language": "R"
    },
    "language_info": {
      "name": "R",
      "codemirror_mode": "r",
      "file_extension": ".r",
      "mimetype": "text/x-r-source",
      "pygments_lexer": "r",
      "version": "4.2.2"
    }
  }
}