MCMC in R

Random Variable Generation

Monte Carlo methods rely on the possibility of producing a supposedly endless flow of random variables, and for well-known or new distributions. Such a simulation is, in turn, based on the production of uniform random variables on the interval (0, 1)

# check the generator
Nsim = 10^4
x=runif(Nsim)
x1=x[-Nsim]
x2=x[-1]
par(mfrow=c(1,3))
hist(x)
plot(x1,x2)
acf(x)

# Inverse Transform
Probability Integral Transform is related to the result that data values that are modelled as being random variables from any given continuous distribution can be converted to random variables having a standard uniform distribution.

# generate exponential
Nsim = 10^4
U=runif(Nsim)
X=-log(U)
Y=rexp(Nsim)
par(mfrow=c(1,2))
hist(X,freq=F)
hist(Y,freq=F)

# generate Chisq(6)
U=runif(3*10^4)
U=matrix(data=U, nrow=3)
X=-log(U)
X=2*apply(X,2,sum)
Y=rchisq(Nsim,6)
par(mfrow=c(1,2))
hist(X,freq=F)
hist(Y,freq=F)

# generate normal by Box-Muller
U=runif(2*10^4)
U=matrix(data=U, ncol=2)
X1=sqrt(-2*log(U[,1]))*cos(2*pi*U[,2])
X2=sqrt(-2*log(U[,1]))*sin(2*pi*U[,2])
Y=rnorm(10^4)
par(mfrow=c(1,3))
hist(X1,freq=F)
hist(X2,freq=F)
hist(Y)

# Poisson for large lambda
Nsim=10^4
lambda=100
spread=3*sqrt(lambda)
t=round(seq(max(0,lambda-spread),lambda+spread,1))
prob=ppois(t,lambda)
X=rep(0,Nsim)
for (i in 1:Nsim){
u=runif(1)
X[i]=t[1]+sum(prob<u)-1
}
hist(X,freq=F)

# Accept-Reject Algorithm to generate beta with Uniform
1 We generate a candidate random variable
2 Only accept it subject to passing a test

Nsim=2500
a=2.7
b=6.3
M=2.67
u=runif(Nsim, max=M)
y=runif(Nsim)
x=y[u<dbeta(y,a,b)]
hist(X,freq=F,breaks=100)

U = runif(100000,0,1)
accept = c()
for(i in 1:length(U)){
U1 = runif(1, 0, 1)
if(dunif(U[i], 0, 1)*3*U1 <= dbeta(U[i], 6, 3)) {
accept[i] = 'Yes'
}
else if(dunif(U[i],0,1)*3*U1 > dbeta(U[i], 6, 3)) {
accept[i] = 'No'
}
}
T = data.frame(U, accept = factor(accept, levels= c('Yes','No')))
hist(T[,1][T$accept=='Yes'], breaks = seq(0,1,0.01), freq = FALSE, main = 'Histogram of X', xlab = 'X')

 Monte Carlo Integration

The validity of Monte Carlo approximations relies on the Law of Large Numbers
The versatility of the representation of an integral as an expectation

We apply probabilistic results of Law of Large Numbers and Central Limit Theorem,