{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# M1399.000200 Homework 1\n",
    "\n",
    "#### Due Sep 26, 2021 @ 11:59PM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q1. Git/GitHub\n",
    "\n",
    " **No handwritten homework reports are accepted for this course.** We work     with Git and GitHub. Efficient and abundant use of Git, e.g., frequent and    well-documented commits, is an important criterion for grading your homework.\n",
    "\n",
    " 1. Apply for the [Student Developer Pack](https://education.github.com/pack)  at GitHub using your `snu.ac.kr` email.\n",
    "\n",
    " 1. A link to join the M1399.000200 Github Classroom and a link to create an individual Github repository for homework is provided in the eTL. First join  the classroom, and then create your own homework repo by accepting these two  invitations in turn.\n",
    "\n",
    " 1. For each homework, the teaching assistant will make a pull request. Merge  each pull request to your homework repo.\n",
    "\n",
    " 1. In this course, you are required to write your homework reports using IJulia. For each homework, you need to submit your Jupyter notebook (e.g, `hw1.ipynb`), html (e.g., `hw1.html`), along with all code and data that are necessary to reproduce the results.\n",
    " \n",
    " 1. Provide your answer *directly* on this notebook. Add your answer right below the question. If the question has subproblems, split the cell at the right location and insert cells for your answer. **Do not modify the questions.**\n",
    "\n",
    " 1. Maintain two branches `master` and `develop`. The `develop` branch will    be your main playground, the place where you develop solution (code) to       homework problems and write up report. The `master` branch will be your       presentation area. Submit your homework files (Jupyter notebook file `ipynb`, `html`  file converted from the notebook, all code and data sets to reproduce results)  in `master` branch.\n",
    "\n",
    " 1. Before each homework's due date, commit your **master** branch. The        teaching assistant and the instructor will check out your committed master    branch for grading. Commit time will be used as your submission time. That    means if you commit your Homework 1 submission after the deadline, penalty    points will be deducted for late submission according to the syllabus.\n",
    " \n",
    " 1. Read the [style guide](https://docs.julialang.org/en/v1/manual/style-guide/index.html) for Julia programming. Following rules in the style guide will be strictly enforced when grading: \"Write functions, not just scripts\", \"Avoid writing overly-specific types\", \"Handle excess argument diversity in the caller\", \"Append ! to names of functions that modify their arguments\",  \"Don't use unnecessary static parameters\", \"Avoid using floats for numeric literals in generic code when possible\". Also it is advised to follow \"Use naming conventions consistent with Julia base/\", \"Write functions with argument ordering similar to Julia Base\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q2. Gibbs sampler in R and Julia\n",
    "\n",
    "The task in this question is to create a Gibbs sampler for the density  \n",
    "$$\n",
    "f(x, y) = k x^2 exp(- x y^2 - y^2 + 2y - 4x), x > 0\n",
    "$$\n",
    "using the conditional distributions\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "  X | Y &\\sim& \\Gamma \\left( 3, \\frac{1}{y^2 + 4} \\right) \\quad \\text{(shape, scale)}\\\\\n",
    "  Y | X &\\sim& N \\left(\\frac{1}{1+x}, \\frac{1}{2(1+x)} \\right).\n",
    "\\end{eqnarray*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following R function `Rgibbs()`.\n",
    "\n",
    "```r\n",
    "library(Matrix)\n",
    "Rgibbs <- function(N, thin) {\n",
    "  mat <- matrix(0, nrow=N, ncol=2)\n",
    "  x <- y <- 0\n",
    "  for (i in 1:N) {\n",
    "    for (j in 1:thin) {\n",
    "      x <- rgamma(1, 3, y * y + 4) # 3rd arg is rate\n",
    "      y <- rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))\n",
    "    }\n",
    "    mat[i,] <- c(x, y)\n",
    "  }\n",
    "  mat\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) Explain what the above function does."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) Write R code that generate a sample from $f$ of size 10,000 with a thinning of 500 **and** measure the time. Your R code should be a one-line. Run your code in Julia using the `RCall.jl` package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) Write a Julia function `jgibbs()` that does the same compuation as `Rgibbs()`. (*Hint*. use `Distributions.jl`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) Run `jgibbs(100, 5)` and then repeat problem 2 for `jgibbs()`. How long does it take? Why do you run `jgibbs(100, 5)` before measuring the time?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5) Compare the speed of `Rgibbs()` and `jgibbs()`. Also compare the programming efforts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q3. JIT\n",
    "\n",
    "Consider Julia function\n",
    "```julia\n",
    "function g(k::Number)\n",
    "    for i in 1:10\n",
    "        k = 2k - 1\n",
    "    end\n",
    "    k\n",
    "end\n",
    "```\n",
    "1) Use `@code_llvm` to find the LLVM bytecode of compiled `g(2)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) Use `@code_llvm` to find the LLVM bytecode of compiled `g(2.0)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) Compare the bytecode from questions 1 and 2. What do you find?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) Repeat questions 1-3 with `@code_native`. What is the difference from `@code_llvm`?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q4. Floats\n",
    "\n",
    "Consider the followiig model for the floating-point representation:\n",
    "$$\n",
    "    \\pm 0.d_1d_2\\dotsb d_p \\times b^e\n",
    "$$\n",
    "with $e_{\\min} \\leq e \\leq e_{\\max}$. Assume that $d_1 \\neq 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) How many floating-point numbers are there? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) What is the smallest positive number?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) What is the smallest number larger than 1? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) What is the smallest number $X$, such that $X + 1 = X$? Assume the round to nearest, ties to even rule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5) Suppose $p=4$ and $b=2$ (and $e_{\\min}$ is very small and $e_{\\max}$ is very large). What is the next number after 20 in this number system?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q5. Sum of squares\n",
    "\n",
    "Typical statistical computing routines would include that for computing the sum of squares of $n$ real number:\n",
    "$$\n",
    "    \\sum_{i=1}^n(x_i - \\bar{x})^2 = \\sum_{i=1}^n x_i^2 - n\\bar{x}^2\n",
    "$$\n",
    "where $\\bar{x} = \\frac{1}{n}\\sum_{i=1}^n x_i$.\n",
    "\n",
    "Either expression in the above equation suggests an algoriihm. The LHS suggests \n",
    "```\n",
    "# Algorithm 1\n",
    "a = x[1]\n",
    "for i = 2,...,n {\n",
    "    a = x[i] + a \n",
    "}\n",
    "a = a/n\n",
    "b = (x[1] − a)^2 \n",
    "for i = 2,...,n {\n",
    "    b=(x[i] − a)^2 + b \n",
    "}\n",
    "```\n",
    "\n",
    "The RHS suggests\n",
    "```\n",
    "# Algorithm 2\n",
    "a = x[1] \n",
    "b = x[1]^2\n",
    "for i = 2,...,n {\n",
    "    a = x[i] + a\n",
    "    b = x[i]^2 + b \n",
    "}\n",
    "a = a/n\n",
    "b = b − n * a^2\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) Compare Algorithms 1 and 2. Which algoritm is more likely to end up with catastrophic cancellation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative, consider the following algorithm.\n",
    "```\n",
    "# Algorithm 3\n",
    "a = x[1]\n",
    "b = 0.0\n",
    "for i = 2,...,n {\n",
    "    d = (x[i] − a) / i \n",
    "    a = d + a \n",
    "    b = i * (i − 1) * d^2 + b\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) Show that Algorithm 3 correctly computes the sum of squares."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) Compare Algorithm 3 with Algorithms 1 and 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) Write a Julia function `wsos(x; w)` that computes the weighted sum of squares $\\sum_{i=1^n}w_i(x_i - \\bar{x})^2$ in a similar way to Algorithm 3. Here, `x` is an array of numbers, and `w` is an array of positive weights."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "153px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}