#Tim Naumovitz #hw22.txt #I give permission to post this ###################Problem 1############ #x(i)=p*x(i+1)+(1-p)*x(i-1) # #(1-p)*x(i-1)-x(i)+p*x(i+1)=0 # #(1-p)/x-1+p*x=0 # #p*x^2-x+(1-p)=0 # #(px-(1-p))(x-1)=0 # #x=1,(1-p)/p # #x(i)=c_1+c_2*((1-p)/p)^i # #c_1+c_2=0 # #c_1+c_2*((1-p)/p)^N=1 # #c_2*(((1-p)/p)^N-1)=1 # #c_2=(((1-p)/p)^N-1)^-1 # #c_1=-(((1-p)/p)^N-1)^-1 # #x(i)=-(((1-p)/p)^N-1)^-1+(((1-p)/p)^N-1)^-1*((1-p)/p)^i ###################Problem 2############ #> A := VW(p, 5); #> B := [seq(simplify(-1/(((1-p)/p)^5-1)+((1-p)/p)^i/(((1-p)/#p)^5-1)), i = 1 .. 4)]; #> evalb(A = B); # true #> A := VW(p, 10); #> B := [seq(simplify(-1/(((1-p)/p)^10-1)+((1-p)/p)^i/(((1-p)/#p)^10-1)), i = 1 .. 9)]; #> evalb(A = B); # true ###################Problem 3############ #x(i)=p*x(i+1)+(1-p)*x(i-1)+1 # #(1-p)*x(i-1)-x(i)+p*x(i+1)=-1 # #(1-p)/x-1+p*x=0 # #p*x^2-x+(1-p)=0 # #(px-(1-p))(x-1)=0 # #x=1,(1-p)/p # #PS(i)=? ################Problem 5############# VDW:=proc(p,N) local Du,eq,unk,i,sol: #Du[i]: expected duration of #a gambler's ruin with max. cap N and porb. p of winning #a dolllar unk:={ seq(Du[i],i=0..N)}: eq:={Du[1]=Du(2)+1, Du[N]=0, seq(Du[i]=p*Du[i+1]+(1-p)*Du[i-1]+1, i=2..N-1)}: sol:=solve(eq,unk): subs(sol, [seq(Du[i],i=1..N-1)]): end: ################Old Stuff############# #April 15, 2013, starting gambling theory Help:=proc(): print(` LC(p) , GRgame(p,s,N) `): print(`ManyGames(p,s,N,HowMany), VW(p,N) , VD(p,N) `): end: #LC(p): simulating a loaded coin whose #that returns true with prob. p, p rational LC:=proc(p) local a,b,ra: a:=numer(p): b:=denom(p): ra:=rand(1..b)(): evalb(ra<=a): end: #GRgame(p,s,N):inputs prob. p of success (a rational number) #starting capital, and exit capital N (if you are a winner) #once you have 0 dollars you must leave the casino #Returns, true (false) if you won (or lost) followed #by the number of rounds it took GRgame:=proc(p,s,N) local ca,du,toss: ca:=s: du:=0: while ca>0 and ca=N), du: end: #ManyGames(p,s,N,HowMany): repeats GRgame(p,s,N) #HowMany times, and returns the percentage of times #the gambler won, and the average duration if it won, #followed by the average duration if it lost followed #by the average duration ManyGames:=proc(p,s,N,HowMany) local DUw, DUl, DU, W,i,jake: W:=0: DUw:=0: DUl:=0: DU:=0: for i from 1 to HowMany do jake:=GRgame(p,s,N): if jake[1] then W:=W+1: DUw:=DUw+jake[2]: DU:=DU+jake[2]: else DUl:=DUl+jake[2]: DU:=DU+jake[2]: fi: od: if W=0 or W=HowMany then print(`This is never going to happen`): RETURN(FAIL): fi: evalf(W/HowMany), evalf(DUw/W), evalf(DUl/(HowMany-W)), evalf(DU/HowMany): end: #VW(p,N): inputs prob. of winning one round, p #and the exit capital N, outputs the List of #size N-1 whose i-th entry is the EXACT prob. #of winning if you currently posses i dollars VW:=proc(p,N) local w,eq,unk,i,sol: #w[i]: prob. of ultimately winning with current capital i #in Gambler's ruin name unk:={ seq(w[i],i=0..N)}: eq:={w[0]=0, w[N]=1, seq(w[i]=p*w[i+1]+(1-p)*w[i-1], i=1..N-1)}: sol:=solve(eq,unk): subs(sol, [seq(w[i],i=1..N-1)]): end: #VD(p,N): inputs prob. of winning one round, p #and the exit capital N, outputs the List of #size N-1 whose i-th entry is the EXACT EXpected Duration #of winning if you currently posses i dollars VD:=proc(p,N) local Du,eq,unk,i,sol: #Du[i]: expected duration of #a gambler's ruin with max. cap N and porb. p of winning #a dolllar unk:={ seq(Du[i],i=0..N)}: eq:={Du[0]=0, Du[N]=0, seq(Du[i]=p*Du[i+1]+(1-p)*Du[i-1]+1, i=1..N-1)}: sol:=solve(eq,unk): subs(sol, [seq(Du[i],i=1..N-1)]): end: