surrogate-models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "include(\"support_code.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "using LinearAlgebra\n",
    "\n",
    "function design_matrix(X)\n",
    "\tn, m = length(X[1]), length(X)\n",
    "\treturn [j==0 ? 1. : X[i][j] for i in 1:m, j in 0:n]\n",
    "end\n",
    "function linear_regression(X, y)\n",
    "\tθ = pinv(design_matrix(X))*y\n",
    "\treturn x -> θ⋅[1; x]\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\tG = GroupPlot(2,2,groupStyle=\"horizontal sep=01cm, vertical sep=1.5cm\", style=\"xlabel=\\$x\\$, ylabel=\\$y\\$, width=6cm, ytick=\\\\empty, xtick=\\\\empty, enlargelimits=0, ymin=-0.2, ymax=0.4, xmin=0, xmax=1\")\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tf = x -> 0.2x + 0.2randn()\n",
    "\n",
    "\tfunction get_plot(X, ys, title)\n",
    "\t\tfhat = linear_regression(X, ys)\n",
    "\t\tplots = Plots.Plot[]\n",
    "\t\tpush!(plots, Plots.Scatter(X, ys, style=\"mark=*, mark size=1, mark options={draw=black, fill=black}\"))\n",
    "\t\tfor (x, y)  in zip(X, ys)\n",
    "\t\t\tdist = abs(y - fhat(x))\n",
    "\t\t\tif dist > 0.01\n",
    "\t\t\t\tpush!(plots, Plots.Linear([x,x], [y, fhat(x)], style=\"solid, gray, mark=none\"))\n",
    "\t\t\tend\n",
    "\t\tend\n",
    "\t\tpush!(plots, Plots.Linear([0,1], [fhat(0), fhat(1)], style=\"solid, pastelBlue, mark=none\"))\n",
    "\t\tAxis(plots, title=title)\n",
    "\tend\n",
    "\n",
    "\tX = [0.2]\n",
    "\tys = f.(X)\n",
    "\tpush!(G, get_plot(X, ys, L\"m < n+1\"))\n",
    "\n",
    "\tpush!(X, 0.6)\n",
    "\tpush!(ys, f(X[end]))\n",
    "\tpush!(G, get_plot(X, ys, L\"m = n+1\"))\n",
    "\n",
    "\tpush!(G, get_plot([X[1], X[1]], [ys[1], f(X[1])], \"nonindependent points\"))\n",
    "\n",
    "\tappend!(X, [0.1, 0.35, 0.4, 0.4, 0.8])\n",
    "\tappend!(ys, f.(X[end-4:end]))\n",
    "\tpush!(G, get_plot(X, ys, L\"m > n+1\"))\n",
    "\n",
    "\tG\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "using LinearAlgebra\n",
    "\n",
    "polynomial_bases_1d(i, k) = [x->x[i]^p for p in 0:k]\n",
    "function polynomial_bases(n, k)\n",
    "\tbases = [polynomial_bases_1d(i, k) for i in 1 : n]\n",
    "\tterms = Function[]\n",
    "\tfor ks in Iterators.product([0:k for i in 1:n]...)\n",
    "\t\tif sum(ks) ≤ k\n",
    "\t\t\tpush!(terms,\n",
    "\t\t\t\tx->prod(b[j+1](x) for (j,b) in zip(ks,bases)))\n",
    "\t\tend\n",
    "\tend\n",
    "\treturn terms\n",
    "end\n",
    "function regression(X, y, bases)\n",
    "    B = [b(x) for x in X, b in bases]\n",
    "    θ = pinv(B)*y\n",
    "    return x -> sum(θ[i] * bases[i](x) for i in 1 : length(θ))\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\ta = 45\n",
    "\tb = 180+123\n",
    "\tG = GroupPlot(2,1, groupStyle=\"horizontal sep = 1cm\", style=\"width=6cm\")\n",
    "\n",
    "\tf = x -> x^2 + 0.1\n",
    "\tRandom.seed!(0)\n",
    "\tX = [0.1, 0.2, 0.35, 0.4, 0.6, 0.8]\n",
    "\ty = f.(X)\n",
    "\n",
    "\tfhat = regression(X, y, polynomial_bases(1,2))\n",
    "\n",
    "\tpush!(G, Axis([\n",
    "\t        Plots.Scatter(X, y, style=\"mark size=1, mark options={draw=black, fill=black}\")\n",
    "\t        Plots.Linear(x->fhat([x]), (0,1), style=\"solid, pastelBlue, mark=none\")\n",
    "\t        ], ylabel=L\"y\", xlabel=L\"x\", style=\"xticklabels=\\\\empty, yticklabels=\\\\empty\"))\n",
    "\n",
    "\tx_arr = collect(range(0, stop=1, length=101))\n",
    "\tx2_arr = x_arr.^2\n",
    "\ty_arr = [fhat(x) for x in x_arr]\n",
    "\n",
    "\tpush!(G, Axis([\n",
    "\t        Plots.Command(\"\\\\addplot3[patch,patch type=rectangle, faceted color=none, color=pastelBlue!40] coordinates {(0,0,0.1) (1,0,0.1) (1,1,1.1) (0,1,1.1)}\"),\n",
    "\t        Plots.Linear3(X, X.^2, y, mark=\"*\", style=\"mark size=1.1, only marks, mark options={draw=black, fill=black}\", legendentry=L\"f(x)\"),\n",
    "\t        Plots.Linear3(x_arr, x2_arr, y_arr, style=\"pastelBlue, solid, mark=none\", legendentry=L\"\\hat{f}(x)\"),\n",
    "\t    ], xlabel=L\"x\", ylabel=L\"x^2\", zlabel=L\"y\", width=\"6cm\", style=\"view={45}{20}, axis equal, xticklabels=\\\\empty, yticklabels=\\\\empty, zticklabels=\\\\empty\", legendPos=\"outer north east\"))\n",
    "\n",
    "\tG\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "function sinusoidal_bases_1d(j, k, a, b)\n",
    "\tT = b[j] - a[j]\n",
    "\tbases = Function[x->1.0]\n",
    "\tfor i in 1 : k\n",
    "\t\tpush!(bases, x->sin(2π*i*x[j]/T))\n",
    "\t\tpush!(bases, x->cos(2π*i*x[j]/T))\n",
    "\tend\n",
    "\treturn bases\n",
    "end\n",
    "function sinusoidal_bases(k, a, b)\n",
    "\tn = length(a)\n",
    "\tbases = [sinusoidal_bases_1d(i, k, a, b) for i in 1 : n]\n",
    "\tterms = Function[]\n",
    "\tfor ks in Iterators.product([0:2k for i in 1:n]...)\n",
    "\t\tpowers = [div(k+1,2) for k in ks]\n",
    "\t\tif sum(powers) ≤ k\n",
    "\t\t\tpush!(terms,\n",
    "\t\t\t\tx->prod(b[j+1](x) for (j,b) in zip(ks,bases)))\n",
    "\t\tend\n",
    "\tend\n",
    "\treturn terms\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\ta = [0.0]\n",
    "\tb = [10.0]\n",
    "\tf = x -> sin(x) + 0.2randn()\n",
    "\n",
    "\tG = GroupPlot(2,2,groupStyle=\"horizontal sep=01cm, vertical sep=1.5cm\", style=\"xlabel=\\$x\\$, ylabel=\\$y\\$, width=6cm, ytick=\\\\empty, xtick=\\\\empty, enlargelimits=0, ymin=-2, ymax=2, xmin=$(a[1]), xmax=$(b[1])\")\n",
    "\n",
    "\tfunction get_plot(X, ys, title)\n",
    "\t\tfhat = regression(X, ys, sinusoidal_bases(1, a, b))\n",
    "\t\tplots = Plots.Plot[]\n",
    "\t\tpush!(plots, Plots.Scatter(X, ys, style=\"mark=*, mark size=1, mark options={draw=black, fill=black}\"))\n",
    "\t\tfor (x, y)  in zip(X, ys)\n",
    "\t\t\tdist = abs(y - fhat(x))\n",
    "\t\t\tif dist > 0.01\n",
    "\t\t\t\tpush!(plots, Plots.Linear([x,x], [y, fhat([x])], style=\"solid, gray, mark=none\"))\n",
    "\t\t\tend\n",
    "\t\tend\n",
    "\t\tpush!(plots, Plots.Linear(x->fhat(x), (a[1],b[1]), style=\"solid, pastelBlue, mark=none\"))\n",
    "\t\tAxis(plots, title=title)\n",
    "\tend\n",
    "\n",
    "\tX = [2.0]\n",
    "\tys = f.(X)\n",
    "\tpush!(G, get_plot(X, ys, L\"m < n+1\"))\n",
    "\n",
    "\tpush!(X, 5.6)\n",
    "\tpush!(ys, f(X[end]))\n",
    "\tpush!(G, get_plot(X, ys, L\"m = n+1\"))\n",
    "\n",
    "\tpush!(G, get_plot([X[1], X[1]], [ys[1], f(X[1])], \"nonindependent points\"))\n",
    "\n",
    "\tappend!(X, [1.1, 9.35, 4, 4.4, 8.8])\n",
    "\tappend!(ys, f.(X[end-4:end]))\n",
    "\tpush!(G, get_plot(X, ys, L\"m > n+1\"))\n",
    "\n",
    "\tG\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\txdom = (0,2)\n",
    "\tydom = (-2,2)\n",
    "\n",
    "\tg = GroupPlot(3,2, groupStyle=\"horizontal sep=1.75cm, vertical sep=1cm, xlabels at=edge bottom, ylabels at=edge left, xticklabels at=edge bottom, yticklabels at=edge bottom\", style=\"xlabel=\\$r\\$, ylabel=\\$\\\\psi\\$, title style={text height=2ex}\")\n",
    "\n",
    "\tfunction add_plot!(f, title)\n",
    "\t    push!(g,\n",
    "\t        Axis(\n",
    "\t            Plots.Linear(r -> f(abs(r)), xdom, style=\"thick, solid, black, mark=none\"),\n",
    "\t            width=\"5cm\", height=\"5cm\", title=title,\n",
    "\t        ),\n",
    "\t    )\n",
    "\tend\n",
    "\n",
    "\tadd_plot!(r -> r, \"linear: \\$r\\$\")\n",
    "\tadd_plot!(r -> r^3, \"cubic: \\$r^3\\$\")\n",
    "\tadd_plot!(r -> r^2 * log(r), \"thin plate spline: \\$r^2 \\\\log r\\$\")\n",
    "\tadd_plot!(r -> exp(-r^2), \"Gaussian: \\$e^{-r^2/2\\\\sigma^2} \\$\")\n",
    "\tadd_plot!(r -> (r^2 + 1)^0.5, \"multiquadratic: \\$(r^2 + \\\\sigma^2)^\\\\frac{1}{2}\\$\")\n",
    "\tadd_plot!(r -> (r^2 + 1)^-0.5, \"inverse multiquadratic: \\$(r^2 + \\\\sigma^2)^{-\\\\frac{1}{2}}\\$\")\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "radial_bases(ψ, C, p=2) = [x->ψ(norm(x - c, p)) for c in C]\n",
    "\n",
    "p = let\n",
    "\n",
    "\tf = x -> x * sin(5x)\n",
    "\n",
    "\tRandom.seed!(6)\n",
    "\tX = [rand() for i in 1 : 4]\n",
    "\ty = (x->f(x[1])).(X)\n",
    "\n",
    "\tAxis([\n",
    "\t        Plots.Linear(f, (0,1), style=\"thick, solid, black, mark=none\", legendentry=L\"x \\sin(5x)\"),\n",
    "\t        Plots.Linear(regression(X, y, radial_bases(r -> exp(-2*r^2), X)), (0,1), style=\"thick, pastelBlue!40\", legendentry=L\"\\psi = \\exp(-2r^2)\"),\n",
    "\t        Plots.Linear(regression(X, y, radial_bases(r -> exp(-5*r^2), X)), (0,1), style=\"thick, pastelBlue!70\", legendentry=L\"\\psi = \\exp(-5r^2)\"),\n",
    "\t        Plots.Linear(regression(X, y, radial_bases(r -> exp(-10*r^2), X)), (0,1), style=\"thick, pastelBlue\", legendentry=L\"\\psi = \\exp(-10r^2)\"),\n",
    "\t        Plots.Scatter(X, y, style=\"black, mark=*, mark size=1, mark options={draw=none, fill=black}\"),\n",
    "\t        ], xlabel=L\"x\", ylabel=L\"y\", style=\"enlarge x limits=0\", width=\"8cm\", legendPos=\"outer north east\"\n",
    "\t    )\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "function regression(X, y, bases, λ)\n",
    "    B = [b(x) for x in X, b in bases]\n",
    "    θ = (B'B + λ^2*I)\\B'y\n",
    "    return x -> sum(θ[i] * bases[i](x) for i in 1 : length(θ))\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\tf = x -> x * sin(5x)\n",
    "\n",
    "\tRandom.seed!(6)\n",
    "\tX = rand(10)\n",
    "\ty = f.(X) + randn(length(X))/10\n",
    "\n",
    "\tAxis([\n",
    "\t        Plots.Linear(f, (0,1), style=\"thick, solid, black, mark=none\", legendentry=L\"x \\sin (5x)\"),\n",
    "\t        Plots.Linear(regression(X, y, radial_bases(r -> exp(-5*r^2), X)), (0,1), style=\"thick, pastelBlue\", legendentry=L\"\\lambda = 0\"),\n",
    "\t        Plots.Linear(regression(X, y, radial_bases(r -> exp(-5*r^2), X), 0.1), (0,1), style=\"thick, pastelPurple\", legendentry=L\"\\lambda = 0.1\"),\n",
    "\t        Plots.Linear(regression(X, y, radial_bases(r -> exp(-5*r^2), X), 0.5), (0,1), style=\"thick, pastelRed\", legendentry=L\"\\lambda = 0.5\"),\n",
    "\t        Plots.Scatter(X, y, style=\"black, mark=*, mark size=1, mark options={draw=none, fill=black}\"),\n",
    "\t        ], xlabel=L\"x\", ylabel=L\"y\", style=\"enlarge x limits=0\", width=\"8cm\", legendPos=\"outer north east\", ymin=-2, ymax=2,\n",
    "\t    )\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\tvariance = x->log(1+exp(5(x-1))) + 0.1\n",
    "\tbias = x->variance(1-x)\n",
    "\ttoterr = x->bias(x) + variance(x)\n",
    "\n",
    "\tAxis([\n",
    "\t        Plots.Linear(variance, (0,1), style=\"solid, pastelBlue, mark=none\", legendentry=\"variance\"),\n",
    "\t        Plots.Linear(bias, (0,1), style=\"solid, pastelRed, mark=none\", legendentry=L\"\\text{bias}^2\"),\n",
    "\t        Plots.Linear(toterr, (0,1), style=\"solid, black, mark=none\", legendentry=\"total error\"),\n",
    "\t        Plots.Linear([0.5,0.5],[0.1,0.9], style=\"solid, gray, mark=none\", legendentry=\"optimal tradeoff\"),\n",
    "\t    ],\n",
    "\t    width=\"9cm\", xlabel=\"complexity\", ylabel=\"error\",\n",
    "\t    ymin=0, ymax=1, xmin=0, xmax=1,\n",
    "\t    style=\"axis lines=middle, xtick=\\\\empty, ytick=\\\\empty, legend cell align=left, legend style={draw=none, fill=none, at={(0.5,-0.25)}, anchor=north, legend columns=1}, x label style={at={(axis description cs:0.75,-0.05)},anchor=north},\",\n",
    "\t)\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "using Distributions\n",
    "using LinearAlgebra\n",
    "\n",
    "function regression(X, y, bases)\n",
    "    B = [b(x) for x in X, b in bases]\n",
    "    θ = pinv(B)*y\n",
    "    return x -> sum(θ[i] * bases[i](x) for i in 1 : length(θ))\n",
    "end\n",
    "\n",
    "f = x->2/(1+exp(-x))\n",
    "r = 4\n",
    "\n",
    "str = \"\"\"\n",
    "\\\\begin{tikzpicture}[x=0.5cm, y=0.5cm]\n",
    "    \\\\coordinate (A) at (2,6);\")\n",
    "    \\\\coordinate (B) at (\\$(A) + (4.5,0)\\$);\n",
    "    \\\\coordinate (C) at (\\$(A) + (0,-4.5)\\$);\n",
    "    \\\\coordinate (D) at (C -| B);\n",
    "\n",
    "    \\\\draw (A) circle (1.75);\n",
    "    \\\\draw[gray] (\\$(A) + (0:2)\\$) -- ++(180:4);\n",
    "    \\\\draw[gray] (\\$(A) + (90:2)\\$) -- ++(270:4);\n",
    "    \\\\draw (B) circle (1.75);\n",
    "    \\\\draw[gray] (\\$(B) + (0:2)\\$) -- ++(180:4);\n",
    "    \\\\draw[gray] (\\$(B) + (90:2)\\$) -- ++(270:4);\n",
    "    \\\\draw (C) circle (1.75);\n",
    "    \\\\draw[gray] (\\$(C) + (0:2)\\$) -- ++(180:4);\n",
    "    \\\\draw[gray] (\\$(C) + (90:2)\\$) -- ++(270:4);\n",
    "    \\\\draw (D) circle (1.75);\n",
    "    \\\\draw[gray] (\\$(D) + (0:2)\\$) -- ++(180:4);\n",
    "    \\\\draw[gray] (\\$(D) + (90:2)\\$) -- ++(270:4);\n",
    "\n",
    "    \\\\node[rotate=90, anchor=south] at (\\$(A) + (-2.25,0)\\$) {low bias\\\\vphantom{g}};\n",
    "    \\\\node[rotate=90, anchor=south] at (\\$(C) + (-2.25,0)\\$) {high bias};\n",
    "    \\\\node[anchor=south] at (\\$(A) + (0,2.25)\\$) {low variance\\\\vphantom{g}};\n",
    "    \\\\node[anchor=south] at (\\$(B) + (0,2.25)\\$) {high variance};\n",
    "    \"\"\"\n",
    "\n",
    "function print_coordinate(str, p, f1, f2, x1, x2)\n",
    "    x = f1(x1) - f(x1)\n",
    "    y = f2(x2) - f(x2)\n",
    "    if hypot(x,y) < 3\n",
    "        str *= \"\\n\\\\fill[red, opacity=0.5] (\\$($p) + ($x, $y)\\$) circle (0.05);\"\n",
    "    end\n",
    "    return str\n",
    "end\n",
    "\n",
    "Random.seed!(0)\n",
    "for i in 1 : 100\n",
    "\n",
    "    X1 = collect(range(-r, stop=r, length=10))\n",
    "    X2 = collect(range(-r, stop=r, length=10))\n",
    "    y1 = (x->f(x) + 0.25randn()).(X1)\n",
    "    y2 = (x->f(x) + 0.25randn()).(X2)\n",
    "\n",
    "    f11 = regression(X1, y1, [x->1.0, x->x, x->x^3])\n",
    "    f12 = regression(X2, y2, [x->1.0, x->x, x->x^3])\n",
    "    f21 = regression(X1, y1, [x->1.0, x->x, x->x^2, x->x^3, x->x^4, x->x^5, x->x^6, x->x^7, x->x^8, x->x^9])\n",
    "    f22 = regression(X2, y2, [x->1.0, x->x, x->x^2, x->x^3, x->x^4, x->x^5, x->x^6, x->x^7, x->x^8, x->x^9])\n",
    "    f31 = regression(X1, y1, [x->x])\n",
    "    f32 = regression(X2, y2, [x->x])\n",
    "    f41 = regression(X1, y1, [x->sin(x/4), x->sin(x/8), x->sin(x/16), x->sin(x/32), x->sin(x/64)])\n",
    "    f42 = regression(X2, y2, [x->sin(x/4), x->sin(x/8), x->sin(x/16), x->sin(x/32), x->sin(x/64)])\n",
    "\n",
    "    x1 = rand(Uniform(-4,4))\n",
    "    x2 = rand(Uniform(-4,4))\n",
    "\n",
    "    str = print_coordinate(str, \"A\", f11, f12, x1, x2)\n",
    "    str = print_coordinate(str, \"B\", f21, f22, x1, x2)\n",
    "    str = print_coordinate(str, \"C\", f31, f32, x1, x2)\n",
    "    str = print_coordinate(str, \"D\", f41, f42, x1, x2)\n",
    "end\n",
    "\n",
    "str *= \"\\\\end{tikzpicture}\"\n",
    "\n",
    "using TikzPictures\n",
    "pic = TikzPicture(str, preamble=\"\\\\usetikzlibrary{calc}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\tg = GroupPlot(2,1)\n",
    "\n",
    "\tf = x -> x + x^2 + randn()/2\n",
    "\tRandom.seed!(0)\n",
    "\tX = [0.1, 0.2, 0.35, 0.4, 0.6, 0.8]\n",
    "\ty = f.(X)\n",
    "\n",
    "\tg = GroupPlot(3,2)\n",
    "\tfor k in 0 : 5\n",
    "\t    fhat = regression(X, y, polynomial_bases(1, k))\n",
    "\t    err = @sprintf(\"%.4f\", norm(fhat.(X) - y,2)^2)\n",
    "\t    push!(g,\n",
    "\t    Axis([\n",
    "\t        Plots.Scatter(X, y, style=\"mark size=1, mark options={draw=pastelBlue, fill=pastelBlue}\"),\n",
    "\t        Plots.Linear(x->fhat([x]), (0,1), style=\"solid, black, mark=none\"),\n",
    "\t        ], title=\"{\\$k=$k, \\\\text{err} = $(err)\\$}\", style=\"xticklabels=\\\\empty, yticklabels=\\\\empty\", legendPos=\"outer north east\", width=\"4.5cm\", ymin=0, ymax=2))\n",
    "\tend\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\tusing Polynomials\n",
    "\tusing Distributions\n",
    "\tusing QuadGK\n",
    "\tusing Random\n",
    "\tusing LinearAlgebra\n",
    "\n",
    "\tp = let\n",
    "\n",
    "\t\tftrue = x -> x + x^2\n",
    "\t\tf = x -> ftrue(x) + randn()/2\n",
    "\t\tRandom.seed!(0)\n",
    "\t\tX = rand(Uniform(0,1), 6)\n",
    "\t\ty = f.(X)\n",
    "\n",
    "\t\tk_arr = collect(0:5)\n",
    "\t\te_train = Array{Float64}(length(k_arr))\n",
    "\t\te_gen = Array{Float64}(length(k_arr))\n",
    "\n",
    "\t\tfor (i,k) in enumerate(k_arr)\n",
    "\t\t\tbases = polynomial_bases(1, k)\n",
    "\t\t\tB = [b(x) for x in X, b in bases]\n",
    "    \t\tθ = pinv(B)*y\n",
    "\t\t    fhat = x -> sum(θ[i] * bases[i](x) for i in 1 : length(θ))\n",
    "\t\t    e_train[i] = norm(fhat.(X) - y,2)^2 / length(y)\n",
    "\t\t    e_gen[i] = QuadGK.quadgk(x->(ftrue(x) - fhat([x]))^2, 0.0, 1.0)[1]\n",
    "\t\tend\n",
    "\n",
    "\t\tAxis([\n",
    "\t\t    Plots.Linear(k_arr[1:5], e_train[1:5], style=\"solid, pastelBlue, mark=*, mark size=1, mark options={draw=pastelBlue, fill=pastelBlue}\", legendentry=\"training error\")\n",
    "\t\t    Plots.Linear(k_arr, e_gen, style=\"solid, pastelRed, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\", legendentry=\"generalization error\")\n",
    "\t\t    ], width=\"6cm\", xlabel=\"polynomial degree \\$k\\$\", ylabel=\"error\", ymode=\"log\", legendPos=\"outer north east\")\n",
    "\tend\n",
    "\n",
    "\tplot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Statistics\n",
    "import QuadGK: quadgk\n",
    "p = let\n",
    "\n",
    "\tfunction regression(X, y, bases)\n",
    "\t    B = [b(x) for x in X, b in bases]\n",
    "\t    θ = pinv(B)*y\n",
    "\t    return x -> sum(θ[i] * bases[i](x) for i in 1 : length(θ))\n",
    "\tend\n",
    "\n",
    "\tpolynomial_bases_1d(i, k) = [x->x[i]^p for p in 0:k]\n",
    "\tfunction polynomial_bases(n, k)\n",
    "\t\tbases = [polynomial_bases_1d(i, k) for i in 1 : n]\n",
    "\t\tterms = Function[]\n",
    "\t\tfor ks in Iterators.product([0:k for i in 1:n]...)\n",
    "\t\t\tif sum(ks) ≤ k\n",
    "\t\t\t\tpush!(terms,\n",
    "\t\t\t\t\tx->prod(b[j+1](x) for (j,b) in zip(ks,bases)))\n",
    "\t\t\tend\n",
    "\t\tend\n",
    "\t\treturn terms\n",
    "\tend\n",
    "\n",
    "\txdom = (-5,5)\n",
    "\tf = x -> x/10 + sin(x)/4 + exp(-x^2)\n",
    "\n",
    "\tx_arr = collect(range(-4, stop=4, length=9))\n",
    "\ty_arr = f.(x_arr)\n",
    "\n",
    "\terr_gen(f, fhat, x_arr) = quadgk(x->(f(x) - fhat(x))^2 * 1/(xdom[2] - xdom[1]), xdom...)[1]\n",
    "\terr_train(f, fhat, x_arr) = mean((f.(x_arr) - fhat.(x_arr)).^2)\n",
    "\n",
    "\tG = GroupPlot(3,3,groupStyle=\"horizontal sep=0.25cm, vertical sep=0.25cm, xlabels at=edge bottom, xticklabels at=edge bottom, ylabels at=edge left, yticklabels at=edge left\",\n",
    "\t              style=\"width=4.75cm, ymin=-1, ymax=1.25, enlarge x limits=0, xlabel=\\$x\\$, ylabel=\\$y\\$, axis on top=true\")\n",
    "\n",
    "\tfor k in 0 : 8\n",
    "\n",
    "\t    fhat = regression(x_arr, y_arr, polynomial_bases(1, k))\n",
    "\n",
    "\t    plots = Plots.Plot[]\n",
    "\t    push!(plots, Plots.Linear(f, xdom, style=\"solid, black, mark=none\"))\n",
    "\t    push!(plots, Plots.Scatter(x_arr, y_arr, style=\"only marks, solid, mark size=1, mark=*, mark options={fill=black, draw=black}\"))\n",
    "\t    push!(plots, Plots.Linear(fhat, xdom, style=\"solid, thick, pastelBlue, mark=none\"))\n",
    "\t    push!(plots, Plots.Node(\"\\\\footnotesize \\$\\\\epsilon_{train} = $(round(err_train(f, fhat, x_arr), digits=3))\\$\", 0.4, 0.3, axis=\"axis description cs\", style=\"right\"))\n",
    "\t    push!(plots, Plots.Node(\"\\\\footnotesize \\$\\\\epsilon_{gen} = $(round(err_gen(f, fhat, x_arr), digits=3))\\$\", 0.4, 0.15, axis=\"axis description cs\", style=\"right\"))\n",
    "\t    push!(plots, Plots.Node(\"\\\\footnotesize \\$k = $k\\$\", 0.05, 0.87, axis=\"axis description cs\", style=\"right\"))\n",
    "\n",
    "\t    push!(G, Axis(plots))\n",
    "\tend\n",
    "\tG\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\n",
    "\tf = x -> sin(x) + randn()*0.1\n",
    "\tRandom.seed!(0)\n",
    "\tX = Vector{Float64}[]\n",
    "\tpush!(X, [0.1])\n",
    "\tpush!(X, [0.2])\n",
    "\tpush!(X, [0.35])\n",
    "\tpush!(X, [0.4])\n",
    "\tpush!(X, [0.6])\n",
    "\tpush!(X, [0.8])\n",
    "\ty = (x->f(x[1])).(X)\n",
    "\n",
    "\tfhat = regression(X, y, polynomial_bases(1, 2))\n",
    "\n",
    "\tX2 = Vector{Float64}[]\n",
    "\tpush!(X2, [2.6])\n",
    "\tpush!(X2, [2.8])\n",
    "\ty2 = (x->f(x[1])).(X2)\n",
    "\n",
    "\tAxis([\n",
    "\t        Plots.Linear(sin, (0,3), style=\"solid, black, mark=none\", legendentry=L\"f\")\n",
    "\t        Plots.Linear(x->fhat([x]), (0,3), style=\"solid, pastelBlue, mark=none\", legendentry=L\"\\hat{f}\")\n",
    "\t        Plots.Scatter([x[1] for x in X], y, style=\"only marks, mark=*, mark size=1, mark options={draw=pastelBlue, fill=pastelBlue}\", legendentry=\"training samples\")\n",
    "\t        Plots.Scatter([x[1] for x in X2], y2, style=\"only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\", legendentry=\"holdout samples\")\n",
    "\t        ], ylabel=L\"y\", xlabel=L\"x\", style=\"xticklabels=\\\\empty, yticklabels=\\\\empty, legend cell align=left, legend style={draw=none, at={(0.5,-0.25)}, anchor=north, legend columns=1}\", width=\"9cm\")\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\tRandom.seed!(0)\n",
    "\txs = range(0, stop=3π, length=10)\n",
    "\ty = sin.(xs)\n",
    "\tX = [[x] for x in xs]\n",
    "\n",
    "\tfit = (X,y) -> regression(X, y, polynomial_bases(1,4))\n",
    "\tmetric = (f, X, y) -> begin\n",
    "\t    m = length(X)\n",
    "\t    return sum((f(X[i]) - y[i])^2 for i in m)/m\n",
    "\tend\n",
    "\n",
    "\tsets = [TrainTest(collect(1:9),collect(10:10)),\n",
    "\t         TrainTest(collect(1:8),collect(9:10)),\n",
    "\t         TrainTest(collect(1:5),collect(6:10)),\n",
    "\t        ]\n",
    "\n",
    "\tmodels = [fit(X[tt.train], y[tt.train]) for tt in sets]\n",
    "\tmetrics = [metric(model, X[tt.test], y[tt.test]) for (tt,model) in zip(sets,models)]\n",
    "\n",
    "\tplots = Plots.Plot[]\n",
    "\tpush!(plots, Plots.Scatter(xs, y, style=\"black, mark=*, mark size=1, mark options={draw=black, fill=black}\",  legendentry=\"dataset\"))\n",
    "\tpush!(plots, Plots.Linear(x->models[1]([x]), (-π/4,3π+π/4), style=\"solid, mark=none, pastelPurple\",  legendentry=\"10-fold\"))\n",
    "\tpush!(plots, Plots.Linear(x->models[2]([x]), (-π/4,3π+π/4), style=\"solid, mark=none, pastelBlue\",   legendentry=\" 5-fold\"))\n",
    "\tpush!(plots, Plots.Linear(x->models[3]([x]), (-π/4,3π+π/4), style=\"solid, mark=none, pastelSeaGreen\", legendentry=\" 2-fold\"))\n",
    "\tAxis(plots, xlabel=L\"x\", ylabel=L\"y\", width=\"10cm\", height=\"4cm\", ymin=-5, ymax=10,\n",
    "\tstyle=\"legend cell align=left, legend style={draw=none, at={(0.5,-0.15)}, anchor=north, legend columns=4}\"\n",
    "\t)\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "struct TrainTest\n",
    "\ttrain\n",
    "\ttest\n",
    "end\n",
    "function k_fold_cross_validation_sets(n, k)\n",
    "    @assert k ≤ n\n",
    "    perm = randperm(n)\n",
    "    sets = TrainTest[]\n",
    "    for i = 1:k\n",
    "        validate = perm[i:k:n];\n",
    "        train = perm[setdiff(1:n, i:k:n)]\n",
    "        push!(sets, TrainTest(train, validate))\n",
    "    end\n",
    "    return sets\n",
    "end\n",
    "function train_and_validate(X, y, tt, fit, metric)\n",
    "    model = fit(X[tt.train], y[tt.train])\n",
    "    return metric(model, X[tt.test], y[tt.test])\n",
    "end\n",
    "function cross_validation_estimate(X, y, sets, fit, metric)\n",
    "\taccuracies = [train_and_validate(X, y, tt, fit, metric) for tt in sets]\n",
    "    return (mean(accuracies), std(accuracies))\n",
    "end\n",
    "\n",
    "p = let\n",
    "\tf = x->sin(2x)*cos(10x)\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tX = rand(10)\n",
    "\ty = f.(X) + randn(length(X))/10\n",
    "\n",
    "\tsets = k_fold_cross_validation_sets(length(X), 3)\n",
    "\tmetric = (f, X, y)->begin\n",
    "\t\tm = length(X)\n",
    "\t\treturn sum((f(X[i]) - y[i])^2 for i in m)/m\n",
    "\tend\n",
    "\n",
    "\tλ_arr = collect(10 .^ range(-4, stop=2, length=101))\n",
    "\te_arr = Array{Float64}(undef, length(λ_arr))\n",
    "\tfor (i,λ) in enumerate(λ_arr)\n",
    "\t    fit = (X,y)->regression(X, y, radial_bases(r->exp(-5*r^2), X), λ)\n",
    "\t    e_arr[i] = cross_validation_estimate(X, y, sets, fit, metric)[1]\n",
    "\tend\n",
    "\n",
    "\tAxis(Plots.Linear(λ_arr, e_arr, style=\"solid, black, mark=none\"), xmode=\"log\", xlabel=L\"\\lambda\", ylabel=\"mean cross validated MSE\", width=\"9cm\")\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "function bootstrap_sets(m, b)\n",
    "\treturn [TrainTest(rand(1:m, m), 1:m) for i in 1 : b]\n",
    "end\n",
    "function bootstrap_estimate(X, y, sets, fit, metric)\n",
    "\tmean(train_and_validate(X, y, tt, fit, metric) for tt in sets)\n",
    "end\n",
    "function leave_one_out_bootstrap_estimate(X, y, sets, fit, metric)\n",
    "\tm, b = length(X), length(sets)\n",
    "\tε = 0.0\n",
    "\tmodels = [fit(X[tt.train], y[tt.train]) for tt in sets]\n",
    "\tfor j in 1 : m\n",
    "\t\tc = 0\n",
    "\t\tδ = 0.0\n",
    "\t\tfor i in 1 : b\n",
    "\t\t\tif j ∉ sets[i].train\n",
    "\t\t\t\tc += 1\n",
    "\t\t\t\tδ += metric(models[i], [X[j]], [y[j]])\n",
    "\t\t\tend\n",
    "\t\tend\n",
    "\t\tε += δ/c\n",
    "\tend\n",
    "\treturn ε/m\n",
    "end\n",
    "function bootstrap_632_estimate(X, y, sets, fit, metric)\n",
    "\tmodels = [fit(X[tt.train], y[tt.train]) for tt in sets]\n",
    "\tϵ_loob = leave_one_out_bootstrap_estimate(X,y,sets,fit,metric)\n",
    "\tϵ_boot = bootstrap_estimate(X,y,sets,fit,metric)\n",
    "    return 0.632ϵ_loob + 0.368ϵ_boot\n",
    "end\n",
    "function holdout_partition(m, h=div(m,2))\n",
    "    p = randperm(m)\n",
    "    train = p[(h+1):m]\n",
    "    holdout = p[1:h]\n",
    "    return TrainTest(train, holdout)\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\tfunction metric(f, X, y)\n",
    "\t    m = length(X)\n",
    "\t    return sqrt(sum((f(X[i]) - y[i])^2 for i in m)/m)\n",
    "\tend\n",
    "\n",
    "\tRandom.seed!(1)\n",
    "\txs = range(-3, stop=3, length=10)\n",
    "\tf = x -> x^2 + randn()/2\n",
    "\ty = f.(xs)\n",
    "\tX = [[x] for x in xs]\n",
    "\n",
    "\tfit = (X,y) -> regression(X, y, polynomial_bases(1,1))\n",
    "\n",
    "\tmins = Float64[]\n",
    "\tmaxs = Float64[]\n",
    "\tmedians = Float64[]\n",
    "\tquartile1s = Float64[]\n",
    "\tquartile3s = Float64[]\n",
    "\n",
    "\tbsets() = bootstrap_sets(length(X), 50)\n",
    "\tcvsets() = k_fold_cross_validation_sets(length(X), 5)\n",
    "\n",
    "\test_methods = [\n",
    "\t    (X,y)->bootstrap_632_estimate(X, y, bsets(), fit, metric),\n",
    "\t    (X,y)->leave_one_out_bootstrap_estimate(X, y, bsets(), fit, metric),\n",
    "\t    (X,y)->bootstrap_estimate(X, y, bsets(), fit, metric),\n",
    "\t    (X,y)->cross_validation_estimate(X, y, cvsets(), fit, metric)[1],\n",
    "\t    (X,y)->train_and_validate(X, y, holdout_partition(length(X), 8), fit, metric),\n",
    "\t]\n",
    "\n",
    "\tfor est_method in est_methods\n",
    "\t    vals = [est_method(X,y) for i in 1 : 100]\n",
    "\t    push!(mins, minimum(vals))\n",
    "\t    push!(maxs, maximum(vals))\n",
    "\t    push!(medians, median(vals))\n",
    "\t    push!(quartile1s, quantile(vals, 0.25))\n",
    "\t    push!(quartile3s, quantile(vals, 0.75))\n",
    "\tend\n",
    "\n",
    "\tcolors = [\"pastelPurple\", \"pastelBlue\", \"pastelSeaGreen\", \"pastelGreen\", \"pastelRed\"]\n",
    "\tplots = Plots.Plot[]\n",
    "\tfor i in 1:length(mins)\n",
    "\t    push!(plots, Plots.Command(\"\"\"\\\\addplot+[$(colors[i]),\n",
    "\t        boxplot prepared={\n",
    "\t          median=$(medians[i]),\n",
    "\t          upper quartile=$(quartile3s[i]),\n",
    "\t          lower quartile=$(quartile1s[i]),\n",
    "\t          upper whisker=$(maxs[i]),\n",
    "\t          lower whisker=$(mins[i])\n",
    "\t        }, ] coordinates {};\"\"\"\n",
    "\t    ))\n",
    "\tend\n",
    "\tAxis(plots, style=\"ytick={1,2,3,4,5}, yticklabels={0.632 Bootstrap, Leave-One-Out Bootstrap, Bootstrap, 5-Fold Cross Validation, Holdout}\", width=\"8cm\", height=\"4cm\", xlabel=\"root mean squared error\")\n",
    "end\n",
    "plot(p)"
   ]
  }
 ],
 "metadata": {
  "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
}