notebook/ch02.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Pkg\n",
    "if isfile(\"../Project.toml\") && isfile(\"../Manifest.toml\")\n",
    "    Pkg.activate(\"..\");\n",
    "end\n",
    "\n",
    "using Random\n",
    "using Distributions\n",
    "using Plots\n",
    "using StatsPlots\n",
    "using prml\n",
    "using prml: pdf\n",
    "Random.seed!(1234);\n",
    "gr();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bernoulli Distributions\n",
    "\n",
    "$$\n",
    "\\mathrm{Bern}(x \\mid \\mu) = x^{\\mu} \\cdot (1 - x)^{1 - \\mu}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = [0.0, 1.0, 1.0, 1.0];\n",
    "bernoulli = BernoulliDist(mu);\n",
    "pdf(bernoulli, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beta Distributions\n",
    "\n",
    "$$\n",
    "\\mathrm{Beta}(\\mu \\mid a, b) = \\dfrac{\\Gamma(a+b)}{\\Gamma(a) \\Gamma(b)} \\mu^{a-1}(1-\\mu)^{b-1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plots = [];\n",
    "l = @layout [a b; c d]\n",
    "\n",
    "x = collect(range(0, stop=1.0, length=100));\n",
    "for (index, values) in enumerate([[0.1, 0.1], [1, 1], [2, 3], [8, 4]])\n",
    "    a, b = values[1], values[2];\n",
    "    beta = BetaDist(a, b);\n",
    "    y = pdf(beta, x);\n",
    "    p = plot(x, y, xlims=(0, 1), ylims=(0, 3), label=\"a=$a, b=$b\", lw=2);\n",
    "    push!(plots, p);\n",
    "end\n",
    "\n",
    "plot(plots[1], plots[2], plots[3], plots[4], layout=l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/ch2/image1.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ML and MAP estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plots = [];\n",
    "l = @layout [a; b]\n",
    "\n",
    "x = collect(range(0, stop=1.0, length=100));\n",
    "\n",
    "beta = BetaDist(2, 2);\n",
    "p = plot(x, pdf(beta, x), xlim=(0, 1), ylim=(0, 2), label=\"prior\", lw=2);\n",
    "push!(plots, p);\n",
    "\n",
    "bern = BernoulliDist(beta);\n",
    "fitting(bern, [1.0])\n",
    "p = plot(x, pdf(bern._beta, x), xlim=(0, 1), ylim=(0, 2), label=\"posterior\", lw=2);\n",
    "push!(plots, p);\n",
    "\n",
    "plot(plots[1], plots[2], layout=l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/ch2/image2.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "println(\"Maximum likehood estimation\");\n",
    "\n",
    "model = BernoulliDist([1.0]);\n",
    "count = draw(model, 10000);\n",
    "println(\"$(count) out of 10000 is label-1\");\n",
    "\n",
    "println(\"Beysian estimation\");\n",
    "model = BernoulliDist(BetaDist(1.0, 1.0));\n",
    "fitting(model, [1.0]);\n",
    "count = draw(model, 10000);\n",
    "println(\"$(count) out of 10000 is label-1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Categorical variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = CategoricalDist([0.0]);\n",
    "fitting(model, [[0. 1. 0.]; [1. 0. 0.]; [0. 1. 0.]]);\n",
    "println(model._mu);\n",
    "\n",
    "mu = DirichletDist(ones(3));\n",
    "model = CategoricalDist(mu);\n",
    "println(\"prior\");\n",
    "println(model._dirichlet._alpha);\n",
    "trials = [[1. 0. 0.]; [1. 0. 0.]; [0. 1. 0.]];\n",
    "fitting(model, trials);\n",
    "println(\"posterior\");\n",
    "println(model._dirichlet._alpha);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dirichlet distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dirichlet1 = DirichletDist([0.1, 0.1, 0.1]);\n",
    "dirichlet2 = DirichletDist([1., 1., 1.]);\n",
    "dirichlet3 = DirichletDist([10., 10., 10.]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Gaussian Distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Central Limit Theorem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_sample = 10000;\n",
    "plots = [];\n",
    "for N in [1, 2, 10]\n",
    "    # a matrix of size (n_samples, N)\n",
    "    samples = rand(Uniform(0, 1), n_sample, N)\n",
    "    mean_samples = mean(samples, dims=2);\n",
    "    p =　StatsPlots.histogram(mean_samples, bins=20, legend=:topright);\n",
    "    push!(plots, p);\n",
    "end\n",
    "\n",
    "plot(plots[1], plots[2], plots[3], layout=(1, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/ch2/image3.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MAP estimation for the mean of Gaussian distribution\n",
    "\n",
    "For a Gaussian distribution $x \\sim \\mathcal{N}(\\mu, \\sigma^2)$ with $\\mu \\sim \\mathcal{N}(\\mu_0, \\sigma^{2}_0)$, given $N$ observations of $x$, denoted as  $\\boldsymbol{x} = [x_0, x_1, \\cdots, x_N]$, the posterior distribution of $\\mu$ is given by $\\mu \\sim \\mathcal{N}(\\mu_N, \\sigma^{2}_N)$ where\n",
    "\n",
    "\n",
    "$$\\begin{align} \\mu_N &= \\dfrac{\\sigma^{2}_N}{N\\sigma_{0}^2 + \\sigma^{2}_N}\\mu_0 + \\dfrac{N \\sigma_{0}^2}{N\\sigma_{0}^2 + \\sigma^{2}_N}\\mu_{\\mathrm{ML}} \\\\ \\dfrac{1}{\\sigma^{2}_N} &= \\dfrac{1}{\\sigma_{0}^2} + \\dfrac{N}{\\sigma^2}\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_prior = GaussianDist(0.0, 0.1);\n",
    "model = GaussianBayesMeanDist(mu_prior, 0.1);\n",
    "\n",
    "x = collect((range(-1, stop=1.0, length=200)));\n",
    "plot(x, pdf(model._mu_gauss, x), label=\"N=0\", xlim=(-1, 1), ylim=(0, 5), lw=2, color=\"steelblue\")\n",
    "\n",
    "fitting(model, [rand(Normal(0.8, 0.01))]);\n",
    "plot!(x, pdf(model._mu_gauss, x), label=\"N=1\", xlim=(-1, 1), ylim=(0, 5), lw=2, color=\"orange\")\n",
    "\n",
    "fitting(model, [rand(Normal(0.8, 0.01))]);\n",
    "plot!(x, pdf(model._mu_gauss, x), label=\"N=2\", xlim=(-1, 1), ylim=(0, 5), lw=2, color=\"green\");\n",
    "\n",
    "fitting(model, rand(Normal(0.8, 0.01), 8));\n",
    "plot!(x, pdf(model._mu_gauss, x), label=\"N=10\", xlim=(-1, 1), ylim=(0, 5), lw=2, color=\"red\", legend=:topleft)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gamma distribution\n",
    "\n",
    "$$\n",
    "\\begin{align} \\mathrm{Gam}(\\lambda \\mid a, b) = (b^a / \\Gamma(a)) \\lambda^{a-1} \\exp(-b\\lambda) \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = collect(range(0, stop=2.0, length=100));\n",
    "plots = [];\n",
    "for (index, values) in enumerate([[0.1, 0.1], [1., 1.], [2., 3.], [4., 6.]])\n",
    "    gamma = GammaDist(values[1], values[2]);\n",
    "    p = plot(x, pdf(gamma, x), xlim=(0, 2), ylim=(0, 2), label=\"a=$(values[1]), b=$(values[2])\", lw=2)\n",
    "    push!(plots, p);\n",
    "end\n",
    "\n",
    "plot(plots[1], plots[2], plots[3], plots[4], layout=(2, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/ch2/image5.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MAP estimation for the variance of Gaussian distribution\n",
    "\n",
    "For a Gaussian distribution $x \\sim \\mathcal{N}(\\mu, \\lambda = \\sigma^2)$ with $\\lambda \\sim \\mathrm{Gam}(\\lambda \\mid a, b)$, given $N$ observations of $x$, denoted as  $\\boldsymbol{x} = [x_0, x_1, \\cdots, x_N]$, the posterior distribution of $\\lambda$ is given by $\\lambda \\sim \\mathrm{Gam}(\\lambda \\mid a_N, b_N)$ where\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{align} a_N &=\\dfrac{N}{2} + a  = \\dfrac{N}{2} + \\dfrac{2a}{2} \\\\ b_N &= b + \\dfrac{N}{2} \\sigma^2_{\\mathrm{ML}}  = \\dfrac{2a}{2}\\dfrac{b}{a} + \\dfrac{N}{2} \\sigma^2_{\\mathrm{ML}} \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = GammaDist(1.0, 1.0);\n",
    "model = GaussianBayesVarDist(tau);\n",
    "println(\"a=$(model._var_gamma._a), b=$(model._var_gamma._b)\");\n",
    "\n",
    "fitting(model, rand(Normal(0, 1.414), 100));\n",
    "\n",
    "println(\"a=$(model._var_gamma._a), b=$(model._var_gamma._b)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate Gaussian Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = rand(Uniform(1, 3), 2, 100);\n",
    "gaussian = MultivariateGaussianDist(2);\n",
    "fitting(gaussian, X);\n",
    "\n",
    "x = range(-5, stop=10, length=100);\n",
    "y = range(-5, stop=10, length=100);\n",
    "grids = [[i, j] for i in x, j in y];\n",
    "vals = [pdf(gaussian, grids[i, j]) for i in 1:100, j in 1:100];\n",
    "contour(x, y, vals)\n",
    "plot!(X[1,:], X[2, :], seriestype=:scatter, legend=nothing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/ch2/image6.png\">"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.0",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}