Chapter 2 Monte Carlo Simulation

Suppose \(\sigma^2\) is unknown. If \({\boldsymbol \beta}\sim \mathcal{N}_p({\boldsymbol 0},\sigma^2\nu{\bf I}_p)\), we know the distribution of \({\boldsymbol \beta}|{\bf y}\) is Eq.(??). According to this known distribution, we can easily compute the \(E({\boldsymbol \beta}|{\bf y})\) and \(V({\boldsymbol \beta}|{\bf y})\). But in practice, the form of density function \(\pi({\boldsymbol \beta}|{\bf y})\) is often an unknown and very complicated distribution, leading to be impossible to solve its integration directly. hence we turn to use numerical methods to solve this problem, such as Monte Carlo methods.

For example, a form density function \(\pi(\theta|{\bf y})\) is an unknown distribution. \(E(\theta|{\bf y}) = \int \theta\pi (\theta|{\bf y})d\theta\). According to lemma , \(\bar X_{n} \rightarrow \mu = E(X)\) as \(n \rightarrow \infty\). Therefore, if we indepedently generate \(\theta^{(1)}, \theta^{(2)}, \ldots,\theta^{(M)}\) from \(\pi(\theta|{\bf y})\), \(M^{-1} \sum_{k=1}^M\theta^{(k)} \rightarrow E(\theta|{\bf y})\) as M \(\rightarrow \infty\). However, there are two problems.

  1. What if we cannot generate sample from \(\pi(\theta|{\bf y})\)

  2. What if they are not identically and independently distributed.

The solutions to these two question are Monte Carlo (MC) method and Markov Chain Monte Carlo (MCMC).

For the first question, we can use importance sampling method.

Lemma 2.1 (Strong Law of Large Number)

Let \(X_{1}, X_{2},...\) be a sequence of i.i.d random variables, each having finite mean m. Then \[\Pr\left(\lim_{n\rightarrow\infty} \frac{1}{n}(X_{1}+X_{2}+...+X_{n} = m)\right) = 1.\]

2.1 \(\sigma^2\) is known

2.1.1 Importance Sampling

To estimate mean value of parameters, usually we randomly sample from \(\pi(\theta|{\bf y})\) and compute their mean value to estimate parameters. However, in practice, it is hard to generate sample directly from \(\pi(\theta|{\bf y})\), because we don’t know what speccific distribution it belongs to. hence we return to importance sampling method. Note that \[\begin{eqnarray*} E(\theta|{\bf y}) &=& \int \theta\pi(\theta|{\bf y})d\theta\\ &=& \int \theta \frac{\pi(\theta|{\bf y})}{g(\theta)}g(\theta)d\theta\\ &=& \int h(\theta)g(\theta)d\theta \\ &=& E\{h(\theta)\} \end{eqnarray*}\] where \(h(\theta)=\theta \frac{\pi(\theta|{\bf y})}{g(\theta)}, g(\theta)\) is a known pdf. \ To estimate the variance of parameters, similarly, \[\begin{eqnarray*} Var(\theta|{\bf y}) &=& E[\theta - E(\theta|{\bf y}) ]^2 = \int (\theta- E(\theta|{\bf y}))^2\pi(\theta|{\bf y}) d\theta\\ &=& \int (\theta- E(\theta|{\bf y}))^2 \frac{\pi(\theta|{\bf y})}{g(\theta)}g(\theta)d\theta\\ &=& \int h^{'}(\theta)g(\theta)d\theta \\ &=& E\{h^{'}(\theta)\} \end{eqnarray*}\] where \(h^{'}(\theta)=(\theta- E(\theta|{\bf y}))^2 \frac{\pi(\theta|{\bf y})}{g(\theta)},\) and \(g(\theta)\) is a known pdf.\ The importance sampling method can be implemented as follows:

Note that \[ \sum_{i=1}^M \frac{h(\theta_{i})}{M} \rightarrow E(h(\theta)) = E(\theta|{\bf y}).\quad as \quad M \rightarrow\infty \] Estimating variance by importance sampling method is similarly.

2.1.2 Importance Sampling Simulation

Setup

Consider a single linear regression model as follows:

\[y_i=\beta_0 +\beta_1 x_i+ \varepsilon_i\] where \(\varepsilon_i \sim N(0,1), x_i\sim N(0,1), \beta_0 = 1, \beta_1 = 2,\) for \(i=1,\ldots,100.\)

We already know that if \({\boldsymbol \beta}\sim \mathcal{N}_p({\boldsymbol 0},\sigma^2\nu{\bf I}_p)\), then the distribution of \({\boldsymbol \beta}|{\bf y}\) is Eq.(??). Assume \(\sigma = 1\) is known and consider a known pdf \(\pi({\boldsymbol \beta}) = \phi({\boldsymbol \beta};(0,0), 10I_2)\), hence in this case, \(\nu = 10\). Our simulation can be broken down into 3 steps in following:

  1. Find exact value of \(E(\beta_{0}|{\bf y}), E(\beta_{1}|{\bf y}), V(\beta_{0}|{\bf y})\) and \(V(\beta_{1}|{\bf y})\)

  2. Use the MC method to simulate the results above with size of 100, 1000 and 5000.

  3. Use Importance Sampling Method to simulate relevant results with the same size in (2). Consider the known pdf of parameters are \(\phi(\beta_0|1, 0.5)\) and \(\phi(\beta_1|2, 0.5).\)

Simulation Results

We use R softeware to do this simulation.

  1. According to Eq.(), compute \(E({\boldsymbol \beta}|{\bf y})=\left({\bf X}^{{\bf T}} {\bf X}+ \frac{1}{\nu}I_p \right)^{-1} {\bf X}^{{\bf T}} {\bf y}\) and \(Var({\boldsymbol \beta}|{\bf y})=\sigma^2 \left({\bf X}^{{\bf T}} {\bf X}+ \frac{1}{\nu}I_p \right)^{-1}\), where \(\sigma =1\) and \(\nu = 10\), then we get the following results:
  1. Sample directly from normal distribution of Eq.(??) with sample size of 100, 1000 and 5000, then compute the mean and variance values of samples. hence we get the following results:
  2. Sample \(\beta_0\) and \(\beta_1\) directly from \(\phi(\beta_0|1, 0.5)\) and \(\phi(\beta_1|2, 0.5)\), respectively, with sample size of 100, 1000 and 5000. And then compute \(h(\theta_i)\) and \(h^{'}(\theta_i)\). Compute mean values of them to get estimates of expectation and variance. The final results are following: