PSTAT 115 Homework4 课业解析

题意：

蒙特卡洛采样之拒绝采样

解析：

给定一个概率分布p(z)=p~(z)/Zp,p~(z)已知，Zp为归一化常数，为未知数。对该分布进行拒绝采样，我们引入一个简单地参考分布，记作q(x)，q(x)分布的采样是易于实现的，比如均匀分布。再引入一个常数k，满足kq(z)>p~(z)。每次采样中首先从q(z)采样一个数值z0，然后在区间[0,kq(z0)]进行均匀采样，得到u0。如果u0p~(z)的基础上尽可能小。

涉及知识点：

拒绝采样

更多可+薇❤讨论：Rainbow890722

pdf

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.