SM-5-Exercise4

4137 days ago by MathFest

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