Lab 2

Attribution

The first exercise is from assignment 4 from Aki Vehtari’s Bayesian Data Analysis course at Aalto university. The course is here: Aalto BDA3 course The second exercise is adapted from a question in Robert & Casella.

Mean estimates and the Monte Carlo standard error (MCSE)

In this exercise, we’ll familiarize ourselves with the concept of the MCSE. As we’ll see, this allows us to determine how many decimal digits should be reported.

We will explore mean estimates and MCSE estimates of two distributions.

First is the gamma distribution, parameterised with shape ($\alpha$) and rate ($\beta$).

Question 1.2

The density function for a gamma distributed random variable (RV) with shape $\alpha$ and rate $\beta$ is defined as: \[ f_X(x) = \frac{\beta^\alpha}{\Gamma(\alpha)}x^{\alpha - 1} e^{\beta x} \] The mean of this distribution is $\alpha / \beta$.

Use the appropriate R function to draw a sample of size 1000 from a gamma($\alpha = 3$, $\beta = 3$) distribution. We know the mean of this distribution is exactly $1$.

Question 1.3

The brute-force way to calculate the Monte Carlo standard error of the mean is to repeatedly take samples of draws from the distribution, calculate the mean of each sample and then take the standard deviation of the means. Let’s use this method to calculate the sampling distribution of the mean of $400$ draws from the gamma($\alpha=3, \beta=3$) distribution.

Simulate 2000 samples, each of 400 draws from the gamma($\alpha = 3$, $\beta = 3$), calculate and save the mean of each sample. Then calculate the standard deviation of the sample of means.

Adjust the following R code to do this:

set.seed(4711)

NUMBER_OF_SAMPLES <- 1
NUMBER_OF_DRAWS <- 1

# create a vector to store the sample means
sample_means <- rep(NA, times = NUMBER_OF_SAMPLES)

# loop over the number of samples
for (i in 1:NUMBER_OF_SAMPLES) {
  # generate a sample from the gamma distribution of size 400
    #sample <- 
  # calculate the sample mean and save it to the vector
    #sample_means[i] <- mean(sample)
}

# calculate the standard deviation of the sample means
 sd(sample_means)

[1] NA

Question 1.5

Enter the MCSE that is calculated from the code and report your answer with 2 decimal digits.

Question 1.6

Now, let’s explore the Cauchy(location = 0, scale = 1) distribution.

What is the mean of the Cauchy(0, 1) distribution?

Change the simulation code you wrote above to take samples of the Cauchy(0, 1) distribution and calculate the MCSE of the mean. Try this a few times.

Question 1.7

What do you notice about the MCSE for the mean of the Cauchy and gamma distributions?

The Central Limit Theorem (CLT) gives us the conditions under which it is also possible to estimate the MCSE from a single sample for means.

Question 1.8

Let’s put the claims of the CLT to test. Assume you’ve taken $n$ independent draws of a random variable and stored them in a vector $x$. Assume also that the random variable follows a distribution with finite mean and variance. What is the correct formula for the MCSE of the mean? $\text{sd}$ denotes the standard deviation and $\text{var}$ the variance of $x$.

Question 1.10

MCSE can be calculated for any estimate, including quantiles and probabilities, and can be used to guide how many digits to report for a given estimate. Suppose you calculate the following summary of a parameter:

mean = 0.483834 (MCSE = 0.003)

Using a 99% confidence level, what is a reasonable number of digits to report the mean? How would you report the confidence interval for the mean?

Self-normalized importance sampling in practice

Recall that self-normalized importance sampling is defined as the estimator for $\ExpA{h(X)}{X \sim f}$:

\[ \tilde{\mu}_{\text{IS}} = \sum_{i=1}^S h(X_i) w(X_i), \, w(X_i) = \frac{f_u(X_i) / g(X_i)}{\sum_{j=1}^S f_u(X_j) / g(X_j)} \]

