R语言代做编程辅导和解答M3S9/M4S9 Stochastic Simulation: Project 2

1. Consider the following density:

f(x) / ( 0 otherwise. x(11-x) exp h- 12 -2 + ln 1-xx2

(a) Devise and implement two efficient algorithms for simulating from f(x).
(b) Estimate the normalizing constant using Monte Carlo integration.
(c) Devise and implement a Metropolis-Hastings sampler for generating variates
from f(x). In particular:
i) You should tune the Metropolis-Hastings algorithm to have acceptance rate
about 20%.
ii) Examine how the rate at which the algorithm reaches equilibrium depends
on the starting value.
iii) Consider carefully the correlation structure of the sequence generated.
iv) Compare the results of the Metropolis-Hastings sampler with the method
implemented in (a).

2. Consider the following bimodal \two-humps" density:

f(x; λ0) / exp -x2 1 + x2 2 + (x1 +2x2)2 - 2λ0x1x2 ; x 2 R2

for some parameter λ0, say λ0 = -4.

(a) Devise and implement a Metropolis-with-Gibbs sampler for generating variates from f(x; λ0).
(b) Devise and implement a Metropolis-Hastings sampler for generating variates from f(x; λ0).
(c) Compare the behavior of the Metropolis-with-Gibbs sampler and MetropolisHastings algorithm when λ0 = -4 and when λ0 = -8.

(a)

 

h=function(x)

{

  options(warn=-1)

  if(x>0 && x<1)v=exp(-((3+log(x/(1-x)))^2)/2)/(x-x^2)

  else v=0

normalfactor =function(n)

  {

ff=function(x){sqrt(f(x))}

fff<-function(x){x*sqrt(f(x))}

 

opt=function (n){#alpha,beta,theta are calculated using optimize function in R

   

    alpha = optimize(ff,c(0,1),maximum=T)$objective

    beta = 0 

    theta = optimize(fff,c(0,1),maximum=T)$objective

    tp <- (nf)/(2 *alpha * (theta - beta))

    factor = 1/((nf)/(2 * alpha * (theta - beta)))

输出前100000个分布的值

 

#envelop function env

env =function(x)

  {

    if(x<=0)v=0

    else if(x<=0.01)v=330*x

    else if(x<=0.03)v=33

黑色代表函数值，绿色代表envelop function的拟合值。

计算envelop function的累计密度函数

mv=optimize(f(x)/env(x),c(0,1),maximum=T)$objective

 

f2 = function(n)

  {

    rand = vector("numeric",0)

B)

nfactor =function(n)

  {

    u = runif(n,0,1)

    theta=mean(f(u))

x=f1(u)

    theta=mean((f(x)/env(x)*a))

    cat("normalising factor   ",theta,"\n")

    f(x)*a/env(x)