#C16.txt: March 23, 2020 virtual class
#The Gambler's ruin's problem

Help:=proc(): print(` GRsim1(N,p,n) , GRsim(N,p,n,M,s,t)  `): end:
with(Statistics):

#GRsim1(N,p,n): Simulating ONE Gambler's ruin game with initial capital of n dollars, 
#At each turn the gambler's win a dollar with probability p and loses a dollar with probability 1-p
#the game ends as soon as it reaches N dollars or 0 dollars (when the gambler's is ruined).
#The output is a pair [0,m] or [1,m], according to whether the gambler's left the casino broke or
#rich, and m is the number of rounds. Try:
#GRsim1(10,1/2,5);
#GRsim1(10,10/21,5);
GRsim1:=proc(N,p,n) local x,U,coin,count:
U:=RandomVariable(Bernoulli(p)):

x:=n:
count:=0:
while x>0 and x<N do

count:=count+1:
coin:=Sample(U,1)[1]:

 if coin=0. then
  x:=x-1:
 else
  x:=x+1:
 fi:
od:

if x=0 then
 RETURN([0,count]):
elif x=N then
 RETURN([1,count]):
else
 print(`Something is wrong`):
 RETURN(FAIL):
fi:


end:

#GRsim(N,p,n,M,s,t): Inputs a positive integer N, a rational number p between 0 and 1, and a pos. integer n<=N
#as well as symbols (variables) s and t, outputs
#The Empirical probability generating function of running GR1(N,p,N) M times
#it returns a polynomial of degree 1 in s, i.e. of the form F(t,s)=P0(t)+s*P1(t) where
#PO(t) records the losing cases and P1(t) the winning cases according to duration.
#In particular the coefficient of s in F(1,s) is the prob. of winning
#the coefficient of s^0 (the constant term) in F(1,s) is the prob. of losing
#and F(1,t) is the prob. generating function for the duration of the game (regardless whether you won or lost)
#Try:
#GRsim(40,1/2,20,100,s,t);
GRsim:=proc(N,p,n,M,s,t) local i,su,G:

su:=0:

for i from 1 to M do
 G:=GRsim1(N,p,n):
 su:=su+s^G[1]*t^G[2]:
od:
su/M:
end: