computer exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set run Julia Box parameter (set to true if you run the notebook on JuliaBox)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_on_juliabox = true "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load packages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Distributions\n",
    "using PyPlot\n",
    "using Random \n",
    "using StatsBase\n",
    "if run_on_juliabox; Pkg.add(\"KernelDensity\"); end \n",
    "using KernelDensity;\n",
    "using LinearAlgebra\n",
    "using Printf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main task: run ABC-RS for the g-and-k distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main task for this exercise is to use the ABC rejection-sampling algorithm to estimate the parameters for the g-and-k distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The g-and-k distribution\n",
    "\n",
    "The g-and-k distribution is a standard problem to analyze using ABC methods, and we will, therefore, first consider the problem of estimating the parameters of the g-and-k distribution using the ABC rejection sampling algorithm. \n",
    "\n",
    "For the g-and-k distribution is the quantile function give by \n",
    "\n",
    "\\begin{align}\n",
    "    F^{-1}(x \\mid A,B,c,g,k) = A + B \\Big( 1 + c \\frac{1- \\exp(-g\\cdot r(x))}{1+ \\exp(-g\\cdot r(x)} \\Big)(1+r^{2}(x))^{k}r(x).\n",
    "\\end{align}\n",
    "\n",
    "Where $A$,$B$,$c$,$g$, and $k$ are the parameters of the g-and-k distribution and $r(x)$ the $x$th standard normal quantile. \n",
    "\n",
    "We will for the g-and-k distribution follow a setting often used in the literatur (see [ [1]](#allingham_etal_09), and [ [3]](#picchini_17)), and we therefore assume that the unknown parameters are $\\theta = (A,B,g,k)$, and set $c = 0.8$. The true parameter value are set to $\\theta_{true} = (3,1,2,0.5)$, and the data set $y$ contains 5000 observations. We will follow [ [3]](#picchini_17) and use a flat prior distribution $\\mathcal{U}(0,10)$ for all parameters, and the summery statistics: $s(y) = (P_{20},P_{40},P_{50},P_{60},P_{80}, \\text{skew}(y))$, where $P_x$ is the $x$th percentile of the data and $\\text{skew}$ is the skewness of the data. We do not use any pre-run to estimate the weighs $w$ for the summery statistics, instead we use the same weighs as in [ [3]](#picchini_17) and set $w = (0.22,0.19,0.53,2.29,1.90)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to implemnt a ```data_generator``` that *fake* data for some set of parameters $\\theta = [A, B, c, g,k]$ with 1000 observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbr_obs = 1000 # We assume that we have 50 observations\n",
    "c = 0.8 # the c parameter is assumed to be known\n",
    "\n",
    "\"\"\"\n",
    "    data_generator(A::Real,B::Real,g::Real, k::Real)\n",
    "\n",
    "Inputs:\n",
    "- `A::Real` - A parameter\n",
    "- `B::Real` - A parameter\n",
    "- `g::Real` - A parameter\n",
    "- `k::Real` - A parameter\n",
    "\n",
    "Output \n",
    "- `F_inv` - `nbr_obs` samples from the g-and-k distribution\n",
    "\"\"\"\n",
    "function data_generator(A::Real,B::Real,g::Real, k::Real)\n",
    "\n",
    "    \n",
    "    # Hint: randn\n",
    "    \n",
    "    F_inv = zeors(nbr_obs) # pre-allocate output \n",
    "    \n",
    "    # Write your code here  \n",
    "  \n",
    "  return F_inv\n",
    "\n",
    "end "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will consider a simulation study, and we therefore first simulate the observed data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_true = 3\n",
    "B_true = 1\n",
    "g_true = 2\n",
    "k_true = 0.5\n",
    "\n",
    "Random.seed!(42)\n",
    "y_obs = data_generator(A_true, B_true, g_true, k_true)\n",
    "\n",
    "# plot simulated data\n",
    "PyPlot.figure(figsize=(10,4));\n",
    "PyPlot.subplot(121);\n",
    "PyPlot.plt[:hist](y_obs,100)\n",
    "PyPlot.ylabel(\"Freq.\");\n",
    "PyPlot.xlabel(\"Data\");\n",
    "PyPlot.title(\"Hist. of observed data\");\n",
    "\n",
    "# compute empercial CDF\n",
    "ecdf_func = ecdf(y_obs);\n",
    "x = collect(LinRange(0, 30, 1000));\n",
    "\n",
    "# plot emperical CDF\n",
    "PyPlot.subplot(122);\n",
    "PyPlot.plot(x,ecdf_func(x));\n",
    "PyPlot.ylabel(\"Percent\");\n",
    "PyPlot.xlabel(\"Data\");\n",
    "PyPlot.title(\"Empirical CDF\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have to implement the function ```S``` that compute the summary statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    S(y::Vector)\n",
    "\n",
    "Input:\n",
    "- `y::Vector` - the vectors with the g-and-k data set\n",
    "\n",
    "Output \n",
    "- `s::Vector` - the vector of summary statistics \n",
    "\"\"\"\n",
    "\n",
    "function S(y::Vector)\n",
    "\n",
    "    # Hint: ?percentile, ?skewness \n",
    "    \n",
    "    s = zeros(5) # pre-allocate output \n",
    "    \n",
    "    # Write your code here  \n",
    "    \n",
    "    \n",
    "    return s\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the summary statistics and convince yourself that the statistics are reasonable (i.e. compare the summary statistics for the observed data with the histogram or the empirical distribution function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S(y_obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we have to implement the function ```Δ``` that computes the distance between two sets of summary statistics. We will compute the distance as \n",
    "\n",
    "\\begin{align}\n",
    "    \\Delta(s^{\\star}, s, w) = (s^{\\star} - s)W^{-1}(s^{\\star} - s). \n",
    "\\end{align}\n",
    "\n",
    "Where $W$ is a diagonal matrix of the weighs $W = \\text{diag}\\{ w_1, \\cdots, w_n \\}$. The weights are obtained from a pre-run. The main advantage of computing the distance using this method is that we now can use summary statistics of different magnitudes without a problem. [ [5]](#prangle_17)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    Δ(s::Vector, s_star::Vector, w::Vector = [0.22; 0.19; 0.53; 2.97; 1.90])\n",
    "\n",
    "Input:\n",
    "- `s::Vector` - vector with summary statistics for a g-and-k data set\n",
    "- `s_star::Vector` - vector with summary statistics for a g-and-k data set\n",
    "- `w::Vector = [0.22; 0.19; 0.53; 2.97; 1.90]` - weigth vector (default input)\n",
    "\n",
    "Output \n",
    "- `dist::Real` - the distance \n",
    "\"\"\"\n",
    "function Δ(s::Vector, s_star::Vector, w::Vector = [0.22; 0.19; 0.53; 2.97; 1.90])\n",
    "    \n",
    "    # Hint: ?Diagonal, ?inv \n",
    "    \n",
    "    dist = NaN # pre-allocate output \n",
    "    \n",
    "    # write your code here \n",
    "    \n",
    "    return dist \n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try ```Δ(S(y_obs),S(y_obs))```. What should you get?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Δ(S(y_obs),S(y_obs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us run the prior predictive analyze, to check the priors before running the inference. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = 20\n",
    "\n",
    "A_prior_pred = rand(Uniform(0,10), samples)\n",
    "B_prior_pred = rand(Uniform(0,10), samples)\n",
    "g_prior_pred = rand(Uniform(0,10), samples)\n",
    "k_prior_pred = rand(Uniform(0,10), samples)\n",
    "\n",
    "y_prior_pred = zeros(length(y_obs), samples)\n",
    "\n",
    "\n",
    "for i in 1:samples; \n",
    "    y_prior_pred[:,i] = data_generator(A_prior_pred[i], B_prior_pred[i], g_prior_pred[i], k_prior_pred[i])\n",
    "end \n",
    "\n",
    "\n",
    "\n",
    "# compute empercial CDF\n",
    "x = collect(LinRange(0, 30, 1000));\n",
    "\n",
    "\n",
    "PyPlot.figure()\n",
    "for i in 1:samples\n",
    "    \n",
    "    ecdf_func = ecdf(y_prior_pred[:,i]);\n",
    "    PyPlot.plot(x,ecdf_func(x));\n",
    "    \n",
    "end \n",
    "ecdf_func = ecdf(y_obs);\n",
    "PyPlot.plot(x,ecdf_func(x), \"k\");\n",
    "\n",
    "PyPlot.ylabel(\"Percent\");\n",
    "PyPlot.xlabel(\"Data\");\n",
    "PyPlot.title(\"Prior predictive\");\n",
    "\n",
    "-- \n",
    "samples = 20\n",
    "\n",
    "# sample parameters from prior \n",
    "\n",
    "# write your code here\n",
    "\n",
    "y_prior_pred = zeros(length(y_obs), samples) # pre-allocate matrix with the g-and-k data sets\n",
    "\n",
    "\n",
    "# generate g-and-k data sets \n",
    "\n",
    "# write your code here\n",
    "\n",
    "\n",
    "# Plot data sets\n",
    "# compute empercial CDF\n",
    "x = collect(LinRange(0, 30, 1000));\n",
    "\n",
    "\n",
    "PyPlot.figure()\n",
    "for i in 1:samples\n",
    "    \n",
    "    ecdf_func = ecdf(y_prior_pred[:,i]);\n",
    "    PyPlot.plot(x,ecdf_func(x));\n",
    "    \n",
    "end \n",
    "ecdf_func = ecdf(y_obs);\n",
    "PyPlot.plot(x,ecdf_func(x), \"k\");\n",
    "\n",
    "PyPlot.ylabel(\"Percent\");\n",
    "PyPlot.xlabel(\"Data\");\n",
    "PyPlot.title(\"Prior predictive\");\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do you think the priors are reasonable?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the ABC-RS algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    abc_rs(N_proposals::Int, ϵ::Real, S::Function, Δ::Function)\n",
    "\n",
    "Input:\n",
    "- `N_proposals::Int` - number of proposals \n",
    "- `ϵ::Real` - ABC threashold\n",
    "- `S::Function` = function to compute the summary statistics\n",
    "- `Δ::Function` = function to compute the ABC distance \n",
    "\n",
    "Output \n",
    "- `abc_posterior_samples[:,1:nbr_accepted_proposals]` - accapted proposals  \n",
    "\"\"\"\n",
    "function abc_rs(N_proposals::Int, ϵ::Real, S::Function, Δ::Function)\n",
    "    \n",
    "    @printf \"Running ABC-RS\\n\"\n",
    "\n",
    "    abc_posterior_samples = zeros(4, N_proposals) # pre allocate output \n",
    "    nbr_accepted_proposals = 1 \n",
    "    \n",
    "    # write your code here\n",
    "    \n",
    "    print_interval = 10000\n",
    "    \n",
    "    for i in 1:N_proposals\n",
    "    \n",
    "        if mod(i,print_interval) == 0; @printf \"Percentage done: %.2f\\n\" 100*(i)/N_proposals; end\n",
    "\n",
    "        # write your code here \n",
    "        \n",
    "    end \n",
    "    \n",
    "    return abc_posterior_samples[:,1:nbr_accepted_proposals]\n",
    "end; "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run ABC-RS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Random.seed!(42) # fix random numbers\n",
    "N_proposals = 2*10^5\n",
    "ϵ=2\n",
    "abc_posterior_samples = @time abc_rs(N_proposals,ϵ, S, Δ)\n",
    "@printf \"Acceptance rate: %.2f %%\" length(abc_posterior_samples)/10^6*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot posterior inference results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calc grid for prior dist\n",
    "x_grid = collect(LinRange(-0.5, 10.5, 100));\n",
    "\n",
    "# calc prior dist\n",
    "priordens = pdf.(Uniform(0, 10), x_grid);\n",
    "\n",
    "# calc kernel density approx. for approx. posterior \n",
    "A_post_kde = kde(abc_posterior_samples[1,:]);\n",
    "B_post_kde = kde(abc_posterior_samples[2,:]);\n",
    "g_post_kde = kde(abc_posterior_samples[3,:]);\n",
    "k_post_kde = kde(abc_posterior_samples[4,:]);\n",
    "\n",
    "\n",
    "PyPlot.figure(figsize=(10,8));\n",
    "PyPlot.subplot(221);\n",
    "PyPlot.plot(A_post_kde.x,A_post_kde.density, \"b\");\n",
    "PyPlot.plot(x_grid,priordens, \"g\");\n",
    "PyPlot.plot((A_true, A_true), (0, maximum(A_post_kde.density)), \"k\");\n",
    "PyPlot.xlabel(L\"$A$\");\n",
    "PyPlot.subplot(222);\n",
    "PyPlot.plot(B_post_kde.x,B_post_kde.density, \"b\");\n",
    "PyPlot.plot(x_grid,priordens, \"g\");\n",
    "PyPlot.plot((B_true, B_true), (0, maximum(B_post_kde.density)), \"k\");\n",
    "PyPlot.xlabel(L\"$B$\");\n",
    "PyPlot.subplot(223);\n",
    "PyPlot.plot(g_post_kde.x,g_post_kde.density, \"b\");\n",
    "PyPlot.plot(x_grid,priordens, \"g\");\n",
    "PyPlot.plot((g_true, g_true), (0, maximum(g_post_kde.density)), \"k\");\n",
    "PyPlot.xlabel(L\"$g$\");\n",
    "PyPlot.subplot(224);\n",
    "PyPlot.plot(k_post_kde.x,k_post_kde.density, \"b\");\n",
    "PyPlot.plot(x_grid,priordens, \"g\");\n",
    "PyPlot.plot((k_true, k_true), (0, maximum(k_post_kde.density)), \"k\");\n",
    "PyPlot.xlabel(L\"$k$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conclusion? What about the $g$ parameter?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the posterior predictive checks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = 20\n",
    "\n",
    "posterior_idx = rand(1:size(abc_posterior_samples,2), samples)\n",
    "\n",
    "A_prior_pred = abc_posterior_samples[1,posterior_idx]\n",
    "B_prior_pred = abc_posterior_samples[2,posterior_idx]\n",
    "g_prior_pred = abc_posterior_samples[3,posterior_idx]\n",
    "k_prior_pred = abc_posterior_samples[4,posterior_idx]\n",
    "\n",
    "y_prior_pred = zeros(length(y_obs), samples)\n",
    "\n",
    "\n",
    "for i in 1:samples; \n",
    "    y_prior_pred[:,i] = data_generator(A_prior_pred[i], B_prior_pred[i], g_prior_pred[i], k_prior_pred[i])\n",
    "end \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# compute empercial CDF\n",
    "\n",
    "x = collect(LinRange(0, 30, 1000));\n",
    "\n",
    "\n",
    "PyPlot.figure()\n",
    "\n",
    "for i in 1:samples\n",
    "    \n",
    "    ecdf_func = ecdf(y_prior_pred[:,i]);\n",
    "    PyPlot.plot(x,ecdf_func(x));\n",
    "    \n",
    "end \n",
    "\n",
    "ecdf_func = ecdf(y_obs);\n",
    "PyPlot.plot(x,ecdf_func(x), \"k\");\n",
    "\n",
    "PyPlot.ylabel(\"Percent\");\n",
    "PyPlot.xlabel(\"Data\");\n",
    "PyPlot.title(\"Posterior predictive\");\n",
    "\n",
    "--\n",
    "\n",
    "samples = 20\n",
    "\n",
    "# sample parameters from posterior \n",
    "\n",
    "y_prior_pred = zeros(length(y_obs), samples)\n",
    "\n",
    "# pre-allocate matrix with the g-and-k data sets\n",
    "\n",
    "\n",
    "# generate g-and-k data sets \n",
    "\n",
    "# write your code here\n",
    "\n",
    "\n",
    "\n",
    "# compute empercial CDF\n",
    "\n",
    "x = collect(LinRange(0, 30, 1000));\n",
    "\n",
    "\n",
    "PyPlot.figure()\n",
    "\n",
    "for i in 1:samples\n",
    "    \n",
    "    ecdf_func = ecdf(y_prior_pred[:,i]);\n",
    "    PyPlot.plot(x,ecdf_func(x));\n",
    "    \n",
    "end \n",
    "\n",
    "ecdf_func = ecdf(y_obs);\n",
    "PyPlot.plot(x,ecdf_func(x), \"k\");\n",
    "\n",
    "PyPlot.ylabel(\"Percent\");\n",
    "PyPlot.xlabel(\"Data\");\n",
    "PyPlot.title(\"Posterior predictive\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conclusions?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional tasks\n",
    "\n",
    "We can conclude that the ABC-RS algorithm does not work very well for the g-and-k distribution. Instead you can try implementing the ABC-MCMC algorithm (see  [ [2]](#marjoram_03), and [ [6]](#scott_10)). Try to replicate the g-and-k example in  [ [3]](#picchini_17) using ABC-MCMC.\n",
    "\n",
    "Another idea is to compare the ABC posterior with the exact analytical posterior, that indeed is available for the g-and-k method. How to compute the analytical posterior for the g-and-k distribution is explained in [[4]](#prangle_17_gk).\n",
    "\n",
    "The solutions for the additional tasks can be found here https://github.com/SamuelWiqvist/introlikelihoodfree (However please note that the binder support does not work)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References \n",
    "\n",
    "1. Allingham, David, R. A. R. King, and Kerrie L. Mengersen. \"Bayesian estimation of quantile distributions.\" Statistics and Computing 19, no. 2 (2009): 189-201. <a name=\"allingham_etal_09\"></a>\n",
    "\n",
    "2. Marjoram, Paul, John Molitor, Vincent Plagnol, and Simon Tavaré. \"Markov chain Monte Carlo without likelihoods.\" Proceedings of the National Academy of Sciences 100, no. 26 (2003). <a name=\"marjoram_03\"></a>\n",
    "\n",
    "3. Picchini, Umberto, and Rachele Anderson. \"Approximate maximum likelihood estimation using data-cloning ABC.\" Computational Statistics & Data Analysis 105 (2017): 166-183. <a name=\"picchini_17\"></a>\n",
    "\n",
    "4. Prangle, Dennis. \"gk: An R Package for the g-and-k and generalised g-and-h Distributions.\" arXiv preprint arXiv:1706.06889 (2017). <a name=\"prangle_17_gk\"></a>\n",
    "\n",
    "5. Prangle, Dennis. \"Adapting the ABC distance function.\" Bayesian Analysis 12.1 (2017): 289-309. <a name=\"prangle_17\"></a>\n",
    "\n",
    "6. Scott A. Sisson, Yanan Fan. \"Likelihood-free Markov chain Monte Carlo.\" arXiv preprint arXiv:1001.2058v1 (2010). <a name=\"scott_10\"></a>"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.0.0",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}