prob-surrogate-models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prob Surrogate Models\n",
    "\n",
    "This notebook was automatically generated from the Algorithms for Optimization source code. Each cell generates a figure from the original text. While this code is not optimized for use in lectures, we provide it here to be adapted for such projects. We hope you find it useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"support_code.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"gp.jl\")\n",
    "\n",
    "using LinearAlgebra\n",
    "\n",
    "p = let\n",
    "\n",
    "\tdomx = (-4,4)\n",
    "\tdomy = (-4,4)\n",
    "\n",
    "\tg = GroupPlot(3,1,groupStyle=\"horizontal sep=0.5cm, xlabels at=edge bottom, xticklabels at=edge bottom, ylabels at=edge left, yticklabels at=edge left\", style=\"width=6.5cm, height=6.5cm, axis equal, view={0}{90}, xlabel=\\$x_1\\$, ylabel=\\$x_2\\$,\")\n",
    "\n",
    "\tNs = [MvNormal(Float64[0,0], Matrix(1.0I, 2, 2)),\n",
    "\tMvNormal(Float64[0,0], Float64[3 0; 0 0.5]),\n",
    "\tMvNormal(Float64[0,0], [1 0.9; 0.9 1])]\n",
    "\ttitles = [L\"\\mat \\Sigma = \\begin{bmatrix}1 & 0 \\\\ 0 & 1\\end{bmatrix}\",\n",
    "\tL\"\\mat \\Sigma = \\begin{bmatrix}3 & 0 \\\\ 0 & 0.5\\end{bmatrix}\",\n",
    "\tL\"\\mat \\Sigma = \\begin{bmatrix}1 & 0.9 \\\\ 0.9 & 1\\end{bmatrix}\"\n",
    "\t]\n",
    "\n",
    "\tfor (N,title) in zip(Ns, titles)\n",
    "\tax = Axis(Plots.Contour((x,y)->pdf(N, [x,y]), domx, domy, levels=[0.001,0.01,0.05,0.1,0.2,0.3], xbins=201, ybins=201),\n",
    "\txmin=-3, xmax=3, ymin=0, ymax=0.6, title = title, style=\"contour/labels=false\",\n",
    "\t)\n",
    "\tpush!(g, ax)\n",
    "\tend\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\tg = GroupPlot(3,1,groupStyle=\"horizontal sep=0.5cm\", style=\"cycle list name = pastelcolors, width=6.5cm\")\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tX = [[x] for x in range(0, stop=10, length=201)]\n",
    "\tfunction get_plot(m, l, title)\n",
    "\t\tGP = GaussianProcess(m=m, k=(x, x′) -> exp(-(norm(x - x′)^2/(2*l^2))))\n",
    "\n",
    "\t\tp = Plots.Plot[]\n",
    "\n",
    "\t\tfor i in 1:5\n",
    "\t\t    y = rand(GP, X)\n",
    "\t\t    push!(p, Plots.Linear([x[1] for x in X], y))\n",
    "\t\tend\n",
    "\t\tax = Axis(p, xlabel=L\"x\",\n",
    "\t\t             xmin=0, xmax=10, ymin=-3.5, ymax=3.5, title=title)\n",
    "\t\tif length(g.axes) == 0\n",
    "\t\t    ax.ylabel = L\"y\"\n",
    "\t\telse\n",
    "\t\t    ax.style=\"yticklabels={,}\"\n",
    "\t\tend\n",
    "\t\tax\n",
    "\tend\n",
    "\n",
    "\tpush!(g, get_plot(x -> 0.0, 0.5, L\"\\ell = \\sfrac{1}{2}\"))\n",
    "\tpush!(g, get_plot(x -> 0.0, 1, L\"\\ell = 1\"))\n",
    "\tpush!(g, get_plot(x -> 0.0, 2, L\"\\ell = 2\"))\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "import SpecialFunctions: gamma, besselk\n",
    "\n",
    "matern(ν,ℓ) = (x,x′) -> begin\n",
    "    r = norm(x-x′)\n",
    "    if isapprox(r, 0.0)\n",
    "        return 1.0\n",
    "    end\n",
    "    g = sqrt(2ν)*r/ℓ\n",
    "    (2^(1-ν))/gamma(ν) * g^ν * besselk(ν,g)\n",
    "end\n",
    "\n",
    "get_nn_kernel(Σ) = (x,x′) -> begin\n",
    "    x = vcat(1.0, x)\n",
    "    x′ = vcat(1.0, x′)\n",
    "    asin((2*(x'*Σ*x′)[1])/sqrt((1+(2*x'*Σ*x)[1])*(1+(2*x′'*Σ*x′)[1])))\n",
    "end\n",
    "\n",
    "p = begin\n",
    "\tkernels = Function[]\n",
    "\ttitles = LaTeXString[]\n",
    "\n",
    "\tpush!(titles, L\"Constant: $\\sigma_0^2$\"); push!(kernels, (x,x′) -> 1.0)\n",
    "\tpush!(titles, L\"Linear: $\\displaystyle \\sum_{d=1}^{n} \\sigma^2_d x_d  x^\\prime_d$\"); push!(kernels, (x,x′) -> x⋅x′)\n",
    "\tpush!(titles, L\"Polynomial: $(x^\\top x^\\prime + \\sigma_0^2)^p$\"); push!(kernels, (x,x′) -> (x⋅x′+1)^2)\n",
    "\tpush!(titles, L\"Exponential: $\\exp(-\\frac{r}{\\ell})$\"); push!(kernels, (x,x′) -> exp(-abs(x[1]-x′[1])))\n",
    "\tpush!(titles, L\"$\\gamma$-Exponential: $\\exp(-(\\frac{r}{\\ell})^\\gamma)$\"); push!(kernels, (x,x′) -> exp(-(norm(x-x′)^0.5)))\n",
    "\tpush!(titles, L\"Squared Exponential: $\\exp(-\\frac{r^2}{2\\ell^2})$\"); push!(kernels, (x,x′) -> exp(-norm(x-x′)^2/2))\n",
    "\tpush!(titles, L\"Mat\\'{e}rn: $\\frac{2^{1-\\nu}}{\\Gamma(\\nu)} \\left(\\sqrt{2\\nu}\\frac{r}{\\ell}\\right)^\\nu K_\\nu \\left(\\sqrt{2\\nu}\\frac{r}{\\ell}\\right)$\"); push!(kernels, matern(0.5,1.0))\n",
    "\tpush!(titles, L\"Rational Quadratic: $\\left(1 + \\frac{r^2}{2\\alpha\\ell^2}\\right)^{-\\alpha}$\"); push!(kernels, (x,x′) -> (1+norm(x-x′)^2)^-0.5)\n",
    "\tpush!(titles, L\"Neural Network: $\\sin^{-1}\\left(\\frac{2\\bar{\\vect{x}}^\\top \\Sigma\\bar{\\vect{x}}^\\prime}{\\sqrt{(1+2\\bar{\\vect{x}}^\\top\\mat \\Sigma\\bar{\\vect{x}})(1+2\\bar{\\vect{x}}^{\\prime\\top}\\mat\\Sigma\\bar{\\vect{x}}^\\prime)}}\\right)$\"); push!(kernels, get_nn_kernel(Matrix(1.0I, 2, 2)))\n",
    "\n",
    "\tg = GroupPlot(3,3,groupStyle=\"horizontal sep=0.75cm, vertical sep=1.25cm, xlabels at=edge bottom, xticklabels at=edge bottom, ylabels at=edge left, yticklabels at=edge left\",\n",
    "\t                    style=\"cycle list name = pastelcolors, height=5cm, width=6cm, xmin=-5, xmax=5, ymin=-3.5, ymax=3.5, xlabel=\\$x\\$, ylabel=\\$y\\$, every axis title/.style={font=\\\\footnotesize, at={(0.5,1.15)}, align=center}\")\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tX = [[x] for x in range(-5, stop=5, length=101)]\n",
    "\n",
    "\tfunction get_plot(title, k)\n",
    "\t    p = Plots.Plot[]\n",
    "\t    GP = GaussianProcess(k=k)\n",
    "\t    for i = 1:5\n",
    "\t        y = rand(GP, X)\n",
    "\t        push!(p, Plots.Linear([x[1] for x in X], y))\n",
    "\t    end\n",
    "\t    Axis(p, title=title)\n",
    "\tend\n",
    "\n",
    "\tfor (title, k) in zip(titles, kernels)\n",
    "\t    push!(g, get_plot(title, k))\n",
    "\tend\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = begin\n",
    "\tg = GroupPlot(3,1,groupStyle=\"horizontal sep=0.75cm, vertical sep=1.5cm, xlabels at=edge bottom, ylabels at=edge left\",\n",
    "\t                       style=\"width=6cm, height=6cm, xlabel=\\$x_1\\$, ylabel=\\$x_2\\$, xtick=\\\\empty, ytick=\\\\empty, contour/labels=false, axis equal, view={0}{90}\")\n",
    "\n",
    "\tx_arr = range(0, stop=10, length=71)\n",
    "\ty_arr = range(0, stop=10, length=71)\n",
    "\n",
    "\tfor l in [0.5, 1.0, 2.0]\n",
    "\n",
    "\t\tGP = GaussianProcess(k = (x, x′) -> exp(-(norm(x - x′)^2/(2*l^2))))\n",
    "\n",
    "\t    # create random \"true\" function\n",
    "\t    Random.seed!(0)\n",
    "\t    X = vec([[x, y] for x in x_arr, y in y_arr])\n",
    "\t    z = rand(GP, X)\n",
    "\t    n = length(x_arr)\n",
    "\t    A = Float64[z[(i-1)*n+j] for i in 1:length(x_arr), j in 1 : length(y_arr)]\n",
    "\n",
    "\t    push!(g, Axis(Plots.Contour(A, x_arr, y_arr)))\n",
    "\tend\n",
    "\n",
    "\tg.axes[1].title = L\"\\ell = \\sfrac{1}{2}\"\n",
    "\tg.axes[2].title = L\"\\ell = 1\"\n",
    "\tg.axes[3].title = L\"\\ell = 2\"\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\n",
    "\t# create random \"true\" function\n",
    "\tRandom.seed!(0)\n",
    "\tGP = GaussianProcess()\n",
    "\tx_arr = range(0, stop=10, length=101)\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = rand(GP, X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true objective function\")\n",
    "\n",
    "\t# create training set\n",
    "\tfor i in [5, 20, 35, 45]\n",
    "\t\tpush!(GP, X[i], y[i])\n",
    "\tend\n",
    "\tp_train = Plots.Scatter([x[1] for x in GP.X], GP.y, style=\"solid, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "\n",
    "\t# predict\n",
    "\tμₚ, νₚ = predict(GP, X)\n",
    "\tp_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted function mean\")\n",
    "\tσₚ = sqrt.(νₚ)\n",
    "\tupperConfidence = μₚ + 1.96*σₚ\n",
    "\tlowerConfidence = μₚ - 1.96*σₚ\n",
    "\tp_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence interval}\")\n",
    "\n",
    "\tAxis([p_pred_hi, p_pred_lo, p_pred, p_true, p_pred_μ, p_train],\n",
    "\t\txmin=0, xmax=8, style=\"axis on top, ytick=\\\\empty, xtick=\\\\empty\",\n",
    "\t\txlabel=L\"x\", ylabel=L\"y\", width=\"7.5cm\", legendPos=\"outer north east\",\n",
    "\t)\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        struct GaussianProcessWithGradient\n",
    "            kff::Function\n",
    "            k∇f::Function\n",
    "            k∇∇::Function\n",
    "            X::Vector{Vector{Float64}}\n",
    "            y::Vector{Float64}\n",
    "            ∇y::Vector{Vector{Float64}}\n",
    "        end\n",
    "        function GaussianProcessWithGradient(α::Float64, γ::Float64)\n",
    "\n",
    "            squaregammanorm = (x, γ) -> (x.*γ)⋅x\n",
    "\n",
    "            kff = (x,x′) -> α*exp(-0.5*squaregammanorm(x-x′, γ))\n",
    "            k∇f = (x,x′,j) -> -γ*α*(x′[j]-x[j])*exp(-0.5*squaregammanorm(x-x′, γ))\n",
    "            k∇∇ = (x,x′,i,j) -> γ*α*((i == j ? 1 : 0) - γ*(x[j]-x′[j])*(x[i]-x′[i]))*exp(-0.5*squaregammanorm(x-x′, γ))\n",
    "\n",
    "            return GaussianProcessWithGradient(\n",
    "                kff,\n",
    "                k∇f,\n",
    "                k∇∇,\n",
    "                Vector{Float64}[],\n",
    "                Float64[],\n",
    "                Vector{Float64}[],\n",
    "            )\n",
    "        end\n",
    "\n",
    "        function Base.push!(GP::GaussianProcessWithGradient, x::Vector{Float64}, y::Real, ∇y::Vector{Float64})\n",
    "            push!(GP.X, x)\n",
    "            push!(GP.y, y)\n",
    "            push!(GP.∇y, ∇y)\n",
    "            return GP\n",
    "        end\n",
    "        function Base.pop!(GP::GaussianProcessWithGradient)\n",
    "            pop!(GP.X)\n",
    "            pop!(GP.y)\n",
    "            pop!(GP.∇y)\n",
    "            return GP\n",
    "        end\n",
    "\n",
    "        K(X::Dset, X′::Dset, k::Function) = [k(x,x′) for x in X, x′ in X′]\n",
    "        K(X::Dset, X′::Dset, k::Function, i::Int) = [k(x,x′,i) for x in X, x′ in X′]\n",
    "        K(X::Dset, X′::Dset, k::Function, i::Int, j::Int) = [k(x,x′,i,j) for x in X, x′ in X′]\n",
    "\n",
    "        function K∇f(X::Dset, X′::Dset, k::Function)\n",
    "            n = length(X)\n",
    "            m = length(X′)\n",
    "            d = length(X[1])\n",
    "            retval = Array{Float64}(undef, n, m*d)\n",
    "            for (i,x) in enumerate(X)\n",
    "                for (j,x′) in enumerate(X′)\n",
    "                    for a in 1 : d\n",
    "                        index = d*(j-1) + a\n",
    "                        retval[i, index] = k(x,x′,a)\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "            return retval\n",
    "        end\n",
    "        function K∇∇(X::Dset, X′::Dset, k::Function)\n",
    "            n = length(X)\n",
    "            m = length(X′)\n",
    "            d = length(X[1])\n",
    "            retval = Array{Float64}(undef, n*d, m*d)\n",
    "            for (i,x) in enumerate(X)\n",
    "                for (j,x′) in enumerate(X′)\n",
    "                    for a in 1 : d\n",
    "                        for b in 1 : d\n",
    "                            retval[d*(i-1) + a, d*(j-1) + b] = k(x,x′,a,b)\n",
    "                        end\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "            return retval\n",
    "        end\n",
    "\n",
    "        function predict(GPN::GaussianProcessWithGradient, X_pred::Dset)\n",
    "            _Kff = K(GPN.X, GPN.X, GPN.kff)\n",
    "            _K∇f = K∇f(GPN.X, GPN.X, GPN.k∇f)\n",
    "            _K∇∇ = K∇∇(GPN.X, GPN.X, GPN.k∇∇)\n",
    "\n",
    "            _Kpf = K(X_pred, GPN.X, GPN.kff)\n",
    "            _K∇p = K∇f(GPN.X, X_pred, GPN.k∇f)\n",
    "            _Kpp = K(X_pred, X_pred, GPN.kff)\n",
    "\n",
    "            v = hcat(_Kpf, _K∇p')'\n",
    "            V = vcat(hcat(_Kff, _K∇f'), hcat(_K∇f, _K∇∇))\n",
    "\n",
    "            μₚ = v'/V * vcat(GPN.y, -vcat(GPN.∇y...))\n",
    "            νₚ = diag(_Kpp - v'/V * v) .+ eps()\n",
    "            return (μₚ, νₚ)\n",
    "        end\n",
    "\n",
    "        p = let\n",
    "            α = 1.0\n",
    "            γ = 1.5\n",
    "\n",
    "            domain = (0,8)\n",
    "            f = x -> 1 - exp(-(x[1]-4)^2)\n",
    "\t\t\t∇f = x -> [2*exp(-(x[1]-4)^2)*(x[1]-4)]\n",
    "\n",
    "            squaregammanorm = (x, γ) -> (x.*γ)⋅x\n",
    "\n",
    "            GP = GaussianProcess(k=(x,x′) -> α*exp(-0.5*squaregammanorm(x-x′,γ)))\n",
    "            GPN = GaussianProcessWithGradient(α, γ)\n",
    "            x_arr = collect(range(domain[1], stop=domain[2], length=151))\n",
    "            X = [[x] for x in collect(x_arr)]\n",
    "            for x in [[1.0], [2.0], [3.5], [4.5]]\n",
    "                push!(GP, x, f(x))\n",
    "                push!(GPN, x, f(x), ∇f(x))\n",
    "            end\n",
    "\n",
    "            p_train = Plots.Scatter([x[1] for x in GPN.X], GPN.y, style=\"draw opacity=0, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "            p_true = Plots.Linear(x_arr, f.(X), style=\"solid, black, mark=none\", legendentry=\"true function\")\n",
    "\n",
    "            # predict w/o gradient\n",
    "            μₚ, νₚ = predict(GP, X)\n",
    "            p_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted mean without gradient\")\n",
    "            σₚ = sqrt.(νₚ)\n",
    "            upperConfidence = μₚ + 1.96*σₚ\n",
    "            lowerConfidence = μₚ - 1.96*σₚ\n",
    "            p_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "            p_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "            p_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence interval without gradient}\")\n",
    "\n",
    "            # predict with gradient\n",
    "            ∇μₚ, ∇νₚ = predict(GPN, X)\n",
    "            p_pred_∇μ = Plots.Linear(x_arr, ∇μₚ, style=\"solid, thick, pastelRed, mark=none\", legendentry=\"predicted mean with gradient\")\n",
    "            ∇σₚ = sqrt.(∇νₚ)\n",
    "            upperConfidence∇ = ∇μₚ + 1.96*∇σₚ\n",
    "            lowerConfidence∇ = ∇μₚ - 1.96*∇σₚ\n",
    "            p_pred_∇hi = Plots.Linear(x_arr, upperConfidence∇, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "            p_pred_∇lo = Plots.Linear(x_arr, lowerConfidence∇, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "            p_pred∇ = Plots.Command(\"\\\\addplot[pastelRed!70] fill between[of=A and B];\\\\addlegendentry{confidence interval with gradient}\")\n",
    "\n",
    "            Axis([p_pred_hi, p_pred_lo, p_pred, p_pred_∇hi, p_pred_∇lo, p_pred∇, p_pred_μ, p_pred_∇μ, p_train, p_true],\n",
    "              style=\"enlarge x limits=0, axis on top, ytick=\\\\empty, xtick=\\\\empty\",\n",
    "              xlabel=L\"x\", ylabel=L\"y\", xmin=0, xmax=8, width=\"7.5cm\", legendPos=\"outer north east\"\n",
    "            )\n",
    "        end\n",
    "        plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\n",
    "\tGP = GaussianProcess(\n",
    "\t\tm = x -> 0.0,\n",
    "\t\tk = (x, x′) -> exp(-(norm(x - x′))^2), # squared exponential\n",
    "\t\tν = 0.05,\n",
    "\t)\n",
    "\n",
    "\t# create random \"true\" function\n",
    "\tRandom.seed!(0)\n",
    "\tx_arr = range(0, stop=10, length=101)\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = rand(GP, X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true objective function\")\n",
    "\n",
    "\t# create training set\n",
    "\tfor i in [5, 20, 35, 45, 45]\n",
    "\t\tpush!(GP, X[i], y[i] + randn().*sqrt(GP.ν))\n",
    "\tend\n",
    "\tp_train = Plots.Scatter([x[1] for x in GP.X], GP.y, style=\"draw opacity=0, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "\n",
    "\t# predict\n",
    "\tμₚ, νₚ = predict(GP, X)\n",
    "\tσₚ = sqrt.(νₚ)\n",
    "\tp_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted function mean\")\n",
    "\tupperConfidence = μₚ + 1.96*σₚ\n",
    "\tlowerConfidence = μₚ - 1.96*σₚ\n",
    "\tp_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence region}\")\n",
    "\n",
    "\tAxis([p_pred_hi, p_pred_lo, p_pred, p_true, p_pred_μ, p_train],\n",
    "\t\txmin=0, xmax=8, style=\"axis on top, ytick=\\\\empty, xtick=\\\\empty\",\n",
    "\t\txlabel=L\"x\", ylabel=L\"y\", width=\"7.5cm\", legendPos=\"outer north east\",\n",
    "\t)\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}