#Binary Markov chain rm(list=ls()) #Legth of chain n<-100 #Initial value x0<-0 #Vector to store values x<-rep(NA,n) #Set first value if (x0==0){x[1]<-rbinom(1,1,0.2)} if (x0==1){x[1]<-rbinom(1,1,0.7)} #We use rbinom(1,1,p) to simulate #p is the probability X=1 #Put this in a loop for (i in 2:n){ if (x[i-1]==0){x[i]<-rbinom(1,1,0.2)} if (x[i-1]==1){x[i]<-rbinom(1,1,0.7)} } plot(1:n,x,"l",main='Sequence of X',xlab = 'i') #Transition Matrix P<-matrix(c(0.8,0.2,0.3,0.7),2,2,byrow = TRUE) #P<-matrix(c(0.9,0.1,0.15,0.85),2,2,byrow = TRUE) #Set pi_1 pi<-c(0.8,0.2) #pi<-c(0.15,0.85) #Set n n<-100 #Repeat Recursion for (i in 1:(n-1)){ pi<-pi%*%P } print(pi) #Code to keep pi and plot it #Initialise all values of pi in a matrix pi<-matrix(NA,n,2) #Set first value of pi pi[1,]<-c(0.8,0.2) #Recursion for (i in 2:n){ pi[i,]<-pi[i-1,]%*%P } #Two Plots par(mfrow=c(2,1)) #First plot plot(1:n,pi[,1],'l',main=expression(X[0]==0),xlab = 'i',ylab=expression(Pr(X[i]==0))) #Repeat with different initial condition pi[1,]<-c(0.3,0.7) for (i in 2:n){ pi[i,]<-pi[i-1,]%*%P } plot(1:n,pi[,1],'l',main=expression(X[0]==1),xlab = 'i',ylab=expression(Pr(X[i]==0))) #Random Walk # set.seed(123) x0<-0 n<-100 x<-rep(NA,n) x[1]<-rnorm(1,x0,1) for (i in 2:n){ x[i]<-rnorm(1,x[i-1],1) } plot(1:n,x,"l") #Metropolis Algorithm (Simple) set.seed(3) x0<-0 n<-110000 x<-rep(NA,n) #Metropolis Algorithm First Step xold<-x0 #Set x old xnew<-rnorm(1,xold,100) #Generate proposal a<-min(1,dt(xnew,5)/dt(xold,5)) #Compute acceptance probability accept<-rbinom(1,1,a) #Accept yes/no if(accept){x[1]<-xnew} #If accept x1=xnew if(!accept){x[1]<-xold} #Otherwise accept x1=xold sd = 100 #Now inside loop for (i in 2:n){ xold<-x[i-1] xnew<-rnorm(1,xold,sd) a<-min(1,dt(xnew,5)/dt(xold,5)) accept<-rbinom(1,1,a) if(accept){x[i]<-xnew} if(!accept){x[i]<-xold} } #pdf('../Pics/sp5.pdf') plot(1:1000,x[1:1000],"l",ylab = 'x',xlab='',main=paste('Initial Value =0, Proposal S.D. =', sd)) #dev.off() #Random Walk Metropolis (More complicated) x0<-0 bi<-10000 #Burn In n<-110000 #Sample Size accrate<-0 #Counter for acceptance rate x<-rep(NA,n) #Initialise vector to store x nu<-2 #Proposal Variance for (i in 1:n){ xold<-ifelse((i==1),x0,x[i-1]) xnew<-rnorm(1,xold,nu) a<-min(1,dt(xnew,5)/dt(xold,5)) accept<-rbinom(1,1,a) if(accept){x[i]<-xnew} if(accept&&(i>bi)){accrate<-accrate+1} if(!accept){x[i]<-xold} } accrate<-accrate/(n-bi) plot(1:n,x,"l") library("coda") geweke.diag(x[10001:110000]) effectiveSize(x[10001:110000])