#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)
|
|
|
|