where $g$ is our importance density and $f_u(X_i)$ is an un-normalized probability density.

Let’s suppose our importance density is a multivariate $t$ distribution, which we can do in R by using the mvtnorm package, and our target density will be the following density:

\[ f_u(x) \propto \norm{x}_2^{1-d} \exp\lp -x^T \Sigma^{-1}x/2 \rp \] and we’re going to estimate the second moment for $X_1$ and the $X_1 X_2$ for $X \in \R^3$.

We’ll define $\Sigma$ in the code below.

set.seed(123)
d <- 2
nu <- 50
n_draws <- 100
Sigma <- diag(d)
Sigma[lower.tri(Sigma)] <- Sigma[upper.tri(Sigma)] <- 0.9

library(mvtnorm)

X <- rmvt(n = n_draws, df = nu, sigma = Sigma)
log_g_draws <- dmvt(X, sigma=Sigma, df = nu)

Your turn

Write a function that calculates $\log f_u(x)$, and apply it to the draws X. We’re going to use the same Sigma we defined above, but in order to make our code efficient, we should only pass in $\Sigma^{-1}$ so we don’t invert the matrix repeatedly.

log_f_u <- function(x, inv_sigma) {
 ## write your code 
}

log_f_u_draws <- rep(NA_real_,n_draws)

Now we need to calculate the importance ratios with the log densities we have.

The formula is:

\[ \log w(X_i) = \log f_u(X_i) - \log g(X_i) - \log\lp \sum_{j=1}^S f_u(X_j)/g(X_j) \rp \]

The first two terms are straightforward, but we need to calculate the log of a sum of the exponentiated log densities. We can use a function called log_sum_exp to do so:

## Code from Gemini
log_sum_exp <- function(x) {
  max_x <- max(x)
  return(max_x + log(sum(exp(x - max_x))))
}

To use this we pass in the vector log_f_u_draws - log_g_draws and we’ll get out the log of the weights denominator.

Your turn

Calculate the log-weights using log_f_u_draws and log_g_draws.

log_weights <- rep(NA_real_,n_draws)

Make a histogram of the weights. What do you notice?

hist(log_weights)

Remember that in order to have well-behaved weights, we want our importance density to have heavier tails than our target distribution.

Try changing nu above to 2 from 10 and see how the distribution of the log weights changes.

Now with these new log-weights we can approximate $\Exp{h(X)}{f}$ using $\tilde{\mu}_{\text{IS}}$.

In order to do so, we need to define a function $h(x)$. Remember, we’re trying to estimate $X_1^2$ and $X_1 X_2$. What is the dimension of the range of $h$?

h <- function(x) {
 ## Your code  
}

Report the estimates for $\ExpA{X_1^2}{X \sim g}$ and $\ExpA{X_1 X_2}{X \sim g}$.

