353 - Topic 06 - Distributions

305 days ago by Professor353

John Travis

Mississippi College

MATH 353 - Introduction to Mathematical Probability and Statistics

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

Uniform, Bernoulli, Binomial, Poisson, Negative Binomial, Hypergeometric Distributions

%hide %auto x = var('x') 
       
# Uniform distribution over 1 .. n pretty_print("Discrete Uniform Distribution over the set 1, 2, ..., n") var('x') @interact def _(n=slider(2,10,1,2)): np1 = n+1 R = range(1,np1) f(x) = 1/n pretty_print(html('Density Function: $f(x) =%s$'%str(latex(f(x)))+' over the space $R = %s$'%str(R))) points((k,f(x=k)) for k in R).show() for k in R: pretty_print(html('$f(%s'%k+') = %s'%f(x=k)+' \\approx %s$'%f(x=k).n(digits=5))) 
       

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

reset() var('x,n') f = 1/n mu = sum(x*f,x,1,n).factor() pretty_print('Mean = ',mu) mu = (1+n)/2 v = sum((x-mu)^2*f, x, 1, n) pretty_print('Variance = ',v.factor()) stand = sqrt(v) pretty_print('Skewness = ',(sum((x-mu)^3*f, x, 1, n)/stand^3)) kurt = sum((x-mu)^4*f, x, 1, n)/stand^4 pretty_print('Kurtosis = ',(kurt-3).factor()) 
       

                                
                            

                                
# Continous uniform distribution reset() var('x,a,b') assume(a<b) f = 1/(b-a) mu = integrate(x*f,x,a,b) pretty_print('Mean = ',mu) v = integrate((x-mu)^2*f, x, a, b) pretty_print('Variance = ',v.factor()) stand = sqrt(v) sk = (integrate((x-mu)^3*f, x, a, b)/stand^3) pretty_print('Skewness = ',sk) kurt = (integrate((x-mu)^4*f, x, a, b)/stand^4) pretty_print('Kurtosis = ',kurt) pretty_print('Several Examples') a1=0 for b1 in range(2,5): pretty_print('Using [',a1,',',b1,']:') pretty_print(' mean = ',mu(a=a1,b=b1)) pretty_print('variance = ',v(a=a1,b=b1)) pretty_print('skewness = ',sk(a=a1,b=b1)) pretty_print('kurtosis = ',kurt(a=a1,b=b1)) 
       

                                
                            

                                
histogram?? 
       
# Binomial distribution over 0 .. n # Probability of success on one independent trial = p must also be given var('x') @interact def _(n=slider(3,50,1,3),p=slider(1/20,19/20,1/20,1/2)): np1 = n+1 R = range(np1) f(x) = factorial(n)/(factorial(x)*factorial(n-x))*p^x*(1-p)^(n-x) pretty_print(html('Density Function: $f(x) =%s$'%str(latex(f(x))))) pretty_print(html('over the space $R = %s$'%str(R))) G = points((k,f(x=k)) for k in R) G.show() R = [k for k in R] probs = [f(x=k) for k in R] # histogram( R, weights = probs, align="mid", linewidth=2, edgecolor="blue", color="yellow").show() # H = histogram([k for k in R], weights=probs,align='mid') # H.show() for k in R: pretty_print(html('$f(%s'%k+') = %s'%latex(f(x=k))+' \\approx %s$'%f(x=k).n(digits=5))) 
       

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

# Binomial calculator @interact def _(p=input_box(0.3,width=15),n=input_box(10,width=15)): R = range(n+1) f(x) = binomial(n,x)*p^x*(1-p)^(n-x) acc = 0 for k in R: prob = f(x=k) acc = acc+prob pretty_print('f(%s) = '%k,' %.8f'%prob,' and F(%s) = '%k,' %.8f'%acc) 
       

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

var('x,n,p') assume(x,'integer') f(x) = binomial(n,x)*p^x*(1-p)^(n-x) mu = sum(x*f,x,0,n) M2 = sum(x^2*f,x,0,n) M3 = sum(x^3*f,x,0,n) M4 = sum(x^4*f,x,0,n) pretty_print('Mean = ',mu) v = (M2-mu^2).factor() pretty_print('Variance = ',v) stand = sqrt(v) sk = ((M3 - 3*M2*mu + 2*mu^3)).factor()/stand^3 pretty_print('Skewness = ',sk) kurt = (M4 - 4*M3*mu + 6*M2*mu^2 -3*mu^4).factor()/stand^4 pretty_print('Kurtosis = ',(kurt-3).factor(),'+3') 
       

                                
                            

                                
var('t') M = sum(e^(t*x)*f(n=n,p=p,k=x), x, 0, n).factor() M 
       
