[R] Stochastic (transition) matrices: how to determine distributions and variance?

[Chris Stubben]

Jonathan Greenberg-2 wrote:
> 
> I am trying to determine the distribution and variance for a classic 
> stochastic (transition) matrix 
> 

I have a few suggestions that I've tried with matrix population models
(so probabilities are survival rates plus dead fates).

# here's a test example

A<-matrix(c(.1,.3,.6,.2,.4,.4, 0,0,1), nrow=3)
n<-c(100, 300, 200)

A %*% n
     [,1]
[1,]   70
[2,]  150
[3,]  380
 
# columns should add to 1, correct?
colSums(A)


#  two possible solutions

# option 1.   sample A from a multinomial distribution 


apply(A, 2, function(x) rmultinom(1, size = 1000, prob=x) /1000)
      [,1]  [,2] [,3]
[1,] 0.090 0.195    0
[2,] 0.322 0.401    0
[3,] 0.588 0.404    1


n1<-matrix(numeric(1000*3), ncol=3)
for(i in 1:1000)
{
A2<-apply(A, 2, function(x) rmultinom(1, size = 1000, prob=x) /1000)
n1[i,] <-A2  %*% n
}

apply(n1, 2, quantile, c(0.025,0.975))
      [,1]  [,2]     [,3]
2.5%  62.3 140.5 370.4975
97.5% 78.3 159.6 389.8000


# option 2.  Since I usually don't know anything about the distribution of
values in A, I
would resample transitions using a bootstrap


# Get number of each transition

a2<-data.frame( round ( as.table(t(A)*n)))
  Var1 Var2 Freq
1    A    A   10
2    B    A   60
3    C    A    0
4    A    B   30
5    B    B  120
6    C    B    0
7    A    C   60
8    B    C  120
9    C    C  200

## expand into 600 rows (so A->A repeats 10 times, B->A repeats 60 times,
etc)

r <- rep(row.names(a2), a2$Freq)
y <- a2[r, 1:2 ]
rownames(y)<-1:nrow(y)

##  this gets the original matrix
prop.table( table(y[,2], y[,1]), 2)
   
# so resample rows in y 

z <- sample(nrow(y), replace = TRUE)
y1<-y[z,]
prop.table( table(y1[,2], y1[,1]), 2)  %*% n

## and repeat 1000 times

n2<-matrix(numeric(1000*3), ncol=3)

for(i in 1:1000)
{
z <- sample(nrow(y), replace = TRUE)
y1<-y[z,]
n2[i,] <-prop.table( table(y1[,2], y1[,1]), 2)  %*% n
}

# the CIs are much wider here

apply(n2, 2, quantile, c(0.025,0.975))
          [,1]     [,2]     [,3]
2.5%  55.89959 132.6868 361.2091
97.5% 85.11097 169.3798 399.7032


Chris Stubben
-- 
View this message in context: http://www.nabble.com/Stochastic-%28transition%29-matrices%3A-how-to-determine-distributions-and-variance--tp25209297p25227303.html
Sent from the R help mailing list archive at Nabble.com.