################################################

# now a function that generates from an arbitrary
# discrete, finite distribution - use the Splus function sample()

# for example, 10 draws from the distribution over the values
# {1,2,3,4} given by p={.1,.5,.3,.1} can be obtained with the
# commands below

values_c(1,2,3,4)
p_c(.1,.5,.3,.1)

x_sample(values,size=10,replace=T,prob=p)


#################################################

# now build a function that samples from
# a MC with TPM P, and starting value x0

rmc_function(n,x0,P){
  # n draws from a MC with initial value x0 and tpm P
  # x0 must be an integer between 1 and length(x0)
 p_dim(P)[1]
 x_rep(NA,n)
 xcur_x0
 for(i in 1:n){
   x[i]_sample(1:p,size=1,prob=P[xcur,])
   xcur_x[i]
 }
 return(x)
}

# now for a 2 state example

tpm_matrix(scan(),ncol=2,nrow=2,byrow=T)
          .9 .1
          .1 .9

x0_1
par(mfcol=c(2,3),mar=c(5,5,2,2))
x_rmc(50,x0,tpm)
plot(x,type="b")
hist(x); mtext("n = 50",side=3,line=.5)
x_rmc(100,x0,tpm)
plot(x,type="b")
hist(x); mtext("n = 100",side=3,line=.5)
x_rmc(500,x0,tpm)
plot(x,type="b")
hist(x); mtext("n = 500",side=3,line=.5)


# now for a 4 state example

tpm_matrix(scan(),ncol=4,nrow=4,byrow=T)
          .2 .2 .2 .4
          .1 .1 .1 .7
          .6 .2 .2  0
          .5 .5 .0 .0

matpower_function(P,pow){
  # a function that computes P^pow
  if(pow==1) return(P)
  else return( P %*% matpower(P,pow-1) )
}

# the limit dist via powering up tpm
> matpower(tpm,20)
          [,1]      [,2]      [,3]      [,4] 
[1,] 0.3111366 0.2666870 0.1111028 0.3110737
[2,] 0.3111782 0.2667203 0.1110891 0.3110124
[3,] 0.3110774 0.2666398 0.1111222 0.3111606
[4,] 0.3110401 0.2666100 0.1111344 0.3112155

# the estimated limit from a sample of 10,000
x_rmc(10000,x0,tpm)
> table(x)/10000
      1      2      3      4 
 0.3083 0.2693 0.1109 0.3115