(p*e^t - p + 1)^n
(p*e^t - p + 1)^n
M.derivative(t,2)(t=0).show() 
       

                                
                            

                                
# Geometric distribution over 0 .. n # Probability of success on one independent trial = p must also be given var('x') # n = 50 by default. actually should be infinite @interact def _(p=input_box(0.1,label='p = '),n=[25,50,75,100,200]): np1 = n+1 R = range(1,np1) f(x) = (1-p)^(x-1)*p F(x) = 1 - (1-p)^x pretty_print(html('Density Function: $f(x) =%s$'%str(latex(f(x)))+' over the space $R = %s$'%str(R))) points((k,f(x=k)) for k in R).show(title="Probability Function") points((k,F(x=k)) for k in R).show(title="Distribution Function") if (n<51): for k in R: pretty_print(html('$f(%s'%k+') = %s'%latex(f(x=k))+' \\approx %s$'%f(x=k).n(digits=5))) 
       

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

var('x,n,p') assume(x,'integer') f(x) = p*(1-p)^(x-1) mu = sum(x*f,x,0,oo).full_simplify() M2 = sum(x^2*f,x,0,oo).full_simplify() M3 = sum(x^3*f,x,0,oo).full_simplify() M4 = sum(x^4*f,x,0,oo).full_simplify() pretty_print('Mean = ',mu) v = (M2-mu^2).factor().full_simplify() pretty_print('Variance = ',v) stand = sqrt(v) sk = (((M3 - 3*M2*mu + 2*mu^3))/stand^3).full_simplify() pretty_print('Skewness = ',sk) kurt = (M4 - 4*M3*mu + 6*M2*mu^2 -3*mu^4).factor()/stand^4 pretty_print('Kurtosis = ',(kurt-3).factor(),'+3') 
       

                                
                            

                                
# Negative Binomial distribution over 0 .. n # Probability of success on one independent trial = p must also be given var('x,r') n = 100 # actually should be infinite @interact def _(p=slider(1/24,23/24,1/24,1/2),r=slider(1,10,1,2)): np1 = n+1 R = range(r,np1) f(x) = (factorial(x-1)/(factorial(r-1)*factorial(x-r)))*(1-p)^(x-r)*p^r html('<font size=+2>Density Function: $f(x) =%s$</font>'%str(latex(f(x)))) html('<font size=+2>over the space $R = %s$</font>'%str(R)) points((k,f(x=k)) for k in R).show() for k in R: html('$f(%s'%k+') = %s'%latex(f(x=k))+' \\approx %s$'%f(x=k).n(digits=5)) 
       

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

# Negative Binomial calculator @interact def _(p=input_box(0.3,width=15),r=slider(1,10,1,2)): n = 4*(floor(r/p)+1) np1 = n+1 R = range(r,np1) f(x) = (factorial(x-1)/(factorial(r-1)*factorial(x-r)))*(1-p)^(x-r)*p^r acc = 0 for k in R: prob = f(x=k) acc = acc+prob pretty_print('f(%s) = '%k,' %.8f'%prob,' and F(%s) = '%k,' %.8f'%acc) 
       

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

# Negative Binomial var('x,n,p,r,alpha') assume(x,'integer') assume(alpha,'integer') assume(alpha>2) assume(p>0) assume(p<1) @interact def _(r=[2,5,10,15,alpha]): f(x) = binomial(x-1,r-1)*p^r*(1-p)^(x-r) mu = sum(x*f,x,r,oo).full_simplify() M2 = sum(x^2*f,x,r,oo).full_simplify() M3 = sum(x^3*f,x,r,oo).full_simplify() M4 = sum(x^4*f,x,r,oo).full_simplify() print M3 print M4 pretty_print('Mean = ',mu) v = (M2-mu^2).full_simplify() pretty_print('Variance = ',v) stand = sqrt(v) # sk = (((M3 - 3*M2*mu + 2*mu^3))/stand^3).full_simplify() sk = (((M3 - 3*M2*mu + 2*mu^3))).full_simplify() pretty_print('Skewness = ',sk) kurt = ((M4 - 4*M3*mu + 6*M2*mu^2 -3*mu^4)/v^2).full_simplify() pretty_print('Kurtosis = ',(kurt-3).factor(),'+3') 
       

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

