PCM_Metropolis_MacKay_Ex29

# Markov Chain (Exercise 29.3 in MacKay, 2003, p.367)
#    Computation of the stationary state distribution
#    R-Code by PCM 2015/05/18

ssqf <- function(v1, v2){ssq <- sum((v1 - v2)^2)} # sum-of-squares of ssqf = (v1-v2)'(v1-v2)

my_var <- function(x, p, xmean)
  # x = vector of state numbers
  # p = vector of probabilities p(x)
  # xmean = sum(x)/n
{
  return(sum(p * (x - xmean)^2))
}

# this function's code borrowed from Matloffs "The Art of R Programming" (2011, p.200)
find_pi1 <- function(p)  
  # p = transition probability matrix of Markov chain
{
  n <- nrow(p)
  imp <- diag(n) - t(p)
  imp[n,] <- rep(1,n)
  rhs <- c(rep(0,n-1),1)
  pivec <- solve(imp, rhs)
  return(pivec)
}

mc_time_series <- function(nt=800, x0=10, eps=1.0E-9)
  # stationary_pmf = stationary probability mass function 
  # nt  = maximum of time steps (600 by default)
  # x0  = start state (10 by default)
  # eps = stop criterion for iteration process (10.0E-9 by default)
{ cat("\n'Drunken' random walk \n    R-code by PCM 2015/05/18\n\n")
  p_target <- vector(mode="numeric",length=21)
  p_target <- find_pi1(T)                           # target distribution = stationary distribution
  print(p_target)
  p <- matrix(0.0, nrow=nt, ncol=21, byrow=TRUE)    # initialization of distribution vector
  t <- 1 # initialization of time index
  states <- c(0:20)
  p[t,] <- c(.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,.0)
  p[t, x0+1] <- 1.0
  xmean <- sum(states*p[t,])
  variance <- my_var(states,p[t,],xmean)
  ssq_old <- 1.0
  ssq_new <- ssqf(p_target, p[t,])
  barplot(p[t,], names.arg=c(0:20), col="steelblue", 
          main=sprintf("state distribution at step t = %i", t),
          xlab="states", ylab="probability p", ylim=c(0,1))
  legend(x=0,y=1.0,legend=sprintf("t=%i,mean = %f, variance = %f", t, xmean, variance))
  while (abs(ssq_old - ssq_new) > eps)
  { tp1 <- t+1
    p[tp1,] <- p[t,] %*% T
    xmean <- sum(states*p[tp1,])
    variance <- my_var(states,p[tp1,],xmean)
    ssq_old <- ssq_new
    ssq_new <- ssqf(p[tp1,],p[t,])
    if((tp1 <= 10) || (tp1 %% 100 == 0) || (tp1 == 121) || (tp1 == 441))
    {barplot(p[tp1,], names.arg=c(0:20), col="steelblue", 
             main=sprintf("state distribution at step t = %i", tp1),
             xlab="states", ylab="probability p", ylim=c(0,1))
     legend(x=0,y=1.0,legend=sprintf("t=%i,mean = %f, variance = %f", tp1, xmean, variance))
     p20 <- p[t,21]
     print(sprintf("p[%i,20] = %f", tp1, p20))} # p of being in state 20
    t <- tp1
  }
  xmean < sum(states*p[t,])
  variance <- my_var(states,p[t,],xmean)
  barplot(p[t,], names.arg=c(0:20), col="steelblue", 
          main=sprintf("state distribution at step t = %i", t),
          xlab="states", ylab="probability p", ylim=c(0,1))
  legend(x=0,y=1.0,legend=sprintf("t=%i,mean = %f, variance = %f", t, xmean, variance))
  cat("\nssq = ", ssq_new, "\n\n")
  p20 <- p[t,21]     # p of being in state 20
  print(sprintf("p[%i,20] = %f", tp1, p20))
}

# T = transition matrix P(State=j at t+1 | State=i at t) = t[i, j]
T <- matrix(c(
  .50,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 01
  .50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 02
  .00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 03
  .00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 04
  .00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 05
  .00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 06
  .00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 07
  .00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 08
  .00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 09
  .00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 10
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,  # row 11
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,  # row 12
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,  # row 13
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,  # row 14
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,  # row 15
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,  # row 16
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,  # row 17
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,  # row 18
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,  # row 19
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,  # row 20
  .00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.50)  # row 21
  , nrow=21, ncol=21, byrow=TRUE)
#

mc_time_series()

## 
## 'Drunken' random walk 
##     R-code by PCM 2015/05/18
## 
##  [1] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
##  [7] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
## [13] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
## [19] 0.04761905 0.04761905 0.04761905

## [1] "p[2,20] = 0.000000"

## [1] "p[3,20] = 0.000000"

## [1] "p[4,20] = 0.000000"

## [1] "p[5,20] = 0.000000"

## [1] "p[6,20] = 0.000000"

## [1] "p[7,20] = 0.000000"

## [1] "p[8,20] = 0.000000"

## [1] "p[9,20] = 0.000000"

## [1] "p[10,20] = 0.000000"

## [1] "p[100,20] = 0.048890"

## [1] "p[121,20] = 0.045327"

## [1] "p[200,20] = 0.048377"

## [1] "p[300,20] = 0.047869"

## [1] "p[400,20] = 0.047700"

## [1] "p[441,20] = 0.047568"

## [1] "p[500,20] = 0.047646"

## [1] "p[600,20] = 0.047628"

## [1] "p[700,20] = 0.047622"

## 
## ssq =  4.359994e-08 
## 
## [1] "p[713,20] = 0.047621"

mc_time_series(x0=0)

## 
## 'Drunken' random walk 
##     R-code by PCM 2015/05/18
## 
##  [1] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
##  [7] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
## [13] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
## [19] 0.04761905 0.04761905 0.04761905

## [1] "p[2,20] = 0.000000"

## [1] "p[3,20] = 0.000000"

## [1] "p[4,20] = 0.000000"

## [1] "p[5,20] = 0.000000"

## [1] "p[6,20] = 0.000000"

## [1] "p[7,20] = 0.000000"

## [1] "p[8,20] = 0.000000"

## [1] "p[9,20] = 0.000000"

## [1] "p[10,20] = 0.000000"

## [1] "p[100,20] = 0.017350"

## [1] "p[121,20] = 0.023023"

## [1] "p[200,20] = 0.037442"

## [1] "p[300,20] = 0.044306"

## [1] "p[400,20] = 0.046541"

## [1] "p[441,20] = 0.046931"

## 
## ssq =  4.390829e-08 
## 
## [1] "p[482,20] = 0.047190"

PCM_Metropolis_MacKay_Ex29_3.R

PCM

Mon Feb 22 11:35:17 2016