# Stat 241/541 fall 2005

# variations on the HHH vs. TTHH example

# setwd("/Users/pollard/Academic/Courses/241/241.fall2005/notes2005")

# create transition matrix for HHH vs TTHH game with coin that lands
  heads with prob  p

transit.H3T2H2 = function(p = 1/2){  # default value is fair coin
               states
               =c("S","H","T","HH","TT","TTH","M.wins","R.wins")
               M = matrix(
               c(0,p,1-p,0,0,0,0,0,
               0,0,1-p,p,0,0,0,0,
               0,p,0,0,1-p,0,0,0,
               0,0,1-p,0,0,0,p,0,
               0,0,0,0,1-p,p,0,0,
               0,0,1-p,0,0,0,0,p,
               0,0,0,0,0,0,1,0,
               0,0,0,0,0,0,0,1
               ), nrow=8,byrow=T,dimnames=list(states,states))
               
               return(M)
}

H3.T2H2 = transit.H3T2H2() # transition matrix for original HHH vs
TTHH game



# find probs that M wins starting from various states 
Mwin = function(tm=H3.T2H2){  # tm is the transition matrix

idM = c(rep(0,6),1,0)
idR = c(rep(0,7),1)

M = diag(8)-tm
M[7,] = idM
M[8,] = idR
b = c(rep(0,6),1,0)
winprob = solve(M,b)
return(list(trans=tm,eqn=cbind(M,rhs=b), probs=winprob))
}



# transition matrix for TTHH
transit.T2H2 = function(p=1/2){
             states =c("S","T","TT","TTH","TTHH")
             M = matrix(
             c(p,1-p,0,0,0,
             p,0,1-p,0,0,
             0,0,1-p,p,0,
             0,1-p,0,0,p,
             0,0,0,0,1
             ), nrow=5, byrow=T, dimnames=list(states,states))
             
             return(M)
}

T2H2 = transit.T2H2()

# transition matrix for HHH
transit.H3 = function(p=1/2){
           states =c("S","H","HH","HHH")
           M = matrix(
           c(1-p,p,0,0,
           1-p,0,p,0,
           1-p,0,0,p,
           0,0,0,1
           ), nrow=4, byrow=T, dimnames=list(states,states))
           
           return(M)
}

H3 = transit.H3()


# find expected time to absorption, starting from various states 
Etime = function(tm=H3.T2H2){  # tm is the transition matrix
      absorbing =  diag(tm) == 1  # find the absorbing states
      b = 1 - absorbing
      nn = dim(tm)[1]  # number of states
      M = diag(nn)- tm
      M[absorbing,] = tm[absorbing,]
      etime = solve(M,b)
      return(list(trans=tm,eqn=cbind(M,rhs=b), expected=etime))
}


# quick way of raising a square matrix to a high power
mpower = function(mm,k){#raise matrix mm to power 2^k
       for (i in 1:k){
           mm = mm %*% mm
           }
           return(mm)
}


# Things to try:

# round(mpower(H3.T2H2,10),3)  
#   Meigen = eigen(H3.T2H2)

# How can we find probs of absorption from Meigen?


# Compare output from
# 
# Etime(H3)
# Etime(T2H2)
# Etime(H3.T2H2)