#DMB10.txt: Maple Code for Dr. Z.'s Dynamical Models in Biology class, Lecture 10, Oct. 3, 2025

with(linalg):

print(`For a list of the Main procedures type: Help10(); for help with a specific procedure type: Help10(ProcedureName); for example Help10(RandomMC);`):
print(`For a list of the other  procedures type: Help10a(); for help with a specific procedure type: Help10(ProcedureName); for example Help10(RandomMC);`):



Help10a:=proc()
if nargs=0 then
print(`The other procedures are: LC, LD  `):
else
Help6(args):
fi:
end:


Help10:=proc()
if nargs=0 then
print(`The available procedures are: EstPD, ExactPD, ForPD, PrLA, Simu, Story `):




elif nargs=1 and args[1]=DuLA then
print(`DuLA(L,p): Inputs the maximal capital L and prob. of  winning a dollar p outputs the list of length L-1 such that its n-th component is the expected duration of the game. Try:`):
print(`DuLA(10,1/2);`):

elif nargs=1 and args[1]=EstPD then
print(`EstPD(n,L,p,K): uses simulation to estimate the prob. of exiting a winner in a Gambler's Ruin with max. capital L dolalrs, initial capital n dolalrs`):
print(`prob. of winning a dollar at any round p, by running Simu(n,L,p) K times and averaging. The output is the pair of numbers`):
print(`[Estimated Prob. of winning, EstinatedDuration], Try:`):
print(`EstPD(5,10,1/2,100); `):

elif nargs=1 and args[1]=ExactPD then
print(`ExactPD(n,L,p): The exact value of  the prob. of winning and expected duration in a Gambler's ruin game with max. capital L and initial captail n`):
print(`with prob. of winning a dollar p, using the ready-made formulas. Try:`):
print(`ExactPD(50,100,49/100);`):

elif nargs=1 and args[1]=ForPD then
print(`ForPD(n,L,p): The formula for the prob. and expected duration in a grambler's ruin with max. capital L and starting capital n and p less than 1. p can be symbolic`):
print(`Try:`):
print(`ForPD(50,100,p):`):

elif nargs=1 and args[1]=LC then
print(`LC(p): inputs a probability (given as rational number of the form a/b) outputs 1 with prob. p and -1 with prob. 1-p. Try:`):
print(`LC(2/3);`):
print(`add(x[C(2/3)],i=1..300));`):

elif nargs=1 and args[1]=LD then
print(`LD(L): Inputs a finite discrete prob. distribution`):
print(`Try: LD([1/2,1/5,3/10]);`):


elif nargs=1 and args[1]=PrLA then
print(`PrLA(L,p): Inputs the maximal capital L and prob. of  winning a dollar p outputs the list of length L-1 such that its n-th component is the prob. of exiting a winner. Try:`):
print(`PrLA(10,1/2);`):

elif nargs=1 and args[1]=Simu then
print(`Simu(n,L,p): Simulating ONE visit to the casino with maximal capital L dollars and initial capital n dollars with prob. of winning a dollar`):
print(`p, terse version. Try:`):
print(`Simu(5,10,1/2);`):

elif nargs=1 and args[1]=Story then
print(`Story(n,L,p): Simulating ONE visit to the casino with maximal capital L dollars and initial capital n dollars with prob. of winning a dollar`):
print(`p, verbose version. Try:`):
print(`Story(5,10,1/2);`):

else
print(`There is no Help for`, args):
fi:

end:


#LC(p): inputs a probability (given as rational number of the form a/b) outputs 1 with prob. p and -1 with prob. 1-p. Try:
#LC(2/3);
#add(x[C(2/3)],i=1..300));
LC:=proc(p) local a,b,i:
if not type(p,fraction) and p>0 and p<1 then
 RETURN(FAIL):
fi:
a:=numer(p):
b:=denom(p):
i:=rand(1..b)():
if i<=a then
1:
else
-1:
fi:
end:


#Story(n,L,p): Simulating ONE visit to the casino with maximal capital L dollars and initial capital n dollars with prob. of winning a dollar
#p, verbose version
Story:=proc(n,L,p) local i,n1,c:

print(`At the start you had`, n, `dollars , the max. capital is`,L, `dollars. `):

n1:=n:

for i from 1 while n1>0 and n1<L do
c:=LC(p):
n1:=n1+c:
print(``):
print(`Now you at round`, i):
print(``):
if c=1 then
print(`Congratulations, you won a dollar, you now have`, n1, `dollars `):
else
print(`Sorry , you lost a dollar, you now have`, n1, `dollars `):
fi:
od:

