#C18.txt: March 28, 2020 Help:=proc(): print(` StatAnalPure(P,t,K), StatAnalPureL(P,t,K), GP1(L,n,d), GP(L,n), MomsGR1(n,N,K) , MomsGR(n,N,K,N1,N2) `): end: ###Stuff from C17.txt #C17.txt: March 26, 2020 virtual class #Exact expressions for the probability of winning expected duration and higher moments of duration Help17:=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 subs(t=1,P)=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: ###End from C17.txt #StatAnalPureL(P,t,K): The expecation, variance, and the UNSCALDED moments about the mean #the first entry is the expectation, the second the variance, and then the third, ..., up to the #K-th moment about the mean. #It uses limit(...,t=1) instead of subs(t=1, ...) MUCH SLOWER #Try: #StatAnalPureL((1+t)^n,t,10); StatAnalPureL:=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(subs(P1),t=1)]: od: L: end: #StatAnalPure(P,t,K): The expecation, variance, and the UNSCALDED moments about the mean #the first entry is the expectation, the second the variance, and then the third, ..., up to the #K-th moment about the mean. Try: #StatAnalPure((1+t)^n,t,10); StatAnalPure:=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 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),subs(t=1,P1)]: od: L: end: #GP1(L,n,d): inputs a list L , a variable name n, and a positive integer d, guesses a polynomial of degree <=d in n, #let's call it P(n) such that P(i)=L[i]. The length of L must be at least d+3. If it fails it returns FAIL. #Try: #GP1([seq(i^3,i=1..10)],n,3); GP1:=proc(L,n,d) local a,i,var,eq,P: #Any d+1 data points can fit a polynomial of degree d. We want at least 2 "confirmation" points if nops(L)<=d+2 then print(`The list must be of length at least`, d+3): RETURN(FAIL): fi: #We define a generic polynomial P of n, of degree d, with UNDETERMINED coefficients. Our goal is to determine them P:=add(a[i]*n^i,i=0..d): #This is the set of unknowns with d+1 unknowns var:={seq(a[i],i=0..d)}: #We set-up a system with d+3 equations fitting the polynomial to the data eq:={seq(subs(n=i,P)-L[i],i=1..d+3)}: #Maple tries to solve the system, we rename var the solution var:=solve(eq,var): #If var=NULL then there is no solution, in other words, no polynomial of degree d fits the data if var=NULL then RETURN(FAIL): fi: #Now we plug the solution var into the generic polynomial P, making it an actual polynomial (that fits the data) P:=subs(var,P): #If nops(L) is larger than d+3, we check that the polynomial fits the remaining data for i from d+4 to nops(L) do if expand(subs(n=i,P))-L[i]<>0 then print(P, `does not agree with n=`, i): RETURN(FAIL): fi: od: #If it agrees we return the polynomial P: end: #GP(L,n): inputs a list L , a variable name n, guesses a polynomial of degree <=nops(L)-3 in n, #let's call it P(n) such that P(i)=L[i] for all i. If it fails it returns FAIL. #Try: #GP([seq(i^3,i=1..10)],n); GP:=proc(L,n) local d,P: for d from 0 to nops(L)-3 do #We start with d=0 and keep trying with higher and higher degree P:=GP1(L,n,d): if P<>FAIL then RETURN(P): fi: od: #If there is no such polynomial of degree <=nopd(L)-3 it returns FAIL FAIL: end: #MomsGR1(n,N,K): inputs a variable name n, and a positive integer N (sufficiently large) and a positive integer K #tries to find polynomial expression in n for the first K moments of the r.v. duration of a Gambler's #Ruin game with a fair coin with exit capital N and starting capital n. If it can only partially do it, #it returns what it can. Otherwise it returns all the K moments #Try: #MomsGR1(n,20,6); MomsGR1:=proc(n,N,K) local t,L,i,M,j,P: #We use the semi-numerical, semi-symbolic procedure (that uses simple linear algebra) #to find the list of N-1 prob. generating functions for the duration stating with n dollars L:=PGF(N,1/2,t): #We rename L to be the list of lists each of length K, describing the first K moments for each of n=1..N-1 L:=[seq(StatAnalPure(L[i],t,K),i=1..nops(L))]: M:=[]: #We try to guess polynomial expressions for the first K moments one by one, each is a polynomial in n for i from 1 to K do P:=GP([seq(L[j][i],j=1..nops(L))],n): #If we failed to get a polynomial expression for the i-th moment, we return what we got if P=FAIL then RETURN(M): fi: #otherwise we append the new moment M:=[op(M),P]: od: #If all goes well we got polynomial expressions in n for the first K moments of the Gambler's ruin with exit capital N (N numeric) M: end: #MomsGR(n,N,K,N1,N2): tries to conjecture polynomial expressions in n and N for the first K moments, #N=N1 to N=N2 is where we collect data MomsGR:=proc(n,N,K,N1,N2) local L1,P,N3,M,i: #L is the list of lists of length K from N1 to N2, we first make sure that with N=N1 we get the full K moments L1:=MomsGR1(n,N1,K): if nops(L1)