353 - Topic 10 - Normal Curve

342 days ago by Professor353

Derivation of Details related to the normal distribution

John Travis

Mississippi College

MATH 353 - Introduction to Mathematical Probability and Statistics

Textbook:  Tanis and Hogg, A Brief Course in Mathematical Statistics

Using the symbolic capabilities of Sage to determine all of the interesting characteristics of the bell curve.

%auto reset() var('x,mu,sigma') f(x) = e^(-((x-mu)/sigma)^2/2)/(sigma*sqrt(2*pi)) f.show() 
       

                                
                            

                                
assume(sigma>0) mean = integrate(x*f,x,-oo,oo).simplify() html('By integrating $x f(x)$ over $(-\infty,\infty)$, we get the mean is $ %s$'%str(latex(mean))) 
       
By integrating  over , we get the mean is 
                                
                            
By integrating  over , we get the mean is 
                                
vari = integrate((x-mu)^2*f,x,-oo,oo).simplify() html('By integrating $x^2 f(x)$ over $(-\infty,\infty)$, we ultimately get variance is $ %s$'%str(latex(vari))) 
       
By integrating  over , we ultimately get variance is 
                                
                            
By integrating  over , we ultimately get variance is 
                                
skew = integrate((x-mu)^3*f/sigma^3,x,-oo,oo).simplify() html('By integrating $x^2 f(x)$ over $(-\infty,\infty)$, we ultimately get skewness is $ %s$'%str(latex(skew))) 
       
By integrating  over , we ultimately get skewness is 
                                
                            
By integrating  over , we ultimately get skewness is 
                                
kurt = integrate((x-mu)^4*f/sigma^4,x,-oo,oo).simplify() html('By integrating $x^2 f(x)$ over $(-\infty,\infty)$, we ultimately get kurtosis is $ %s$'%str(latex(kurt))) 
       
By integrating  over , we ultimately get kurtosis is 
                                
                            
By integrating  over , we ultimately get kurtosis is 
                                
# to find the maximum of f(x) fp = diff(f,x) soln = solve(fp==0,x) # can't make soln = mu work right here. Going ahead and just plugging in. Grrr. big = f(mu) html('Taking the derivative equal to zero gives a critical value of $%s$'%str(latex(soln[0]))) html('By glancing at the function, then there is a max at $%s$'%str(latex(big))) 
       
Taking the derivative equal to zero gives a critical value of 
By glancing at the function, then there is a max at 
                                
                            
Taking the derivative equal to zero gives a critical value of 
By glancing at the function, then there is a max at 
                                
# to find possible points of inflection of f(x) fpp = diff(f,x,2) soln = solve(fpp==0,x) # again, can't plug these solns in yet. turn = f(mu+sigma).simplify() html('Taking the second derivative equal to zero gives possible inflection points at $%s$'%str(latex(soln[0]))+' and $%s$'%str(latex(soln[1]))) html('Plugging into the function gives inflection points $%s$'%str(latex((mu+sigma,turn)))+' and $%s$'%str(latex((mu-sigma,turn)))) 
       
Taking the second derivative equal to zero gives possible inflection points at  and 
Plugging into the function gives inflection points  and 
                                
                            
Taking the second derivative equal to zero gives possible inflection points at  and 
Plugging into the function gives inflection points  and 
                                
html('A nice graph of the function f(x) with maximum at the red point') html('and green inflection points.') var('x,mu,sigma') f(x) = e^(-((x-mu)/sigma)^2/2)/(sigma*sqrt(2*pi)) f.show() @interact def _(m=slider(-10,10,1,0,label='$\mu$'),s=slider(1/5,5,1/10,1,label='$\sigma$')): titletext = "Normal Curve with mean "+str(m)+" and standard deviation "+str(s) G = plot(f(mu=m,sigma=s),(x,m-5*s,m+5*s)) G += point((0,1),size=1)+point((12,0),size=1)+point((-12,0),size=1) G += point((m,f(x=m,mu=m,sigma=s)),color='red',size=20) G += point((m+s,f(x=m+s,mu=m,sigma=s)),color='green',size=20) G += point((m-s,f(x=m-s,mu=m,sigma=s)),color='green',size=20) show(G,figsize=(5,3),title=titletext,ymin=0,ymax=1,xmin=-15,xmax=15) 
       
A nice graph of the function f(x) with maximum at the red point and green inflection points.
$\mu$ 
$\sigma$ 

Click to the left again to hide and once more to show the dynamic interactive window

# Let's redo this using Moment Generating Functions var('t') M = integrate(e^(t*x)*f,x,-oo,oo).simplify() M.show() 
       

                                
                            

                                
M(t=0) 
       
1
1
MP = derivative(M,t) MP(t=0).show() 
       

                                
                            

                                
MPP = derivative(M,t,2) MPP(t=0).show() 
       

                                
                            

                                
MPPP = derivative(M,t,3) MPPP(t=0).show() skewness = ((MPPP(t=0)-3*mu*MPP(t=0)+2*mu^3)/sigma^3).simplify() skewness.show() # Note this is really just zero! 
       

                                
                            

                                
MPPPP = derivative(M,t,4) MPPPP(t=0).show() kurtosis = ((MPPPP(t=0)-4*mu*MPPP(t=0)+6*mu^2*MPP(t=0)-3*mu^4)/sigma^4).simplify() kurtosis.factor() 
       
3
3
@interact(layout=dict(top=[['a', 'b']],bottom=[['mu','sigma']])) def _(a=input_box(-2,width=10,label='a = '),b=input_box(2,width=10,label='b = '),mu=input_box(0,width=8,label='$\mu = $'),sigma=input_box(1,width=8,label='$\sigma = $')): f = e^(-((x-mu)/sigma)^2/2)/(sigma*sqrt(2*pi)) P = integral_numerical(f,a,b)[0] print "P("+str(a)+" < X < "+str(b)+") ~= "+str(P) 
       
a =   b =  
$\mu = $  $\sigma = $ 

Click to the left again to hide and once more to show the dynamic interactive window