markov=function(x,theta,lambda,k,run)
{
n=length(x)
Q=0
output=matrix(0,ncol=3,nrow=run)
for(i in 1:run)
{
if(k>1)
theta=rbeta(1,1+sum(x[1:k]),1+k-sum(x[1:k]))
if(k<n)
lambda=rbeta(1,1+sum(x[(k+1):n]),1+n-k+sum(x[(k+1):n]))
for(j in 1:n)
{
Q[j]=((1/k)/(n-k)/(n-1)((1-lambda)*theta/lambda/(1-theta))**sum(x[1:j]))*((1-theta)/(1-lambda))**sum(x[1:j])
}
k=sample(n,size=1,prob=Q)
output[i,]=c(theta,lambda,k)
}
output
}
z=markov(x,0.5,0.5,1,200)
