{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Probability and Distributions in Julia\n",
        "\n",
        "## Overview\n",
        "\n",
        "This tutorial is a reference for the probability you need to do Monte\n",
        "Carlo analysis in this course. It is written for students who have taken\n",
        "a probability or statistics course but may not have used it in a while,\n",
        "and it assumes no prior experience with `Distributions.jl`.\n",
        "\n",
        "Everything here is covered or used in lecture, but lecture moves\n",
        "quickly. Come back to this when you need to look something up.\n",
        "\n",
        "If you want to *plot* distributions, see the “Plotting Distributions”\n",
        "section of [Julia Plotting](julia-plots.qmd). This tutorial is about\n",
        "computing with them.\n",
        "\n",
        "## Setup\n",
        "\n",
        "One package does almost all of the work."
      ],
      "id": "23849c02-1f74-41d0-afa8-96f46f685637"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "using Distributions\n",
        "using Random\n",
        "using Statistics"
      ],
      "id": "2"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Random Variables and Distributions\n",
        "\n",
        "A **random variable** is a quantity whose value we do not know. A\n",
        "**distribution** describes how plausible its possible values are.\n",
        "\n",
        "The distinction that matters in practice is whether the variable is\n",
        "discrete or continuous.\n",
        "\n",
        "### Discrete variables and probability mass\n",
        "\n",
        "A discrete variable takes values you can list: the number of storms in a\n",
        "year, the number of treatment plants that fail. Its distribution is\n",
        "described by a **probability mass function** (PMF), which gives the\n",
        "probability of each individual value:\n",
        "\n",
        "$$p(x) = \\mathbb{P}(X = x)$$\n",
        "\n",
        "These probabilities are genuinely probabilities — they are between 0 and\n",
        "1, and they sum to 1."
      ],
      "id": "5110e646-8e9f-4419-aadc-cbbb70d3fdc5"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.22404180765538775"
            ]
          }
        }
      ],
      "source": [
        "storms = Poisson(3.0)      # average of 3 storms per year\n",
        "pdf(storms, 2)             # probability of exactly 2 storms"
      ],
      "id": "4"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that `Distributions.jl` uses `pdf` for both cases. For a discrete\n",
        "distribution it returns the mass; for a continuous one it returns the\n",
        "density.\n",
        "\n",
        "### Continuous variables and probability density\n",
        "\n",
        "A continuous variable takes any value in a range: a concentration, a\n",
        "flow rate, a temperature. Here the probability of any *exact* value is\n",
        "zero, so we work with a **probability density function** (PDF) instead:\n",
        "\n",
        "$$\\mathbb{P}(a < X < b) = \\int_a^b f(x)\\,dx$$\n",
        "\n",
        "A density is **not** a probability. It can be larger than 1. Only the\n",
        "area under it is a probability."
      ],
      "id": "081e3d75-c603-46b2-a905-4a61773811a8"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.2659615202676218"
            ]
          }
        }
      ],
      "source": [
        "load = Normal(9.0, 1.5)    # CBOD load, mg/L\n",
        "pdf(load, 9.0)             # density at the mean -- not a probability"
      ],
      "id": "6"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Narrow the distribution and the density climbs straight past 1, which no\n",
        "probability can do:"
      ],
      "id": "00a34eac-e3ce-450a-946c-0a21aa5f372e"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "1.9947114020071635"
            ]
          }
        }
      ],
      "source": [
        "pdf(Normal(9.0, 0.2), 9.0)"
      ],
      "id": "8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The Cumulative Distribution Function\n",
        "\n",
        "The **CDF** answers the question you usually actually have: how much\n",
        "probability is at or below some value?\n",
        "\n",
        "$$F(x) = \\mathbb{P}(X \\le x)$$"
      ],
      "id": "fb29e23a-4fed-471a-9982-e332608516ae"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.9087887802741321"
            ]
          }
        }
      ],
      "source": [
        "cdf(load, 11.0)            # probability the load is at or below 11 mg/L"
      ],
      "id": "10"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It works the same way for discrete distributions:"
      ],
      "id": "aedbe935-f3ac-4778-ae3c-6a79c969ff68"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.42319008112684353"
            ]
          }
        }
      ],
      "source": [
        "cdf(storms, 2)             # probability of 2 or fewer storms"
      ],
      "id": "12"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Exceedance probabilities\n",
        "\n",
        "Standards are usually written in terms of being *worse* than a\n",
        "threshold, which is the complement of the CDF:\n",
        "\n",
        "$$\\mathbb{P}(X > x) = 1 - F(x)$$\n",
        "\n",
        "You can write this directly, but `ccdf` (“complementary CDF”) is more\n",
        "accurate for small probabilities, because it avoids subtracting two\n",
        "nearly-equal numbers:"
      ],
      "id": "7f7fea01-bd54-4f56-a319-0bc36030fe06"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.003830380567589775"
            ]
          }
        }
      ],
      "source": [
        "1 - cdf(load, 13.0)        # works"
      ],
      "id": "14"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.003830380567589736"
            ]
          }
        }
      ],
      "source": [
        "ccdf(load, 13.0)           # better"
      ],
      "id": "16"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Quantiles\n",
        "\n",
        "The **quantile function** is the CDF run backwards. Give it a\n",
        "probability, it returns the value that sits there:"
      ],
      "id": "3d9b3e2e-55cb-4aba-aeeb-6081b5a83d97"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "9.0"
            ]
          }
        }
      ],
      "source": [
        "quantile(load, 0.5)        # the median"
      ],
      "id": "18"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "11.467280440427208"
            ]
          }
        }
      ],
      "source": [
        "quantile(load, 0.95)       # the value exceeded 5% of the time"
      ],
      "id": "20"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is how you get the endpoints of an interval. For a 95% interval,\n",
        "take the 2.5% and 97.5% quantiles:"
      ],
      "id": "84b5bc0a-06ee-4d42-a164-8ad7c5848302"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "(6.060054023189911, 11.939945976810087)"
            ]
          }
        }
      ],
      "source": [
        "(quantile(load, 0.025), quantile(load, 0.975))"
      ],
      "id": "22"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Summarizing a Distribution\n",
        "\n",
        "The two summaries you will use constantly:"
      ],
      "id": "50639bcc-929e-4845-9115-a0b91e1ba9ad"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "(9.0, 2.25, 1.5)"
            ]
          }
        }
      ],
      "source": [
        "mean(load), var(load), std(load)"
      ],
      "id": "24"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### The average of a function is not the function of the average\n",
        "\n",
        "This one causes more trouble than anything else on this page.\n",
        "\n",
        "If you push a random variable through a model $h$, the average *output*\n",
        "is generally **not** the output at the average *input*:\n",
        "\n",
        "$$\\mathbb{E}[h(X)] \\neq h(\\mathbb{E}[X])$$\n",
        "\n",
        "They are equal only when $h$ is linear. Here is a quick demonstration\n",
        "with $h(x) = x^2$:"
      ],
      "id": "b5c94ccd-7367-4727-8f18-e18b923b868c"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "(83.25567014845863, 81.0165967503216)"
            ]
          }
        }
      ],
      "source": [
        "Random.seed!(4750)\n",
        "draws = rand(load, 100_000)\n",
        "\n",
        "mean(draws .^ 2), mean(draws)^2"
      ],
      "id": "26"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The two differ by the variance, and for a model as nonlinear as a\n",
        "dissolved oxygen sag curve the gap can be much larger. **This is the\n",
        "reason we sample instead of running a model once at the mean.**\n",
        "\n",
        "## A Short Catalogue\n",
        "\n",
        "You do not need to memorize these. You do need to be able to say why you\n",
        "picked one.\n",
        "\n",
        "| Distribution | Use it when | Constructor |\n",
        "|:-----------------------|:-----------------------|:-----------------------|\n",
        "| Normal | A quantity varies symmetrically about a typical value | `Normal(μ, σ)` |\n",
        "| LogNormal | A positive quantity varies multiplicatively; right-skewed | `LogNormal(μ, σ)` |\n",
        "| Uniform | Every value in a range is equally plausible | `Uniform(a, b)` |\n",
        "| Exponential | Waiting time between independent events | `Exponential(θ)` |\n",
        "| Poisson | Count of independent events in a fixed window | `Poisson(λ)` |\n",
        "| Binomial | Number of successes in a fixed number of trials | `Binomial(n, p)` |\n",
        "| Beta | A proportion, bounded between 0 and 1 | `Beta(α, β)` |\n",
        "| Gamma | A positive, right-skewed quantity such as a sum of waiting times | `Gamma(α, θ)` |\n",
        "\n",
        "Two cautions worth internalizing.\n",
        "\n",
        "**`LogNormal` takes the parameters of the underlying normal**, not the\n",
        "mean and standard deviation of the distribution itself.\n",
        "`LogNormal(log(2), 0.5)` has a median of 2, not a mean of 2.\n",
        "\n",
        "**A `Normal` is unbounded.** If you use it for a concentration, it will\n",
        "eventually hand you a negative one. Either check, or use a distribution\n",
        "that cannot:"
      ],
      "id": "26375220-fcc8-44c2-9cbc-ac2d94bb10b5"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "-1.5881740321384474"
            ]
          }
        }
      ],
      "source": [
        "Random.seed!(4750)\n",
        "minimum(rand(Normal(2.0, 1.0), 10_000))    # negative concentrations"
      ],
      "id": "28"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "Truncated(Distributions.Normal{Float64}(μ=2.0, σ=1.0); lower=0.0)"
            ]
          }
        }
      ],
      "source": [
        "truncated(Normal(2.0, 1.0); lower=0.0)     # one way to fix it"
      ],
      "id": "30"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Drawing Samples\n",
        "\n",
        "`rand` draws from any distribution. One value, or many:"
      ],
      "id": "c4297fdb-d600-45a1-a621-354d787fc203"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "9.741525584287155"
            ]
          }
        }
      ],
      "source": [
        "Random.seed!(4750)\n",
        "rand(load)"
      ],
      "id": "32"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "5-element Vector{Float64}:\n",
              " 9.741525584287155\n",
              " 8.551739787702976\n",
              " 9.881899367462028\n",
              " 6.926693118258402\n",
              " 9.293861462481525"
            ]
          }
        }
      ],
      "source": [
        "Random.seed!(4750)\n",
        "rand(load, 5)"
      ],
      "id": "34"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is the operation at the heart of Monte Carlo: draw inputs, push\n",
        "each one through the model, and look at the distribution of outputs."
      ],
      "id": "d2e16488-94cb-4854-8ad7-0a81cd24f29f"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "(9.008633740674425, 1.4944133097855608)"
            ]
          }
        }
      ],
      "source": [
        "Random.seed!(4750)\n",
        "samples = rand(load, 10_000)\n",
        "mean(samples), std(samples)"
      ],
      "id": "36"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compare those to the true `mean(load)` and `std(load)` above. They are\n",
        "close but not exact, and that gap is the Monte Carlo error.\n",
        "\n",
        "### Estimating a probability\n",
        "\n",
        "To estimate the probability of an event, count the fraction of samples\n",
        "for which it happens:"
      ],
      "id": "07b75921-8fe5-4d6e-9624-1828a000bebb"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.0886"
            ]
          }
        }
      ],
      "source": [
        "mean(samples .> 11.0)        # estimated P(load > 11)"
      ],
      "id": "38"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.09121121972586788"
            ]
          }
        }
      ],
      "source": [
        "ccdf(load, 11.0)             # the exact answer, for comparison"
      ],
      "id": "40"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Taking the `mean` of a vector of `true`/`false` values gives the\n",
        "fraction that are `true`. That idiom appears throughout the course.\n",
        "\n",
        "### Seeds and reproducibility\n",
        "\n",
        "`rand` uses a **pseudorandom** number generator: the sequence looks\n",
        "random but is completely determined by a starting seed. Setting the seed\n",
        "makes your results reproducible, by you and by whoever is grading them."
      ],
      "id": "58e18004-6bd1-47e8-92f1-39e847183e61"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "true"
            ]
          }
        }
      ],
      "source": [
        "Random.seed!(1)\n",
        "a = rand(load, 3)\n",
        "Random.seed!(1)\n",
        "b = rand(load, 3)\n",
        "a == b"
      ],
      "id": "42"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Set a seed in every assignment that uses randomness.** It is worth a\n",
        "point on the Monte Carlo rubric, and without it nobody — including you —\n",
        "can reproduce what you did.\n",
        "\n",
        "## Conditional Probability\n",
        "\n",
        "The probability of $A$ given that $B$ happened:\n",
        "\n",
        "$$\\mathbb{P}(A \\mid B) = \\frac{\\mathbb{P}(A \\cap B)}{\\mathbb{P}(B)}$$\n",
        "\n",
        "Two events are **independent** when conditioning changes nothing, so\n",
        "$\\mathbb{P}(A \\mid B) = \\mathbb{P}(A)$ and\n",
        "$\\mathbb{P}(A \\cap B) = \\mathbb{P}(A)\\mathbb{P}(B)$.\n",
        "\n",
        "This matters in this course mostly through what it lets you multiply. If\n",
        "annual maxima are independent from year to year, the probability of no\n",
        "exceedance in $n$ years is $(1-p)^n$, so the probability of **at least\n",
        "one** is\n",
        "\n",
        "$$1 - (1 - p)^n.$$"
      ],
      "id": "2afd96c2-f2d5-4e2c-a9c3-0140fcdf9c64"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.2602996266117198"
            ]
          }
        }
      ],
      "source": [
        "p = 0.01                     # annual exceedance probability, a 100-year event\n",
        "n = 30                       # design life, years\n",
        "1 - (1 - p)^n"
      ],
      "id": "44"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A “100-year event” has better than a one-in-four chance of showing up\n",
        "during a 30-year design life. Independence is doing real work in that\n",
        "calculation — if exceedances cluster, the answer changes.\n",
        "\n",
        "## Getting Help\n",
        "\n",
        "- `?Normal` at the `julia>` prompt gives the constructor and its\n",
        "  parameters.\n",
        "- The [`Distributions.jl`\n",
        "  documentation](https://juliastats.org/Distributions.jl/stable/) lists\n",
        "  every distribution and the functions each supports.\n",
        "- For plotting any of this, see [Julia Plotting](julia-plots.qmd)."
      ],
      "id": "48f36aaa-2786-4ca8-9d7d-2757c5ee5bc1"
    }
  ],
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernel_info": {
      "name": "julia"
    },
    "kernelspec": {
      "name": "julia",
      "display_name": "Julia",
      "language": "julia"
    },
    "language_info": {
      "name": "julia",
      "codemirror_mode": "julia",
      "version": "1.11.5"
    }
  }
}