#Ross Berkowitz firstdraft.txt
#VERY Rough first draft of project on Shepp Urns.  More to come after work for
#the semester subsides.  Currently only contains the procedures tested.  Tests
#so far have been about effects on tweaks of the urn to the number of positive/negative balls
#needed to entice a player to play, and changes in variance of the various versions.
#Results and more procedures to come later.  
#Posting Permission Granted


HelpShepp:=proc():  print("RLarry(m,p,n), RLarryE(m,p,n), RLarryPos(m,p,n), RLarryPosE(m,p,n)"):
print("GLarry(m,p,x,y), GLarryE(m,p,x,y), R1LarryFastE1(p)"):  end:
#RLarry samples from an urn with replacement of negative balls only.  m negative, n negative
#simulates for N selections
RLarry:=proc(m,p,N) 
option remember:

if N=0 then
	RETURN(0):
fi:

if m=0 then
 RETURN(min(p,N)):
fi:

if p=0 then
 RETURN(0):
fi:

   
max(0,p/(m+p)*(1+RLarryE(m,p-1,N-1))+ m/(m+p)*(-1+RLarryE(m,p,N-1))):



end:

#companion to RLarryE, ensures no silly negatives taken.
RLarryE:=proc(m,p,N) 
option remember:
max(RLarry(m,p,N),0):
end:

#RLarryPos samples from urns with replacement of positive balls only (m neg, p pos, N trials)
RLarryPos:=proc(m,p,N) 
option remember:

if N=0 then
	RETURN(0):
fi:

if m=0 then
 RETURN(N):
fi:

if p=0 then
 RETURN(0):
fi:
0,p/(m+p)*(1+RLarryPosE(m,p,N-1))+ m/(m+p)*(-1+RLarryPosE(m-1,p,N-1)):
end:

#companion program to RLarryPos
RLarryPosE:=proc(m,p,N) 
option remember:
max(RLarryPos(m,p,N),0):
end:




RLarryE1Fast:=proc(p::posint) local x,y:
if p=1 then
	RETURN(0):
fi:

(1+RLarryE1Fast(p-1))-1/p:
end:


#GLarrry samples from urns where minus balls are worth -x, plus balls are worth y
GLarryE:=proc(m,p,x,y) 
option remember:
max(GLarry(m,p,x,y),0):
end:

#Companion program to GLarry
GLarry:=proc(m,p,x,y) 
option remember:

if m=0 then
 RETURN(p):
fi:

if p=0 then
 RETURN(0):
fi:
   
p/(m+p)*(y+GLarryE(m,p-1,x,y))+ m/(m+p)*(-x+GLarryE(m-1,p,x,y)):



end:











#March 11, 2013, Larry Shepp Urn
Help:=proc(): print(` Larry(m,p), LarryE(m,p) `): 
print(`LarryTale(N),beta(p) , LarryPGF(m,p,x)`):
print(`PrDist(m,p) `):
print(` Stat(P,x,K) `):
end:

#LarryE(m,p): input non.neg. integers m and p
#representing a box (urn) with m (-1)dollar bills
#and p 1-one dollar bills, outputs the
#expected gain of the game (where at any given time)
#you pick at random one or the other dollars bills
#Pr. of picking a dollar bill is p/(m+p) and 
#Pr. of picking a -1 dollar bill m/(m+p)
#And you quit whenever you want,
LarryE:=proc(m,p) 
option remember:
max(Larry(m,p),0):
end:

#Larry(m,p): input non.neg. integers m and p
#representing a box (urn) with m (-1)dollar bills
#and p 1-one dollar bills, outputs the
#expected gain of the game (where at any given time)
#you pick at random one or the other dollars bills
#Pr. of picking a dollar bill is p/(m+p) and 
#Pr. of picking a -1 dollar bill m/(m+p)
#And you quit whenever you want 
Larry:=proc(m,p) 
option remember:

if m=0 then
 RETURN(p):
fi:

if p=0 then
 RETURN(-1):
fi:
   
p/(m+p)*(1+LarryE(m,p-1))+ m/(m+p)*(-1+LarryE(m-1,p)):



end:

#LarryTable(N):
LarryTable:=proc(N) local m,p:
for p from N by -1 to 0 do
 lprint(seq(evalf(LarryE(m,p),3),m=0..N)):
od:




end:


#beta(p): The largest m for which LarryE(m,p) is >0
#beat(p)-p is asymtotically sqrt(p) times a (explicit) Constant
beta:=proc(p) local m:

for m from p while LarryE(m,p)>0 do od:
m-1:

end:


#LarrtPGF(m,p,x): The generalized "polynomial"
#whose exponets are numbers (rational) and
#where the coeff. of x^a is the prob. of exactly getting
#a dollars when you end the game (but following the
#stupid max. exp. strategy)
LarryPGF:=proc(m,p,x) 
option remember:

if LarryE(m,p)=0 then
  RETURN(1):
fi:

if m=0 then
 RETURN(x^p):
fi:


   
expand(p/(m+p)*x*LarryPGF(m,p-1,x)+ 
m/(m+p)/x*LarryPGF(m-1,p,x)):



end:


#PrDist(m,p): the list of pairs [i,p(i)] where
#the prob. of getting i is p(i)
PrDist:=proc(m,p) local olivia,x,i:

olivia:=LarryPGF(m,p,x):

[seq([i, coeff(olivia,x,i)],
i=ldegree(olivia,x)..degree(olivia,x))]:

end:

#Stat(P,x,K): inputs a prob. g.f. P in x and a pos. integer K
#outputs a list of length K where the first entry if the
#average, the second the variance, and the r-th is the s-called
#alpha coefficiient m[r]/m[2]^(r/2) where m[r] is the
#r-th moment about the mean
#try:
#evalf(Stat(((1+x)/2)^1000,x,10));
#Taken from homework problem, not done in class
Stat:=proc(P,x,K) local av,va,L,P1,r:

P1:=P/subs(x=1,P):

av:=subs(x=1,diff(P1,x)):
P1:=expand(P1/x^av):
P1:=expand(x*diff(P1,x)):
P1:=expand(x*diff(P1,x)):
va:=subs(x=1,P1):

L:=[av,va]:

for r from 3 to K do
P1:=expand(x*diff(P1,x)):
L:=[op(L),subs(x=1,P1)/va^(r/2)]:
od:

L:
end: