#C16.txt: March 23, 2020 virtual class
##Discrete Prob. generating functions, and starting The Gambler's ruin's problem using simulation

Help:=proc(): print(`  StatAnal(P,t,K), GRsim1(N,p,n) , GRsim(N,p,n,M,s,t) ,  PH(N,p,n), Roullete(L)  `): end:
with(Statistics):


#StatAnal(P,t,K): Inputs a polynomial (or even a function whose Taylor expansion has positive coefficients)
#describing the (not-necessarily normalized) probability generating function of some discrete random variable
#outputs  a list of length K whose
#(i) First entry is the expectation (alias mean, alias average)
#(ii) Second entry is the variance (aka as second moment about the mean)
#(iii) third-through K entries are the scaled moments about the mean
#in particular the third entry is the skewness and the fourth moment is the kurtosis. For example, try:
#StatAnal((1+t)^1000,t,6);
StatAnal:=proc(P,t,K) local P1,L,k,mu:

#First let's normalize it (just in case)
if subs(t=1,P)=0 then
 RETURN(FAIL):
else
 P1:=P/subs(t=1,P):
fi:

#mu is the average
mu:=subs(t=1,diff(P1,t)):
L:=[mu]:

#Now let's CENTRALIZE, get the prob. gen. function for the "deviation from the mean"
P1:=P1/t^mu:

#We will now apply the operation td/dt repeatedly and plug-in t=1
P1:=t*diff(P1,t):

#The expectation of the "deviation from the mean" should be 0, let's check it
if subs(t=1,P1)<>0 then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

for k from 2 to K do
P1:=t*diff(P1,t):
#We append the current moment-about the mean
L:=[op(L),subs(t=1,P1)]:
od:

#Finally we adjust the third-through the K-th moment by dividing the k-th momend by
#the L[2]^(k/2) (the k-th power of the standard deviation)

if L[2]=0 then
 print(`The variance is 0, the unscaled moments are`):
 print(L):
 RETURN(FAIL):
else
[L[1],L[2],seq(L[k]/L[2]^(k/2),k=3..K)]:
fi:

end:

#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:


#Added after class, code by Victoria Chayes
#PH(N,p,n): Drawing  ONE Gambler's ruin game with initial capital of n dollars, 
#PH(10,1/2,5);
#PH(10,10/21,5);
PH:=proc(N,p,n) local x,U,coin,count,L:
U:=RandomVariable(Bernoulli(p)):

x:=n:
count:=0:

L:=[]:

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:
 L:=[op(L),[count,x]]:
od:
plot(L):
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 GRsim1(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:

#Added after class, code by Victoria Chayes
#PH(N,p,n): Drawing  ONE Gambler's ruin game with initial capital of n dollars, 
#PH(10,1/2,5);
#PH(10,10/21,5);
PH:=proc(N,p,n) local x,U,coin,count,L:
U:=RandomVariable(Bernoulli(p)):

x:=n:
count:=0:

L:=[]:

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:
 L:=[op(L),[count,x]]:
od:
plot(L):
end:


#ADDED AFTER CLASS
#Roulette(L): Inputs a list of pairs [ProbOfGain,Gain], returns Gain with probability ProbOfGain.
#To get the Bernoulli distribution with prob. p of winning a dollar and prob. 1-p of losing a dollar
# do Roullete([[1/2,-1],[1/2,1]]);
Roulette:=proc(L) local Probs,Vals, U,i:
Probs:=[seq(L[i][1],i=1..nops(L))]:

if not add(Probs)=1 then
 RETURN(FAIL):
fi:

Vals:=[seq(L[i][2],i=1..nops(L))]:

U:=RandomVariable(ProbabilityTable(Probs)):
Vals[round(Sample(U,1)[1])]:
end: