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

全文链接：http://tecdat.cn/?p=30829

1. Consider the following density:

image.png

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

2

(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).

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

image.png

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

image.png

image.png

image.png

normalfactor =function(n)

  {

image.png

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)))

image.png

输出前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

image.png

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

计算envelop function的累计密度函数

image.png

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

 

f2 = function(n)

  {

    rand = vector("numeric",0)

image.png

B)

nfactor =function(n)

  {

    u = runif(n,0,1)

    theta=mean(f(u))

image.png

x=f1(u)

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

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

    f(x)*a/env(x)

image.png

image.png

image.png

image.png

image.png