--- title: "Lab 2" format: html: html-math-method: method: mathjax url: "https://cdn.jsdelivr.net/npm/mathjax@4.1/tex-mml-chtml.js" include-in-header: text: | <script> window.MathJax = { startup: { ready() { MathJax.startup.defaultReady(); const {STATE} = MathJax._.core.MathItem; MathJax.tex2mml(String.raw` \newcommand{\lp}{\left(} \newcommand{\rp}{\right)} \newcommand{\R}{\mathbb{R}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\lb}{\left[} \newcommand{\rb}{\right]} \newcommand{\ind}[1]{\mathbbm{1}\lp#1\rp} \newcommand{\Exp}[1]{\mathbb{E} \lb #1 \rb} \newcommand{\ExpA}[2]{\mathbb{E}_{#2} \lb #1 \rb} \newcommand{\Prob}[2]{\mathbb{P}_{#2} \lp #1 \rp} \newcommand{\indy}{\mathrel{\unicode{x2AEB}}} \newcommand{\comp}{^\mathsf{c}} \newcommand{\widebar}[1]{\overline{\! #1}} \newcommand{\Var}[1]{\mathrm{Var}\lp #1 \rp} \newcommand{\Cov}[2]{\mathrm{Cov}\lp #1,#2 \rp} \newcommand{\VarA}[2]{\mathrm{Var}_{#2}\lp #1 \rp} \newcommand{\norm}[1]{\left\lVert #1 \right\rVert} \newcommand{\ddtA}[1]{\frac{\partial^2}{\partial #1^2}} \newcommand{\ddtB}[2]{\frac{\partial^2}{\partial #1 \partial #2}} \newcommand{\diff}[1]{{\mathrm{d}}\!\lp #1 \rp} \newcommand{\veco}[1]{{\mathrm{vec}}\lp #1 \rp} \newcommand{\tr}[1]{{\mathrm{tr}}\lp #1 \rp} `); } }, loader: {load: ['[tex]/bbm','[tex]/upgreek']}, tex: {packages: {'[+]': ['bbm','upgreek']}} }; </script> code-tools: source: true toggle: false caption: "Download Source" --- # Attribution The first exercise is from assignment 4 from Aki Vehtari's Bayesian Data Analysis course at Aalto university. The course is here: [Aalto BDA3 course](https://avehtari.github.io/BDA_course_Aalto/Aalto2026.html) The second exercise is adapted from a question in Robert & Casella. # Mean estimates and the Monte Carlo standard error (MCSE) In this exercise, we'll familiarize ourselves with the concept of the MCSE. As we'll see, this allows us to determine how many decimal digits should be reported. We will explore mean estimates and MCSE estimates of two distributions. First is the gamma distribution, parameterised with shape ($\alpha$) and rate ($\beta$). ## Question 1.2 The density function for a gamma distributed random variable (RV) with shape $\alpha$ and rate $\beta$ is defined as: $$ f_X(x) = \frac{\beta^\alpha}{\Gamma(\alpha)}x^{\alpha - 1} e^{\beta x} $$ The mean of this distribution is $\alpha / \beta$. Use the appropriate R function to draw a sample of size 1000 from a gamma($\alpha = 3$, $\beta = 3$) distribution. We know the mean of this distribution is exactly $1$. ## Question 1.3 The brute-force way to calculate the Monte Carlo standard error of the mean is to repeatedly take samples of draws from the distribution, calculate the mean of each sample and then take the standard deviation of the means. Let's use this method to calculate the sampling distribution of the mean of $400$ draws from the gamma($\alpha=3, \beta=3$) distribution. Simulate 2000 samples, each of 400 draws from the gamma($\alpha = 3$, $\beta = 3$), calculate and save the mean of each sample. Then calculate the standard deviation of the sample of means. Adjust the following R code to do this: ```{r} set.seed(4711) NUMBER_OF_SAMPLES <- 1 NUMBER_OF_DRAWS <- 1 # create a vector to store the sample means sample_means <- rep(NA, times = NUMBER_OF_SAMPLES) # loop over the number of samples for (i in 1:NUMBER_OF_SAMPLES) { # generate a sample from the gamma distribution of size 400 #sample <- # calculate the sample mean and save it to the vector #sample_means[i] <- mean(sample) } # calculate the standard deviation of the sample means sd(sample_means) ``` ## Question 1.5 Enter the MCSE that is calculated from the code and report your answer with 2 decimal digits. ## Question 1.6 Now, let's explore the Cauchy(location = 0, scale = 1) distribution. What is the mean of the Cauchy(0, 1) distribution? Change the simulation code you wrote above to take samples of the Cauchy(0, 1) distribution and calculate the MCSE of the mean. Try this a few times. ## Question 1.7 What do you notice about the MCSE for the mean of the Cauchy and gamma distributions? The Central Limit Theorem (CLT) gives us the conditions under which it is also possible to estimate the MCSE from a single sample for means. ## Question 1.8 Let's put the claims of the CLT to test. Assume you've taken $n$ independent draws of a random variable and stored them in a vector $x$. Assume also that the random variable follows a distribution with finite mean and variance. What is the correct formula for the MCSE of the mean? $\text{sd}$ denotes the standard deviation and $\text{var}$ the variance of $x$. ## Question 1.10 MCSE can be calculated for any estimate, including quantiles and probabilities, and can be used to guide how many digits to report for a given estimate. Suppose you calculate the following summary of a parameter: - mean = 0.483834 (MCSE = 0.003) Using a 99% confidence level, what is a reasonable number of digits to report the mean? How would you report the confidence interval for the mean? # Self-normalized importance sampling in practice Recall that self-normalized importance sampling is defined as the estimator for $\ExpA{h(X)}{X \sim f}$: $$ \tilde{\mu}_{\text{IS}} = \sum_{i=1}^S h(X_i) w(X_i), \, w(X_i) = \frac{f_u(X_i) / g(X_i)}{\sum_{j=1}^S f_u(X_j) / g(X_j)} $$ where $g$ is our importance density and $f_u(X_i)$ is an un-normalized probability density. Let's suppose our importance density is a multivariate $t$ distribution, which we can do in R by using the `mvtnorm` package, and our target density will be the following density: $$ f_u(x) \propto \norm{x}_2^{1-d} \exp\lp -x^T \Sigma^{-1}x/2 \rp $$ and we're going to estimate the second moment for $X_1$ and the $X_1 X_2$ for $X \in \R^3$. We'll define $\Sigma$ in the code below. ```{r} set.seed(123) d <- 2 nu <- 50 n_draws <- 100 Sigma <- diag(d) Sigma[lower.tri(Sigma)] <- Sigma[upper.tri(Sigma)] <- 0.9 library(mvtnorm) X <- rmvt(n = n_draws, df = nu, sigma = Sigma) log_g_draws <- dmvt(X, sigma=Sigma, df = nu) ``` ## Your turn Write a function that calculates $\log f_u(x)$, and apply it to the draws `X`. We're going to use the same `Sigma` we defined above, but in order to make our code efficient, we should only pass in $\Sigma^{-1}$ so we don't invert the matrix repeatedly. ```{r} log_f_u <- function(x, inv_sigma) { ## write your code } log_f_u_draws <- rep(NA_real_,n_draws) ``` Now we need to calculate the importance ratios with the log densities we have. The formula is: $$ \log w(X_i) = \log f_u(X_i) - \log g(X_i) - \log\lp \sum_{j=1}^S f_u(X_j)/g(X_j) \rp $$ The first two terms are straightforward, but we need to calculate the log of a sum of the exponentiated log densities. We can use a function called `log_sum_exp` to do so: ```{r} ## Code from Gemini log_sum_exp <- function(x) { max_x <- max(x) return(max_x + log(sum(exp(x - max_x)))) } ``` To use this we pass in the vector `log_f_u_draws - log_g_draws` and we'll get out the log of the weights denominator. ## Your turn Calculate the log-weights using `log_f_u_draws` and `log_g_draws`. ```{r} log_weights <- rep(NA_real_,n_draws) ``` Make a histogram of the weights. What do you notice? ```{r} #| eval: false hist(log_weights) ``` Remember that in order to have well-behaved weights, we want our importance density to have heavier tails than our target distribution. Try changing `nu` above to 2 from 10 and see how the distribution of the log weights changes. Now with these new log-weights we can approximate $\Exp{h(X)}{f}$ using $\tilde{\mu}_{\text{IS}}$. In order to do so, we need to define a function $h(x)$. Remember, we're trying to estimate $X_1^2$ and $X_1 X_2$. What is the dimension of the range of $h$? ```{r} #| eval: false h <- function(x) { ## Your code } ``` Report the estimates for $\ExpA{X_1^2}{X \sim g}$ and $\ExpA{X_1 X_2}{X \sim g}$.