SM-5-Exercise3

4135 days ago by MathFest

#Strodt Generating Function generalized to finite probability distribution #Input p is an array listing locations of point masses #Input w is an array listing weights of point masses #Note: Make sure P and w have the same length x,t=var('x,t') def FiniteStGF(p,w,x,t): Q=0 L=len(p) for j in range(L): Q=Q+exp(p[j]*t)*w[j] return exp(x*t)/Q 
       
#Pulls the Strodt polynomials out from the generating function #Again p and w are lists of locations and weights of point masses, resp. def P(p,w,n,x): g=taylor(FiniteStGF(p,w,x,t),t,0,n) out=g.coefficient(t^n)*factorial(n) if n==0: out=1 return out 
       
#Listing the Strodt polynomials for the first few values of n with a particular symmetric distribution #Note: The denominators are all powers of 2---Why? for n in range(11): P([0,1/2,1],[1/4,1/2,1/4],n,x) 
       
1
x - 1/2
x^2 - x + 1/8
x^3 - 3/2*x^2 + 3/8*x + 1/16
x^4 - 2*x^3 + 3/4*x^2 + 1/4*x - 1/16
x^5 - 5/2*x^4 + 5/4*x^3 + 5/8*x^2 - 5/16*x - 1/32
x^6 - 3*x^5 + 15/8*x^4 + 5/4*x^3 - 15/16*x^2 - 3/16*x + 17/256
x^7 - 7/2*x^6 + 21/8*x^5 + 35/16*x^4 - 35/16*x^3 - 21/32*x^2 + 119/256*x
+ 17/512
x^8 - 4*x^7 + 7/2*x^6 + 7/2*x^5 - 35/8*x^4 - 7/4*x^3 + 119/64*x^2 +
17/64*x - 31/256
x^9 - 9/2*x^8 + 9/2*x^7 + 21/4*x^6 - 63/8*x^5 - 63/16*x^4 + 357/64*x^3 +
153/128*x^2 - 279/256*x - 31/512
x^10 - 5*x^9 + 45/8*x^8 + 15/2*x^7 - 105/8*x^6 - 63/8*x^5 + 1785/128*x^4
+ 255/64*x^3 - 1395/256*x^2 - 155/256*x + 691/2048
1
x - 1/2
x^2 - x + 1/8
x^3 - 3/2*x^2 + 3/8*x + 1/16
x^4 - 2*x^3 + 3/4*x^2 + 1/4*x - 1/16
x^5 - 5/2*x^4 + 5/4*x^3 + 5/8*x^2 - 5/16*x - 1/32
x^6 - 3*x^5 + 15/8*x^4 + 5/4*x^3 - 15/16*x^2 - 3/16*x + 17/256
x^7 - 7/2*x^6 + 21/8*x^5 + 35/16*x^4 - 35/16*x^3 - 21/32*x^2 + 119/256*x + 17/512
x^8 - 4*x^7 + 7/2*x^6 + 7/2*x^5 - 35/8*x^4 - 7/4*x^3 + 119/64*x^2 + 17/64*x - 31/256
x^9 - 9/2*x^8 + 9/2*x^7 + 21/4*x^6 - 63/8*x^5 - 63/16*x^4 + 357/64*x^3 + 153/128*x^2 - 279/256*x - 31/512
x^10 - 5*x^9 + 45/8*x^8 + 15/2*x^7 - 105/8*x^6 - 63/8*x^5 + 1785/128*x^4 + 255/64*x^3 - 1395/256*x^2 - 155/256*x + 691/2048
#The next few plots show a symmetric distribution, where the asymptotic works, #and an asymmetric distribution, where it doesn't work plot(P([0,1/2,1],[1/4,1/2,1/4],20,x),xmin=-5,xmax=5,ymin=-10^4,ymax=10^4) 
       
plot(P([0,1/2,1],[1/4,1/2,1/4],40,x),xmin=-5,xmax=5,ymin=-10^18,ymax=10^18) 
       
plot(P([0,2/5,1],[1/4,1/2,1/4],40,x),xmin=-5,xmax=8,ymin=-10^22,ymax=10^22) 
       
#Creates a point plot of the roots of the (scaled) Strodt polynomial #Note: If the polynomial were not squarefree, this would loop forever from sage.rings.polynomial.complex_roots import complex_roots def PlotScaledRootsP(p,w,n): L=len(complex_roots(P(p,w,n,n*x),skip_squarefree=True)) l=[] for j in range(L): l.append((real(complex_roots(P(p,w,n,n*x),skip_squarefree=True)[j][0]).center(), imaginary(complex_roots(P(p,w,n,n*x),skip_squarefree=True)[j][0]).center())) return point(l,rgbcolor=hue(1),size=30) 
       
#Showing the roots of Strodt Polynomials for symmetric and asymmetric distributions PlotScaledRootsP([0,1/2,1],[1/4,1/2,1/4],20) 
       
PlotScaledRootsP([0,2/5,1],[1/4,1/2,1/4],20)