Monte-Carlo algorithms

Task brief

When analysing a methylation data, it may be needed to simulate values from mixtures of beta distributions. In this assessment, you are requested to simulate values from a mixture of beta distributions using accept/reject sampling and importance sampling.

Define X a random variable describing the frequency of methylation. \(X_1, \ldots, X_n\) are i.i.d. distributed according to \(f(x)\). Define a mixture of three beta distributions as following:

\[f(x) = \frac{1}{3} Be(\alpha_1=1,\beta_1=5) + \frac{1}{3} Be(\alpha_2=3,\beta_2=5) + \frac{1}{3} Be(\alpha_3=10,\beta_3=5)\]

1. Plot the density \(f(x)\).

2. Implement an accept/reject algorithm, with \(Unif(0,1)\) proposal distribution.

We first set the seed at 1234 in order to allow reproducability

We then define the value of K, which is the maximum of the density. From looking at the graph above of the mixed beta distribution, we can see that the function has a maximum at x=0. Therefore to find the maximum value of \(f(x)\), we solve \(f(x=0)\) and find the associated theta value is 1.67. We can alternatively use an \(R\) function max() which returns the same value.

We then use runif() to simulate 10,000 values from a uniform distribution.

The acceptance probability is defined as \[ p = \frac{\pi(\theta|\mathbf{x})}{K*g(\theta)} \] where \(\pi(\theta|\mathbf{x})\) is our mixed beta function f(x), \(K\) is the maximum of that function, and \(g(\theta)\) is the uniform distribution in (0, 1), so it is equal to 1 for all values of theta.

Therefore we can define our acceptance probability as simply f(theta)/K.

We then run our acceptance algorithm. First we define another variable u with 10000 values from another \(Unif(0, 1)\) distribution. Then we make a filter ind that identifies only the values in u which are less than the acceptance probability.

Finally we apply the filter to our proposed thetas, and plot the accepted thetas on a histogram, with the target distribution added for comparison.

3. Compute the observed acceptance rate and compare it with the theoretical one.

The observed acceptance rate is the proportion of points in proposed theta that are accepted, ie length(theta[ind])/N >> [1] 0.5978.

The theoretical acceptance rate is \(1/K\), because \(K*1\) represents the entire area of proposed thetas, while 1 is the area under the pdf. In our case \(1/K = 0.6\).

The theoretical and observed acceptance rates are nearly the same.

[1] 0.5978
[1] 0.6

4. Implement an importance sampling algorithm.

We use the same variable theta with 10000 values from \(Unif(0,1)\) for the implementation of the accept/reject algorithm.

The weights are computed by \[ w^{(i)}\propto\frac{\pi(\theta^{(i)}|\mathbf{x})}{g(\theta^{(i)})} \] In our case \(g(\theta) == 1\) for all \(\theta\) because we are using the \(Unif(0,1)\) distribution. So we can omit \(g(\theta)\) from our calculations.

We then need to normalise the weights to \((0, 1)\).

Finally, we can derive the density of the importance sample directly with density(), setting weights=W to apply the weights to each value of proposed thetas.

5. Compare on the same plot the target density, the density of the accepted values for the accept/reject algorithm, and the density of the values weighted by importance sampling.

6. Which of the two algorithms seems to approximate the target distribution better? Which one you would use and why?

Both algorithms show a very similar density. Unfortunately both algorithms deviate from the target distribution towards x=0, for an unknown reason.

To compare further I checked the mean of the target distribution by weighting the averages of each component distribution by \(1/3\). I found the mean = 0.4028.

\[\mathbb{E}[\theta|\mathbf{x}] = \frac {1}{3} (\frac{1}{1+5}) + \frac {1}{3}(\frac{3}{3+5}) + \frac {1}{3}(\frac{10}{10+5}) = 0.4028\] The mean of the rejection sampling distribution can be checked with mean(theta[ind]).

\[\mathbb{E}_R[\theta|\mathbf{x}] = 0.4018\]

Finally, the mean of the importance sampling distribution can be checked with

\[\mathbb{E}_I[\theta|\mathbf{x}] \leftarrow \frac {1}{N} {\sum w^{(i)}\theta^{(i)}} = 0.4056\]

We can see that the mean of the rejection sampling algorithm comes closer to the true value than the importance sampling algorithm, though only slightly.

[1] 0.4017751
[1] 0.4055834

I would use the rejection sampling algorithm because the mean is closer to the target value, and I find it to be more intuitive and elegant. However it does depend on our being able to compute \(K\), which can easily be done in this case, but may sometimes be difficult. If we can’t compute \(K\), or if we can’t find a distribution \(g(\theta)\) such that \(K*g(\theta) \geqslant \theta \ , \forall \ \theta\), then we would need to use importance sampling.