#Nathan Fox #Homework 22 #I give permission for this file to be posted online ##Read old files read(`C22.txt`): #Help procedure Help:=proc(): print(` VWS(N) , VWcomp(N) , VWcheck(N) , CheckProblem6() `): print(` VDS(N) , VDcomp(N) , VDcheck(N) , VWD(N) , VWDinterp(N) `): end: ##Problem 1 #See hw22.pdf ##Problem 2 #Use my formula from problem 1 for all i from 0 to N VWS:=proc(p,N) local i,alpha,beta: alpha:=(1+sqrt(1+4*p*(p-1)))/(2*p): beta:=(1-sqrt(1+4*p*(p-1)))/(2*p): [seq((alpha^i-beta^i)/(alpha^N-beta^N),i=1..N-1)]: end: #Compare VW(p,N) and VWS(N) VWcomp:=proc(N) local i,vw,vws,p,vwdiff: vw:=VW(p,N): vws:=VWS(p,N): vwdiff:={seq(simplify(vw[i]-vws[i]),i=1..nops(vw))}: evalb(vwdiff={0}): end: #Run VWcomp from 2 to N VWcheck:=proc(N) local i: evalb({seq(VWcomp(i),i=2..N)}={true}): end: #VWcheck(100) was taking way too long #VWcheck(20) returned true, as it should have ##Problem 3 #See hw22.pdf ##Problem 4 #Use my formula from problem 3 for all i from 0 to N VDS:=proc(p,N) local i,alpha,beta: alpha:=(1+sqrt(1+4*p*(p-1)))/(2*p): beta:=(1-sqrt(1+4*p*(p-1)))/(2*p): [seq((N/(2*p-1))*(alpha^i-beta^i)/(alpha^N-beta^N)+i/(1-2*p), i=1..N-1)]: end: #Compare VD(p,N) and VDS(N) VDcomp:=proc(N) local i,vd,vds,p,vddiff: vd:=VD(p,N): vds:=VDS(p,N): vddiff:={seq(simplify(vd[i]-vds[i]),i=1..nops(vd))}: evalb(vddiff={0}): end: #Run VDcomp from 2 to N VDcheck:=proc(N) local i: evalb({seq(VDcomp(i),i=2..N)}={true}): end: #VDcheck(100) was taking way too long #VDcheck(20) returned true, as it should have ##Problem 5 #VWD(p,N): inputs a positive integer N and a prob. of winning #(a single coin toss) p (numeric or symbolic) and outputs the #expected duration of a game conditioned on winning the game. VWD:=proc(p,N) local du,eq,unk,i,sol,probi,vw: vw:=VW(p,N): vw:=[op(vw),1]: #Conditional probability of increasing with capital i #given that you will win (computed using Bayes' Theorem) probi:=proc(i) vw[i+1]*p/vw[i]: end: unk:={seq(du[i],i=1..N)}: eq:={du[1]=1+du[2], du[N]=0, seq(du[i]=1+probi(i)*du[i+1]+(1-probi(i))*du[i-1], i=2..N-1)}: sol:=solve(eq,unk): subs(sol, [seq(du[i],i=1..N-1)]): end: ##Problem 6 with(CurveFitting): #Interpolate VWD(1/2,N) VWDinterp:=proc(N) local vwd,i: vwd:=[op(VWD(1/2,N)),0]: PolynomialInterpolation([seq([i,vwd[i]],i=1..nops(vwd))],i): end: #Based on this, it appears that VWD(1/2,N) = (N^2-i^2)/3 ##Problem 7 #Note that the conditional probability of increasing with capital #i given that you will win=probi(i)=(i+1)/(2*i) if p=1/2 (using #the fact that VW(p,N)[i]=i/N). Also, it is easy to see that #(N^2-i^2)/3 will satisfy the desired boundary conditions. So, #all we need to check is that it satisfies the non-boundary #conditions. In other words, the following procedure must #return true CheckProblem6:=proc() local N,i: evalb((N^2-i^2)/3=simplify(1+((i+1)/(2*i))*((N^2-(i+1)^2)/3) +((i-1)/(2*i))*((N^2-(i-1)^2)/3))): end: #It returns true. QED