题意：
蒙特卡洛采样之拒绝采样
解析：
给定一个概率分布p(z)=p~(z)/Zp,p~(z)已知，Zp为归一化常数，为未知数。对该分布进行拒绝采样，我们引入一个简单地参考分布，记作q(x)，q(x)分布的采样是易于实现的，比如均匀分布。再引入一个常数k，满足kq(z)>p~(z)。每次采样中首先从q(z)采样一个数值z0，然后在区间[0,kq(z0)]进行均匀采样，得到u0。如果u0<p~(z0),则保留该采样值，否则丢弃该采样值。最后得到的数据就是一个对该分布的近似采样。为了提高接受效率，防止舍弃过多的采样值而导致采样效率低下，k值应该满足在kq(z)>p~(z)的基础上尽可能小。
涉及知识点：
拒绝采样
更多可加微信讨论
微信号：aaaapple1208
Homework 4
PSTAT 115, Fall 2019
Due on November 3, 2019 at 11:59 pm
Note: If you are working with a partner, please submit only one homework per group with both names
and whether you are taking the course for graduate credit or not. Submit your Rmarkdown (.Rmd) and the
compiled pdf on Gauchospace.
1. Rejection Sampling the Beta distribution. Assume we did not have access to the rbeta function for
sampling from a Beta, but we were able to evaluate the density, dbeta. This is a very common setting
in Bayesian statistics, since we can always evaluate the (proportional) posterior density p(θ | y) ∝ p(y |
θ)p(θ) but we don’t have immediate access to a method for sampling from this distribution.
(a) Let p(x) be a Beta(3, 9) density, q1(x) a Uniform(0, 1) density, and q2(x) a Normal(µ = 0.25, σ =
0.15) density.
(b) Use rejection sampling to sample from p(x) by proposing samples from q1(x). To do so, first find
M1 = max
x
p(x)/q1(x) using the optimize function and set lower=0, upper=1, and maximum =
TRUE (since we are maximizing not minimizing, the default). M will be the value in the objective
argument returned by optimize (maximum tells us where the maximum occurs, but not what height
it achieves). Propose 10000 samples and keep only the accepted samples.
(c) Use rejection sampling to sample from p(x) by proposing samples from q2(x). To do this you
need to find M2 = max
x
p(x)/q2(x) as above. Propose 10000 samples and keep only the accepted
samples.
(d) Plot the p(x), M1q1(x) and M2q2(x) all on the same plot and verify visually that the scaled
proposal densities “envelope” the target, p(x). Set the xlimits of the plot from 0 to 1. Use different
color lines for the various densities so are clearly distinguishable.
(e) Which rejection sampler had the higher rejection rate? Why does this make sense given the plot
from the previous part? This means when proposing 10000 samples from each proposal, the Monte
Carlo error of our approximation will be higher when proposing from ____ (choose q1 or q2).
(f) Report the variance of Beta(3, 9) distribution by computing the variance of the beta samples. How
does this compare to the theoretical variance (refer to the probability cheatsheet).
2. Interval estimation with rejection sampling.
(a) Use rejection sampling to sample from the following density:
p(x) = 1
4
|sin(x)| × I{x ∈ [0, 2π]}
Use a proposal density which is uniform from 0 to 2π and generate at least 1000 true samples from
p(x). Compute and report the Monte Carlo estimate of the upper and lower bound for the 50%
quantile interval using the quantile function on your samples. Compare this to the 50% HPD
region calculated on the samples. What are the bounds on the HPD region? Report the length of
the quantile interval and the total length of the HPD region. What explains the difference? Hint:
to compute the HPD use the hdi function from the HDInterval package. As the first argument
pass in density(samples), where samples is the name of your vector of true samples from the
density. Set the allowSplit argument to true and use the credMass argument to set the total
probability mass in the HPD region to 50%.
(b) Plot p(x) using the curve function (base plotting) or stat_function (ggplot). Add lines corresponding to the intervals / probability regions computed in the previous part to your plot using
1
them segments function (base plotting) or geom_segements (ggplot). To ensure that the lines
don’t overlap visually, for the HPD region set the y-value of the segment to 0 and for the quantile
interval set the y-value to to 0.01. Make the segments for HPD region and the segment for quantile
interval different colors. Report the length of the quantile interval and the total length of the HPD
region, verifying that indeed the HPD region is smaller.
2
