题意： 蒙特卡洛采样之拒绝采样解析： 给定一个概率分布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
Homework 4PSTAT 115, Fall 2019Due on November 3, 2019 at 11:59 pmNote: If you are working with a partner, please submit only one homework per group with both namesand whether you are taking the course for graduate credit or not. Submit your Rmarkdown (.Rmd) and thecompiled pdf on Gauchospace.1. Rejection Sampling the Beta distribution. Assume we did not have access to the rbeta function forsampling from a Beta, but we were able to evaluate the density, dbeta. This is a very common settingin 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 findM1 = maxxp(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 objectiveargument returned by optimize (maximum tells us where the maximum occurs, but not what heightit 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 youneed to find M2 = maxxp(x)/q2(x) as above. Propose 10000 samples and keep only the acceptedsamples.(d) Plot the p(x), M1q1(x) and M2q2(x) all on the same plot and verify visually that the scaledproposal densities “envelope” the target, p(x). Set the xlimits of the plot from 0 to 1. Use differentcolor 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 plotfrom the previous part? This means when proposing 10000 samples from each proposal, the MonteCarlo 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. Howdoes 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) = 14|sin(x)| × I{x ∈ [0, 2π]}Use a proposal density which is uniform from 0 to 2π and generate at least 1000 true samples fromp(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% HPDregion calculated on the samples. What are the bounds on the HPD region? Report the length ofthe 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 argumentpass in density(samples), where samples is the name of your vector of true samples from thedensity. Set the allowSplit argument to true and use the credMass argument to set the totalprobability 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 using1them segments function (base plotting) or geom_segements (ggplot). To ensure that the linesdon’t overlap visually, for the HPD region set the y-value of the segment to 0 and for the quantileinterval set the y-value to to 0.01. Make the segments for HPD region and the segment for quantileinterval different colors. Report the length of the quantile interval and the total length of the HPDregion, verifying that indeed the HPD region is smaller.2
