R語言代做編程輔導(dǎo)和解答M3S9/M4S9 Stochastic Simulation: Project 2
全文鏈接：http://tecdat.cn/?p=30829
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).
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個(gè)分布的值
?#envelop function envenv =function(x)
? {
??? if(x<=0)v=0??? else if(x<=0.01)v=330*x
??? else if(x<=0.03)v=33
黑色代表函數(shù)值，綠色代表envelop function的擬合值。
?
計(jì)算envelop function的累計(jì)密度函數(shù)
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)
