Summary of Markov Chain Demo Script

EXAMPLE 1. BIRTH-DEATH MODEL -- for queueing with M servers

X(n+1) = X(n) + 1 wp  lam / (lam + min(X(n),M)*mu) = p(X(n),X(n)+1)
              - 1 wp  1 - "                        = p(X(n),X(n)-1)

Explain this by saying: if next event occurring is at minimum of
min(x,M) indep Expon(mu) service times or an Expon(lam) arrival time
then Prob that the arrival comes first is lam/(lam + min(x,M)*mu)

The only way to achieve balance is for     lam < M*mu and

             P(X(T)=k) * p(k,k+1) = 

pi(k+1)*min(k+1,M)*mu = pi(k)*lam


## R steps -- simulation

> Xnew = function(Xold, U, lam, mu, M) 
        #### M-server queueing discrete-time birth-death
         Xold + (2*(U < lam/(lam+min(Xold,M)*mu))-1)

> Xvals = numeric(1e5+1)    ## initializes 1e5+1 - vector to 0
  Uvec = runif(1e5)

> for (i in 1:1e5) Xvals[i+1] = Xnew(Xvals[i], Uvec[i], 3, 1, 5)


> table(Xvals[-1]) 

    0     1     2     3     4     5     6     7     8     9    10    11    12
 2367  9554 17983 21377 18344 12380  7384  4397  2589  1508   887   517   276
   13    14    15    16    17    18    19    20
  157    96    62    52    39    18     9     4
 
tmp = table(Xvals[-1])
round(tmp[1:17]/tmp[2:18],3)

    0     1     2     3     4     5     6     7     8     9    10    11    12
0.248 0.531 0.841 1.165 1.482 1.677 1.679 1.698 1.717 1.700 1.716 1.873 1.758
   13    14    15    16
1.635 1.548 1.192 1.333

### Compare with "theoretical" ratios

> round(pmin(1:17,5)*(3+pmin(0:16,5))/(3*(3+pmin(1:17,5))),3)
 [1] 0.250 0.533 0.833 1.143 1.458 1.667 1.667 1.667 1.667 1.667 1.667 1.667
[13] 1.667 1.667 1.667 1.667 1.667
                                         ### VERY CLOSE AGREEMENT !
> plot(Xvals[2:1001])

> vec1 = table(Xvals[2:50001])/50000       ## values 0:18
  vec2 = table(Xvals[50002:1e5+1])/50000   ## values 0:20

> c(vec1)
      0       1       2       3       4       5       6       7       8       9
0.02348 0.09504 0.17824 0.21370 0.18738 0.12728 0.07436 0.04330 0.02516 0.01434
     10      11      12      13      14      15      16      17      18
0.00774 0.00446 0.00258 0.00128 0.00074 0.00046 0.00030 0.00014 0.00002

> c(vec2)
      0       1       2       3       4       5       6       7       8       9
0.02386 0.09604 0.18142 0.21382 0.17950 0.12032 0.07332 0.04464 0.02662 0.01582
     10      11      12      13      14      15      16      17      18      19
0.01000 0.00588 0.00294 0.00186 0.00118 0.00078 0.00074 0.00064 0.00034 0.00018
     20
0.00008

### Here are two ways to check that the process has "settled down" to equilibrium

> summary(c(vec1)-c(vec2)[1:19])
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-3.180e-03 -1.380e-03 -4.400e-04  1.474e-05 -3.200e-04  7.880e-03


> plot(cumsum(Xvals[-1] <=4)/(1:1e5), xlab="Time-steps", ylab="Fraction <=4",
    main="Cumulative fraction of steps with X(t)<=4", ylim=c(.6,1))



EXAMPLE 2. Lifetime Model Age increases successively +1 until 
year T (age at death) when it reverts instantaneously to 0.

"Stationary Population Model"

### Suppose that the successive probabilities of T >= t for integers t
###    are given by   exp(-.01*t^1.3) , t=1,...,105 (max age 105)


ages = numeric(1+1.1e6)
Uvec = runif(3e4)

LifNxt = function(U,scal, powr) 
          trunc((-log(U)/scal)^(1/powr))
      
### This works for S(k) = P(T>=k) = P(U <= S(k))
###   where we calculate the largest k such that U <= S(k)

ctr=1
for(i in 1:3e4) { 
     k = LifNxt(Uvec[i],.01,1.3)
     ages[ctr+(1:k)] = if(k>1) c(1:(k-1),0)  else 0
     ctr=ctr+k }
> ctr
[1] 946164

> mean(ages[1:(ctr-1)])
[1] 25.14916
> mean(ages[1:1e4])
[1] 25.0439
> mean(ages[(1e4+1):2e4])
[1] 23.2936
> mean(ages[1:1e4]^2)
[1] 1101.07
> mean(ages[(1e4+1):2e4]^2)
[1] 947.942
> mean(ages[1:(ctr-1)]^2)
[1] 1155.251

### These calculations are intended to show that the ages entries
are fluctuating in a very stable way over long times.