if n1=0 then
 print(`Poor you, now you are broke, and it your stay in the casion took`, i-1, `rounds `):
 elif n1=L then
  print(`Congratulations, now you are rich, and it your stay in the casion took`, i-1, `rounds `):
  else
  print(`Something is wrong`):
   RETURN(FAIL):
fi:

if n1=0 then
 RETURN([0,i-1]):
 else
 RETURN([1,i-1]):
fi:
end:




#Simu(n,L,p): Simulating ONE visit to the casino with maximal capital L dollars and initial capital n dollars with prob. of winning a dollar
#p, terse version
Simu:=proc(n,L,p) local i,n1,c:

n1:=n:

for i from 1 while n1>0 and n1<L do
c:=LC(p):
n1:=n1+c:
od:

if n1=0 then
 RETURN([0,i-1]):
 else
 RETURN([1,i-1]):
fi:
end:



#EstPD(n,L,p,K): uses simulation to estimate the prob. of exiting a winner in a Gambler's Ruin with max. capital L dolalrs, initial capital n dolalrs
#prob. of winning a dollar at any round p, by running Simu(n,L,p) K times and averaging. The output is the pair of numbers
#[Estimated Prob. of winning, EstinatedDuration], Try:
#EstPD(5,10,1/2,100):

EstPD:=proc(n,L,p,K) local i:
evalf(add(Simu(n,L,p),i=1..K)/K):
end:

#ExactPD(n,L,p): The exact value of  the prob. of winning and expected duration in a Gambler's ruin game with max. capital L and initial captail n
#with prob. of winning a dollar p, using the ready-made formulas. Try:
#ExactPD(100,49/100);
ExactPD:=proc(n,L,p) local lu:

if not (type(p,fraction) and p>0 and p<1 and type(n,integer) and type(L, integer) and n<=L and n>=0 ) then
RETURN(FAIL):
fi:
 
if p=1/2 then
RETURN([n/L,n*(L-n)]):
elif p<1/2 then
RETURN(evalf([ (1-((1-p)/p)^n)/ (1-((1-p)/p)^L), n/(1-2*p)- L/(1-2*p)*(1-((1-p)/p)^n)/(1-((1-p)/p)^L)])):
else
lu:=ExactPD(L-n,L,1-p):
RETURN([1-lu[1],lu[2]]):
:
fi:

end:

#ForPD(n,L,p): The formula for the prob. and expected duration in a grambler's ruin with max. capital L and starting capital n and p less than 1. p can be symbolic
#Try:
#ForPD(50,100,p):
ForPD:=proc(n,L,p)
[ (1-((1-p)/p)^n)/ (1-((1-p)/p)^L), n/(1-2*p)- L/(1-2*p)*(1-((1-p)/p)^n)/(1-((1-p)/p)^L)]:
end:


#PrLA(L,p): Inputs the maximal capital L and prob. of  winning a dollar p outputs the list of length L-1 such that its n-th component is the prob. of exiting a winner. Try:
#PrLA(10,1/2);
PrLA:=proc(L,p) local a,var,eq,i:

var:=[seq(a[i],i=1..L-1)]:
eq:={seq(a[i]=p*a[i+1]+(1-p)*a[i-1],i=1..L-1)}:
eq:=subs({a[0]=0,a[L]=1},eq):
subs(solve(eq,{op(var)}),var):
end:



#DuLA(L,p): Inputs the maximal capital L and prob. of  winning a dollar p outputs the list of length L-1 such that its n-th component is the expected duration of the game.  Try:
#DuLA(10,1/2);
DuLA:=proc(L,p) local a,var,eq,i:

var:=[seq(a[i],i=1..L-1)]:
eq:={seq(a[i]=p*a[i+1]+(1-p)*a[i-1]+1,i=1..L-1)}:
eq:=subs({a[0]=0,a[L]=0},eq):
subs(solve(eq,{op(var)}),var):
end:



Nor1:=proc(L): L/add(L):end:

#LD(L): Inputs a finite discrete prob. distribution
LD:=proc(L) local L1,N,i,r,j:
L1:=Nor1(L):


L1:=L1*lcm(seq(denom(L1[i]),i=1..nops(L1)) ):
N:=add(L1):
r:=rand(1..N)():

for i from 1 while add(L1[j],j=1..i)<r do od:
i:

end: