## function to generate the data and plot it
generate.and.plot = function(a,b,M,N,x1){
x=numeric(N)
x[1] = x1
for(i in 2:N){
x[i]= (a*x[i-1]+b)%%(M+1)
}
u=x/(M+1)
## plot it
op=par(mfrow=c(2,2), pty='s', mar=rep(1,4), pch=20, cex=.4)
# plot pairs (u_i, u_{i+1}) generated by the congruential generator
plot(u[-N], u[-1])
u2=runif(N)
# plot pairs (u_i, u_{i+1}) of a "true" uniform
plot(u2[-N],u2[-1], col=1)
#acf(u)
#pacf(u)
#par(op)
}

N=1e4

## very bad choices of (a,b,M)
M=1e4-1
a=2
b=0
x1=1
generate.and.plot(a,b,M,N,x1)

##  bad choices of (a,b,M)
M=1e4-1
a=336
b=1
x1=5041
generate.and.plot(a,b,M,N,x1)

## not too bad choices of (a,b,M)
M=1e4-1
a=69069
b=1
x1=41
generate.and.plot(a,b,M,N,x1)


## better choices of (a,b)
M=15485863-1
a=69069
b=0
x1=41
generate.and.plot(a,b,M,N,x1)


## better choices of (a,b)
M=2^31-2
a=7^5
b=0
x1=41
generate.and.plot(a,b,M,N,x1)







##############################
# Accept-reject algorithm of a Beta(2.7,6.3) based on uniforms
N=1000
alpha=2.7
beta=6.3
ys=runif(N)
M=dbeta((alpha-1)/(alpha + beta -2),alpha,beta)
us=M*runif(N)
val=seq(0,1.,.01)
valf=dbeta(val,shape1=alpha,shape2=beta)
plot(val,valf,type="l",xlab="x",ylab="density",lwd=2)
polygon(c(val,rev(val)),c(valf,0*valf),col="gold2")
points(ys,us,cex=.4,pch=20)
sum(us<=dbeta(ys, alpha, beta))/N


##############################
## Accept-reject method based on Gaussian   
# function 
dev.off()
N=1000
f=function(x){
exp(-x^2/2)*(1+sin(6*x)^2+3*cos(x)^2*sin(4*x)^2)}
# ratio
rat=function(x){
.2* (1+sin(6*x)^2+3*cos(x)^2*sin(4*x)^2)}
#
# Accept-Reject
refn=rnorm(N)
refu=runif(N)
ind.accept=refu[]<rat(refn[])
## accepted points
x.accept=refn[ind.accept]
y.accept=refu[ind.accept]
## refused points
x.refused=refn[!ind.accept]
y.refused=refu[!ind.accept]
#
polb=f(seq(-4,4,.01))
polo=5*exp(-seq(-4,4,.01)^2/2)
plot(seq(-4,4,.01),polo,type="l",xlab="",ylab="", main=paste("Acceptance rate=", sum(ind.accept)/N))
lines(seq(-4,4,.01),polb)
polygon(c(seq(-4,4,.01),seq(4,-4,-.01)),c(polb,0*seq(4,-4,-.01)),col="sienna3")
polygon(c(seq(-4,4,.01),seq(4,-4,-.01)),c(polo,f(seq(4,-4,-.01))),col="wheat")
points(x.accept,y.accept*f(x.accept)/rat(x.accept),pch=20,cex=.3)
points(x.refused,y.refused*f(x.refused)/rat(x.refused),pch=18,cex=.9)





##############################
## Accept-reject method to generate a Gaussian, based on a double exponential proposal distribution 
# function 
N=2000
f=dnorm
g=function(x){
    .5*exp(-abs(x))
}
M=2
# ratio
rat=function(x){
    f(x)/(M*g(x))
}
#
# Accept-Reject
refn=(2*rbinom(N, 1, .5) - 1)*rexp(N,1)
refu=runif(N)
ind.accept=refu[]<rat(refn[])
## accepted points
x.accept=refn[ind.accept]
y.accept=refu[ind.accept]
## refused points
x.refused=refn[!ind.accept]
y.refused=refu[!ind.accept]
#
polb=f(seq(-4,4,.01))
polo=M*g(seq(-4,4,.01))
plot(seq(-4,4,.01),polo,type="l",xlab="",ylab="", main=paste("Acceptance rate=", sum(ind.accept)/N))
lines(seq(-4,4,.01),polb)
polygon(c(seq(-4,4,.01),seq(4,-4,-.01)),c(polb,0*seq(4,-4,-.01)),col="sienna3")
polygon(c(seq(-4,4,.01),seq(4,-4,-.01)),c(polo,f(seq(4,-4,-.01))),col="wheat")
points(x.accept,y.accept*f(x.accept)/rat(x.accept),pch=20,cex=.3)
points(x.refused,y.refused*f(x.refused)/rat(x.refused),pch=18,cex=.9)
sum(ind.accept)/N

qqnorm(x.accept)