#Maple Code for Dr. Z.'s Dynamcial Methods in Biology, Lecture 7

Help7:=proc(): print(` GR(p,i,N), GRt(p,i,N),  GRm(N,p), OneStepMarkov(P,i), MarkovTrip(P,K), StSa(P,K) ,  StS(P), StSp(P,K), RandSM(N)  `): end:

with(Statistics):
with(LinearAlgebra):

#GR(p,i,N): Simulating a Gambler's Ruin problem. Ypu enter a casino with i dollars. At every round you win a dollar with prob. p and lose a
#dollar with prob. 1-p. You exit as soon as you are broke or got the maximum allowed  that is N. It  returbs [0,c] or [1,c]
#if you lost or won respectively, where c is the number of rounds. Try: GR(3/5,5,10);
GR:=proc(p,i,N) local X,x,d,c:
X:=RandomVariable(Bernoulli(p)):
c:=0:
x:=i:
 print(`I enter the casino with `, i, `dollars `):
while x>0 and x<N do
 c:=c+1:
 d:=trunc(Sample(X,1)[1]):
  if d=0 then
    x:=x-1:
   print(`I lost a dollar, now I have`, x, `dollars `):
   else
  x:=x+1:
   print(`I lost a dollar, now I have`, x, `dollars `):
  fi:
 od:
 if x=0 then
  print(`I am so sad, I am broke`):
  print(`It took`, c, `rounds `):
   RETURN([0,c]):
  elif x=N then
   print(`Yea, I am rich, I now have the max`, N):
    print(`It took`, c, `rounds `):
   RETURN([1,c]):
  fi:
end:


#GRt(p,i,N): A terse version of GR(p,i,N). Returns [0,c] or [1,c], where c is the number of rounds
GRt:=proc(p,i,N) local X,x,d,c:
X:=RandomVariable(Bernoulli(p)):
c:=0:
x:=i:
while x>0 and x<N do
 c:=c+1:
 d:=trunc(Sample(X,1)[1]):
  if d=0 then
    x:=x-1:
   else
  x:=x+1:
  fi:
 od:
 if x=0 then
   RETURN([0,c]):
  elif x=N then
   RETURN([1,c]):
  fi:
end:


#GRm(N,p): The (N+1)x(N+1) Transition matrix for a Gambler's ruin with max. capital N
GRm:=proc(N,p) local i:
[
[1,0$N],
seq([0$(i-1),(1-p),0,p,0$(N-1-i)],i=1..N-1),
[0$N,1]
]:
end:


##OneStepMarkov(P,i): Suppose that there are n sites, let P be an n by n  transition matrix given as a list of lentgh n 
#where P[i][j] is the probability of getting from site i to site j. Outputs the state visited next

OneStepMarkov:=proc(P,i)
trunc(Sample(RandomVariable(ProbabilityTable(P[i])),1)[1]):
end:

#MarkovTrip(P,K): The trajectory of a Markov chain with transition matrix P with K steps
#Try: MakovTrip([[0.2,0.8],[0.4,0.6]],1000);
MarkovTrip:=proc(P,K) local x,L,i:

x:=1:
L:=[1]:

for i from 1 to K do
x:=OneStepMarkov(P,x):
L:=[op(L),x]:
od:

L:

end:

#StSa(P,K): Given a transition matrix, ESTIMATES the steady-state probability distribution of the Markov chain using  a walk with K steps
StSa:=proc(P,K) local L,x,f,i:
L:=MarkovTrip(P,K):
f:=add(x[L[i]],i=1..K):
evalf([seq(coeff(f,x[i],1)/K,i=1..nops(P))]):

end:





#StS(P,K): Given a transition matrix, FIND  the EXACT steady-state probability distribution of the Markov chain using  linear algebra
StS:=proc(P) local x,n,eq,var,i,j:
n:=nops(P):
var:={seq(x[i],i=1..n)}:
eq:={seq( add(x[i]*P[i][j],i=1..n)=x[j],j=1..n), add(x[i],i=1..n)=1}:
subs(solve(eq,var),[seq(x[i],i=1..n)]):
end:


#StSp(P,K): Given a transition matrix, FIND  the approximate steady-state probability distribution  using raising the matrix to  the K-th power
StSp:=proc(P,K) local Q,i:
Q:=evalm(P&^K):
[seq(Q[i,1],i=1..nops(P))]:
end:


#RandPV(N): A radnom stochastc  vector of length N
RandPV:=proc(N) local ra,i,v,s:
ra:=rand(1..1000):
v:=[seq(ra(),i=1..N)]:
s:=add(v[i],i=1..N):
v/s:
end:

#Added after class
#RandSM(N): A random N by N stochastic matrix
RandSM:=proc(N) local i:
[seq(RandPV(N),i=1..N)]:
end: