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.
Apply for the Student Developer Pack at GitHub using your snu.ac.kr email.
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.
For each homework, the teaching assistant will make a pull request. Merge each pull request to your homework repo.
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.
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.
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.
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.
Read the style guide 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".
The task in this question is to create a Gibbs sampler for the density
$$
f(x, y) = k x^2 exp(- x y^2 - y^2 + 2y - 4x), x > 0
$$
using the conditional distributions
$$
\begin{eqnarray*}
X | Y &\sim& \Gamma \left( 3, \frac{1}{y^2 + 4} \right) \quad \text{(shape, scale)}\\
Y | X &\sim& N \left(\frac{1}{1+x}, \frac{1}{2(1+x)} \right).
\end{eqnarray*}
$$
Consider the following R function Rgibbs().
library(Matrix)
Rgibbs <- function(N, thin) {
mat <- matrix(0, nrow=N, ncol=2)
x <- y <- 0
for (i in 1:N) {
for (j in 1:thin) {
x <- rgamma(1, 3, y * y + 4) # 3rd arg is rate
y <- rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))
}
mat[i,] <- c(x, y)
}
mat
}
1) Explain what the above function does.
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.
3) Write a Julia function jgibbs() that does the same compuation as Rgibbs(). (Hint. use Distributions.jl)
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?
5) Compare the speed of Rgibbs() and jgibbs(). Also compare the programming efforts.
Consider Julia function
function g(k::Number)
for i in 1:10
k = 2k - 1
end
k
end
1) Use @code_llvm to find the LLVM bytecode of compiled g(2).
2) Use @code_llvm to find the LLVM bytecode of compiled g(2.0).
3) Compare the bytecode from questions 1 and 2. What do you find?
4) Repeat questions 1-3 with @code_native. What is the difference from @code_llvm?
Consider the followiig model for the floating-point representation: $$ \pm 0.d_1d_2\dotsb d_p \times b^e $$ with $e_{\min} \leq e \leq e_{\max}$. Assume that $d_1 \neq 0$.
1) How many floating-point numbers are there?
2) What is the smallest positive number?
3) What is the smallest number larger than 1?
4) What is the smallest number $X$, such that $X + 1 = X$? Assume the round to nearest, ties to even rule.
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?
Typical statistical computing routines would include that for computing the sum of squares of $n$ real number: $$ \sum_{i=1}^n(x_i - \bar{x})^2 = \sum_{i=1}^n x_i^2 - n\bar{x}^2 $$ where $\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i$.
Either expression in the above equation suggests an algoriihm. The LHS suggests
# Algorithm 1
a = x[1]
for i = 2,...,n {
a = x[i] + a
}
a = a/n
b = (x[1] − a)^2
for i = 2,...,n {
b=(x[i] − a)^2 + b
}The RHS suggests
# Algorithm 2
a = x[1]
b = x[1]^2
for i = 2,...,n {
a = x[i] + a
b = x[i]^2 + b
}
a = a/n
b = b − n * a^21) Compare Algorithms 1 and 2. Which algoritm is more likely to end up with catastrophic cancellation?
As an alternative, consider the following algorithm.
# Algorithm 3
a = x[1]
b = 0.0
for i = 2,...,n {
d = (x[i] − a) / i
a = d + a
b = i * (i − 1) * d^2 + b
}2) Show that Algorithm 3 correctly computes the sum of squares.
3) Compare Algorithm 3 with Algorithms 1 and 2.
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.