## y|theta ~ B(100,theta)
## theta ~ Beta(2,3)
## target marginale de y (beta-binomial distribution)
Niter <- 10000
theta <- y <- rep(0,Niter)
theta[1] <- 0.5
n <- 100
for (i in 2:Niter)
{
y[i] <- rbinom(1,n,theta[i-1])
theta[i] <- rbeta(1,y[i]+2,n-y[i]+3)
}
plot(theta,type="l",lwd=2)
plot(acf(theta))
hist(theta,prob=TRUE,nclass=30,col="blue")
curve(dbeta(x,2,3),add=TRUE,lwd=2,col="green")
proby <- 12*choose(n,0:n)*gamma((0:n)+2)*gamma(n-(0:n)+3)/
gamma(n+5)
barplot(table(as.factor(y))/Niter,add=TRUE,col="red")
barplot(proby,add=TRUE,col="green")
## y|theta ~ B(100,theta)
## theta|alpha ~ Beta(alpha,1)
## alpha ~ U(0,10)
## y is observed and equal to 30
## target theta|y
Niter <- 10000
theta <- alpha <- rep(0,Niter)
n <- 100
y <- 30
alpha[1] <- 5
dcalpha <- function(alpha,theta)
{
alpha*theta^alpha
}
for (i in 2:Niter)
{
theta[i] <- rbeta(1,y+alpha[i-1],n-y+1)
alphatilde <- rnorm(1,alpha[i-1],1)
if (alphatilde<0 || alphatilde>10) rho <- 0
else rho <-dcalpha(alphatilde,theta[i])/
dcalpha(alpha[i-1],theta[i])
if (runif(1)<=rho) alpha[i] <- alphatilde
else alpha[i] <- alpha[i-1]
}
plot(theta,type="l",lwd=2)
plot(acf(theta))
hist(theta,prob=TRUE,nclass=30,col="blue")