# Poisson distribution over 0 .. inf # Average number of successes must be given var('x') n = 40 # actually should be infinite @interact def _(mu=slider(1/4,20,1/4,1,label='$\mu$')): np1 = n+1 R = range(np1) f(x) = mu^x*exp(-mu)/factorial(x) html('<font size=+2>Density Function: $f(x) =%s$</font>'%str(latex(f(x)))) html('<font size=+2>over the space $R = %s$</font>'%str(R)) points((k,f(x=k)) for k in R).show() for k in R: html('$f(%s'%k+') = %s'%latex(f(x=k))+' \\approx %s$'%f(x=k).n(digits=5)) 
       

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

var('x,mu') assume(x,'integer') f(x) =e^(-mu)*mu^x/factorial(x) mu = sum(x*f,x,0,oo).factor() M2 = sum(x^2*f,x,0,oo).factor() M3 = sum(x^3*f,x,0,oo).factor() M4 = sum(x^4*f,x,0,oo).factor() pretty_print('Mean = ',mu) v = (M2-mu^2).factor() pretty_print('Variance = ',v) stand = sqrt(v) sk = ((M3 - 3*M2*mu + 2*mu^3)).factor()/stand^3 pretty_print('Skewness = ',sk) kurt = (M4 - 4*M3*mu + 6*M2*mu^2 -3*mu^4).factor()/stand^4 pretty_print('Kurtosis = ',(kurt-3).factor(),'+3') 
       

                                
                            

                                
# Hypergeometric distribution over 0 .. N # Size of classes N1 and N2 must be given as well as subset size r var('x') @interact def _(N1=slider(1,40,1,10,label='$N_1$'), N2=slider(1,40,1,10,label='$N_2$'), r=slider(1,40,1,10,label='$r$')): N = N1 + N1 R = range(r+1) if (r>N1)|(r>N2): pretty_print('When r is bigger than N1 or N2, special consideration must be made') else: f(x) = binomial(N1,x)*binomial(N2,r-x)/binomial(N,r) pretty_print(html('Density Function: $f(x) =%s$'%str(latex(f(x))))) pretty_print(html('over the space $R = %s$'%str(R))) points((k,f(x=k)) for k in R).show() for k in R: print (html('$f(%s'%k+') = %s'%latex(f(x=k))+' \\approx %s$'%f(x=k).n(digits=5))) 
       

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

reset() var('x,r,N1,N2') f(x,N1=10, N2=10, r=5) = binomial(N1,x)*binomial(N2,r-x)/binomial(N,r) # mu = sum(x*f,x,0,r).factor() # pretty_print('Mean = ',mu) 
       
Traceback (click to the left of this block for traceback)
...
ValueError: The name "N1=_sage_const_10" is not a valid Python
identifier.
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_96.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cmVzZXQoKQp2YXIoJ3gscixOMSxOMicpCgpmKHgsTjE9MTAsIE4yPTEwLCByPTUpID0gYmlub21pYWwoTjEseCkqYmlub21pYWwoTjIsci14KS9iaW5vbWlhbChOLHIpCiMgbXUgPSBzdW0oeCpmLHgsMCxyKS5mYWN0b3IoKQojIHByZXR0eV9wcmludCgnTWVhbiA9ICcsbXUp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmppyhVmS/___code___.py", line 6, in <module>
    __tmp__=var("x,N1=_sage_const_10,N2=_sage_const_10,r=_sage_const_5"); f = symbolic_expression(binomial(N1,x)*binomial(N2,r-x)/binomial(N,r)).function(x,N1=_sage_const_10,N2=_sage_const_10,r=_sage_const_5)
  File "sage/calculus/var.pyx", line 140, in sage.calculus.var.var (build/cythonized/sage/calculus/var.c:977)
  File "sage/symbolic/ring.pyx", line 613, in sage.symbolic.ring.SymbolicRing.var (build/cythonized/sage/symbolic/ring.cpp:8883)
  File "sage/symbolic/ring.pyx", line 666, in sage.symbolic.ring.SymbolicRing.var (build/cythonized/sage/symbolic/ring.cpp:8557)
ValueError: The name "N1=_sage_const_10" is not a valid Python identifier.