#May 1, 2014, C27.txt, more gambler's ruin
Help:=proc(): print(` PGFam(N,p,t), Moms(N,p,R) `): 
print(`Momsi(N,R,i) , MomsiN(N,R,i) , AStmomsiN(R,x)`):
 end:


###old stuff
#C26.txt, Gambler's ruin, April 28, 2014 
Help26:=proc(): print(`Die1(p), GR(i,N,p), ExpD(N,p)`): 
print(`ExpDf(N,p) , PGF(N,p,t), PGFf(N,p,t), GuessP(L,x) `):
end:

#Die1: inputs p rational between 0 and 1
#and outputs 1 with prob. p and 0 with prob. 1-p
Die1:=proc(p) local a,b,ra:
a:=numer(p):
b:=denom(p):
ra:=rand(1..b)():
if ra<=a then
 1:
else
 -1:
fi:

end:


#GR(i,N,p): simulates a gambler's ruin session
#staring with i dollars, ending with 0 or N dollars
#and with prob. of winning a dollar =p
GR:=proc(i,N,p) local x,j:

x:=i:

for j from 0 while (x>0 and x<N) do
 x:=x+Die1(p):
od:

x,j:

end:

#ExpD(N,p): Let f(i) be the expected duration of the
#game starting at i, outputs the vector
#[f(1),f(2),..., f(N-1)] of expected durations with
#prob. p of winning one dollar
#f(i)=p*f(i+1)+(1-p)*f(i-1)+1, i=1..N-1, f(0)=0, f(N)=0
ExpD:=proc(N,p) local eq,var,f,i,var1:
var:=[seq(f[i],i=1..N-1)]:
eq:={seq(f[i]=p*f[i+1]+(1-p)*f[i-1]+1,i=1..N-1)}:
eq:=subs({f[0]=0,f[N]=0},eq):
var1:=solve(eq,convert(var,set)):
subs(var1,var):
end:


#ExpD(N,p): fast version;
#Let f(i) be the expected duration of the
#game starting at i, outputs the vector
#[f(1),f(2),..., f(N-1)] of expected durations with
#prob. p of winning one dollar
#f(i)=(f(i-1)-(1-p)*f(i-2)-1)/p, i=1..N-1, f(0)=0, f(N)=0
ExpDf:=proc(N,p) local f,i,a,a1:
f[1]:=a:
f[0]:=0:
for i from 2 to N do
f[i]:=normal((f[i-1]-(1-p)*f[i-2]-1)/p):
od:

a1:=solve(f[N]=0,a):

[seq(subs(a=a1,f[i]),i=1..N-1)]:

end:


#PGF(N,p,t): Let f(i,t) be the prob. generating function of the
#a rational function in t whose Maclaurin expansion's
#coeff. of t^j
#gives the exact prob. of finishing in j steps
#game starting at i, outputs the vector
#[f(1),f(2),..., f(N-1)] of expected durations with
#prob. p of winning one dollar
#f(i,t)=t*(p*f(i+1,t)+(1-p)*f(i-1,t)), i=1..N-1, f(0)=0, f(N)=0
PGF:=proc(N,p,t) local eq,var,f,i,var1:
var:=[seq(f[i],i=1..N-1)]:
eq:={seq(f[i]=t*(p*f[i+1]+(1-p)*f[i-1]),i=1..N-1)}:
eq:=subs({f[0]=1,f[N]=1},eq):
var1:=solve(eq,convert(var,set)):
subs(var1,var):
end:


#PGFf(N,p,t): fast version;
#f(i,t)=t*(p*f(i+1,t)+(1-p)*f(i-1,t))
#f(i,t)= (f(i-1,t)/t-(1-p)*f(i-2,t))/p
PGFf:=proc(N,p,t) local f,i,a,a1:
f[1]:=a:
f[0]:=1:
for i from 2 to N do
f[i]:= normal((f[i-1]/t-(1-p)*f[i-2])/p):
od:

a1:=solve(f[N]=1,a):

[seq(subs(a=a1,f[i]),i=1..N-1)]:

end:

# begin Added after class
GuessP1:=proc(L,x,d) local i,P,a,var,eq:


if nops(L)-d<2 then
RETURN(FAIL):
fi:

P:=add(a[i]*x^i,i=0..d):

var:={seq(a[i],i=0..d)}:
eq:={seq(subs(x=i,P)=L[i],i=1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

subs(var,P):

end:

#GuessP(L,x): guesses a polynomial expression in x for the list L  such that L[i]=P(i)
GuessP:=proc(L,x)  local d,gu:

for d from 0 to nops(L)-3 do
gu:=GuessP1(L,x,d):
if gu<>FAIL then
  RETURN(gu):
fi:
od:

FAIL:

end:
#end Added after class
#end of old stuff


#PGFam(N,p,t): the list of pgf for duration of
#a gambler's ruin game with exit capital N
#about the mean
PGFam:=proc(N,p,t) local Li, Liav:
Li:=PGFf(N,p,t):

Liav:=subs(t=1,diff(Li,t)):

[seq(Li[i]/t^Liav[i],i=1..nops(Li))]:

end:


#Moms(N,p,R): the list of length N-1 whose
#ith item is the R-th moment
Moms:=proc(N,p,R) local t,Li,i1:
Li:=PGFam(N,p,t):

for i1 from 1 to R do
 Li:=expand(t*diff(Li,t)):
od:

subs(t=1,Li):

end:

#Momsi(N,R,i): The explicit polynomial expression in
#the symbol i for the R-th moment of the r.v.
#remainign duration of the FAIR game
Momsi:=proc(N,R,i) :

GuessP(Moms(N,1/2,R),i):



end:

#MomsiN(N,R,i): a guessed polynomial in symbolic i and N
#for the R-th moment of the r.v. "duration of a game
#staring at i dollars and ending at either 0 or N, with
#a FAIR coin
MomsiN:=proc(N,R,i) local N1:
option remember:
expand(subs(N=N-19,GuessP([seq(Momsi(N1,R,i),N1=20..30+3*R)],N))):



end:

#AStmomsiN(R,x): Asympotic standarized R-th moment
#if you start with x*N dollars to order k
AStmomsiN:=proc(R,x) local N,f,i,y:

f:=MomsiN(N,R,i)/MomsiN(N,2,i)^(R/2):

f:=subs(i=x*N,f):

if R mod 2=0 then

f:=normal(subs(N=1/y,f)):

coeff(taylor(f,y=0,7),y,0):

else

f:=simplify(normal(subs(N=1/y,f))^2):

coeff(taylor(f,y=0,7),y,0)^(1/2):

fi:




end: