#Strodt Generating Function generalized to continuous probability distribution on a real interval I
#Input f(x) is the continuous density function on I
#Inputs a and b are the lower and upper limits of the interval I (a could be -oo, b could be oo)
#Note: Make sure f is given as a function of 'x' which is cts. on [a,b]
#Note: I'm not sure how well Sage does symbolic integration
x,t=var('x,t')
def CtsStGF(f,a,b,x,t):
Q=integral(f*exp(x*t),(x,a,b))
return exp(x*t)/Q
|
|
#Pulls the Strodt polynomials out from the generating function
#Again f must be a function of 'x' cts. on [a,b].
def P(f,a,b,n,x):
g=taylor(CtsStGF(f,a,b,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 Gaussian distribution
#Note: These are (scaled) Hermite polynomials
f=exp(-x^2)/sqrt(pi)
for n in range(10):
P(f,-oo,oo,n,x)
1 x x^2 - 1/2 x^3 - 3/2*x x^4 - 3*x^2 + 3/4 x^5 - 5*x^3 + 15/4*x x^6 - 15/2*x^4 + 45/4*x^2 - 15/8 x^7 - 21/2*x^5 + 105/4*x^3 - 105/8*x x^8 - 14*x^6 + 105/2*x^4 - 105/2*x^2 + 105/16 x^9 - 18*x^7 + 189/2*x^5 - 315/2*x^3 + 945/16*x 1 x x^2 - 1/2 x^3 - 3/2*x x^4 - 3*x^2 + 3/4 x^5 - 5*x^3 + 15/4*x x^6 - 15/2*x^4 + 45/4*x^2 - 15/8 x^7 - 21/2*x^5 + 105/4*x^3 - 105/8*x x^8 - 14*x^6 + 105/2*x^4 - 105/2*x^2 + 105/16 x^9 - 18*x^7 + 189/2*x^5 - 315/2*x^3 + 945/16*x |
#Plotting Strodt Polynomials, in the Gaussian case
f=exp(-x^2)/sqrt(pi)
plot(P(f,-oo,oo,20,x),xmin=-5,xmax=5,ymin=-10^7,ymax=10^7)
|
|
f=exp(-x^2)/sqrt(pi)
plot(P(f,-oo,oo,40,x),xmin=-10,xmax=10,ymin=-10^20,ymax=10^20)
|
|
#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(f,a,b,n):
L=len(complex_roots(P(f,a,b,n,n*x),skip_squarefree=True))
l=[]
for j in range(L):
l.append((real(complex_roots(P(f,a,b,n,n*x),skip_squarefree=True)[j][0]).center(), imaginary(complex_roots(P(f,a,b,n,n*x),skip_squarefree=True)[j][0]).center()))
return point(l,rgbcolor=hue(1),size=30)
|
|
#Showing the roots of Strodt Polynomials for gaussian distribution (Hermite)
f=exp(-x^2)/sqrt(pi)
PlotScaledRootsP(f,-oo,oo,20)
|
|
|
|