#C17.txt: March 26, 2020 virtual class
#Exact expressions for the probability of winning  expected duration and higher moments of duration
Help:=proc(): print(`  StatAnal(P,t,K) , ProbWinning(N,p), ExpDuration(N,p), PGF(N,p,t), `): 
print(`  ProbWinningCF(N,n,p)  , ExpDurationCF(N,n,p), PGFcf(N,n,p,t) `):
end:
with(Statistics):

#From C16.txt:
#modified March 26, to get around "division by 0" by taking the limit as t goes to 1
#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 limit(P,t=1)=0 then
 RETURN(FAIL):
else
 P1:=P/limit(P,t=1):
fi:

#mu is the average
mu:=limit(diff(P1,t),t=1):
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 normal(limit(P1,t=1))<>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),limit(P1,t=1)]:
od:

#Finally we adjust the third-through the K-th moment by dividing the k-th moment 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:


#ProbWinning(N,p): inputs a positive integer N and a prob. p (either numeric or symbolic),
#outputs a list of length N-1 whose i-th entry is the probability of exiting a winner in
#Gambler's ruin game where you start with a capital of i dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ProbWinning(10,1/2);
#ProbWinning(10,1/3);
#ProbWinning(10,2/3);
#ProbWinning(10,p);
ProbWinning:=proc(N,p) local eq,var,x,i:

#We call the prob. of winning if you are currently with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently in possession of i dollars, then with probability p next round you would have i+1 dollars
#and with prob. 1-p next time you would have i-1 dollars. By CONDITIONING on the event that you
#currently have i dollars, we get the equation
#x[i]=p*x[i+1]+(1-p)*x[i-1]. 
#Also if right now you have 0 dollars your prob. of winning is 0
#and if you currently have N dollars than your prob. of winning the game is 1 (THESE ARE THE BOUNDARY CONDITIONS)

eq:={seq(x[i]=p*x[i+1]+(1-p)*x[i-1],i=1..N-1),x[0]=0,x[N]=1}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

[seq(subs(var,x[i]),i=1..N-1)]:

end:




#ExpDuration(N,p): inputs a positive integer N and a prob. p (either numeric or symbolic),
#outputs a list of length N-1 whose i-th entry is the expected duration of the a
#Gambler's ruin game where you start with a capital of i dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ExpDuration(10,1/2);
#ExpDuration(10,1/3);
#ExpDuration(10,2/3);
#ExpDuration(10,p);
ExpDuration:=proc(N,p) local eq,var,x,i:

#We call the expected remaining duration if you currently have  with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently in possession of i dollars, then with probability p next round you would have i+1 dollars
#and with prob. 1-p next time you would have i-1 dollars. By CONDITIONAL EXPECTATION on the event that you
#currently have i dollars, we get the equation
#x[i]=p*x[i+1]+(1-p)*x[i-1]+1
#(Note that the +1 accounts that the duration today is one more than the duration tomorrow)
#Also if right now you have 0 dollars or N dollars your duration (and hence expected duration) is 0
eq:={seq(x[i]=p*x[i+1]+(1-p)*x[i-1]+1,i=1..N-1),x[0]=0,x[N]=0}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

[seq(subs(var,x[i]),i=1..N-1)]:

end:




#PGF(N,p,t): inputs a positive integer N and a prob. p (either numeric or symbolic),
#and a variable name t,
#outputs a list of length N-1 whose i-th entry is the PROBABILITY GENERATING FUNCTION FOR DURATION
#(the rational function in t whose coefficient of t^i is the probability that the game will end exactly after i round)
#In a Gambler's ruin game where you start with a capital of i dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#PGF(10,1/2,t);
#PGF(10,1/3,t);
#PGF(10,2/3,t);
#PGF(10,p,t);
PGF:=proc(N,p,t) local eq,var,x,i:

#We call the prob. generating functions for the remaining duration if you currently have  with i dollars, x[i], these are the unknowns
var:={seq(x[i],i=0..N)}:

#If you currently have i dollars, then with probability p next round you would have i+1 dollars
#and with prob. 1-p next time you would have i-1 dollars. By CONDITIONAL EXPECTATION on the event that you
#currently have i dollars, we get the equation
#x[i]=t*(p*x[i+1]+(1-p)*x[i-1])
#(Note that the multiplication by t  accounts that the duration today is one more than the duration tomorrow)
#Also if right now you have 0 dollars or N dollars your duration the prob. generating function is t^0=1 
#(I will definitely take you 0 rounds to finish the game)
eq:={seq(x[i]=t*(p*x[i+1]+(1-p)*x[i-1]),i=1..N-1),x[0]=1,x[N]=1}:

#We now have N+1 equations and N+1 unknowns. Let's solve this linear system, and rename the
#solutions var

var:=solve(eq,var):

#var is a set of assignments. To get the solutions  yous substitute into x[i]. This is the output of
#the procedure

normal([seq(subs(var,x[i]),i=1..N-1)]):

end:



#ProbWinningCF(N,n,p): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic),
#outputs the CLOSED-FORM EXPRESSION for the prob. of winning in
#Gambler's ruin game where you start with a capital of n dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ProbWinningCF(N,n,p);
ProbWinningCF:=proc(N,n,p) local x,F,A,A0:
#We let Maple solve the recurrence equation x(n)=(1-p)*x(n-1)+p*x(n+1), but it can only handle
#INITIAL conditions, and so far we have no clue what x(1) is, so we call it A

F:=rsolve({x(n)=(1-p)*x(n-1)+p*x(n+1),x(0)=0, x(1)=A},x(n)):

#Now we use the fact that x(N)=1, and solve for A
A0:=solve(subs(n=N,F)=1,A):

#Now we plug that expression for A0 into F

normal(simplify(subs(A=A0,F))):

end:


#ExpDurationCF(N,n,p): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic),
#outputs the CLOSED-FORM EXPRESSION for the Expected duration
#Gambler's ruin game where you start with a capital of n dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#ExpDurationCF(N,n,p);
ExpDurationCF:=proc(N,n,p) local x,F,A,A0:
#We let Maple solve the recurrence equation x(n)=(1-p)*x(n-1)+p*x(n+1)+1, but it can only handle
#INITIAL conditions, and so far we have no clue what x(1) is, so we call it A

F:=rsolve({x(n)=(1-p)*x(n-1)+p*x(n+1)+1,x(0)=0, x(1)=A},x(n)):

#Now we use the fact that x(N)=0, and solve for A
A0:=solve(subs(n=N,F)=0,A):

#Now we plug that expression for A0 into F

normal(simplify(subs(A=A0,F))):

end:


#PGCcf(N,n,p,t): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic), and a variable name (symbol) t;
#outputs the CLOSED-FORM EXPRESSION for the Prob. Generating function in t, of the duration
#Gambler's ruin game where you start with a capital of n dollars, and then at each round
#you win a dollar with prob. p and lose a dollar with prob. 1-p. Try:
#PGFcf(N,n,p,t);
PGFcf:=proc(N,n,p,t) local x,F,A,A0:
#We let Maple solve the recurrence equation x(n)=t*((1-p)*x(n-1)+p*x(n+1)), but it can only handle
#INITIAL conditions, and so far we have no clue what x(1) is, so we call it A


F:=rsolve({x(n)=t*((1-p)*x(n-1)+p*x(n+1)),x(0)=1, x(1)=A},x(n)):

#Now we use the fact that x(N)=1, and solve for A
A0:=solve(subs(n=N,F)=1,A):

#Now we plug that expression for A0 into F

normal(simplify(subs(A=A0,F))):

end: