353 - Topic 06 - Distributions -- Sage

%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())

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Mean|\phantom{\verb!x!}\verb|=| \frac{1}{2} \, n + \frac{1}{2}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Variance|\phantom{\verb!x!}\verb|=| \frac{1}{12} \, {\left(n + 1\right)} {\left(n - 1\right)}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Skewness|\phantom{\verb!x!}\verb|=| 0$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Kurtosis|\phantom{\verb!x!}\verb|=| -\frac{6 \, {\left(n^{2} + 1\right)}}{5 \, {\left(n + 1\right)} {\left(n - 1\right)}}$

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

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Mean|\phantom{\verb!x!}\verb|=| \frac{a^{2} - b^{2}}{2 \, {\left(a - b\right)}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Variance|\phantom{\verb!x!}\verb|=| \frac{1}{12} \, {\left(a - b\right)}^{2}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Skewness|\phantom{\verb!x!}\verb|=| 0$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Kurtosis|\phantom{\verb!x!}\verb|=| \frac{9 \, {\left(a^{5} - 5 \, a^{4} b + 10 \, a^{3} b^{2} - 10 \, a^{2} b^{3} + 5 \, a b^{4} - b^{5}\right)} {\left(a - b\right)}}{5 \, {\left(a^{3} - 3 \, a^{2} b + 3 \, a b^{2} - b^{3}\right)}^{2}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Several|\phantom{\verb!x!}\verb|Examples|$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Using|\phantom{\verb!x!}\verb|[| 0 \verb|,| 2 \verb|]:|$

$\newcommand{\Bold}[1]{\mathbf{#1}}\phantom{\verb!xxxx!}\verb|mean|\phantom{\verb!x!}\verb|=| 1$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|variance|\phantom{\verb!x!}\verb|=| \frac{1}{3}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|skewness|\phantom{\verb!x!}\verb|=| 0$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|kurtosis|\phantom{\verb!x!}\verb|=| \frac{9}{5}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Using|\phantom{\verb!x!}\verb|[| 0 \verb|,| 3 \verb|]:|$

$\newcommand{\Bold}[1]{\mathbf{#1}}\phantom{\verb!xxxx!}\verb|mean|\phantom{\verb!x!}\verb|=| \frac{3}{2}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|variance|\phantom{\verb!x!}\verb|=| \frac{3}{4}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|skewness|\phantom{\verb!x!}\verb|=| 0$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|kurtosis|\phantom{\verb!x!}\verb|=| \frac{9}{5}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Using|\phantom{\verb!x!}\verb|[| 0 \verb|,| 4 \verb|]:|$

$\newcommand{\Bold}[1]{\mathbf{#1}}\phantom{\verb!xxxx!}\verb|mean|\phantom{\verb!x!}\verb|=| 2$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|variance|\phantom{\verb!x!}\verb|=| \frac{4}{3}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|skewness|\phantom{\verb!x!}\verb|=| 0$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|kurtosis|\phantom{\verb!x!}\verb|=| \frac{9}{5}$

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

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Mean|\phantom{\verb!x!}\verb|=| n p$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Variance|\phantom{\verb!x!}\verb|=| -n {\left(p - 1\right)} p$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Skewness|\phantom{\verb!x!}\verb|=| \frac{n {\left(2 \, p - 1\right)} {\left(p - 1\right)} p}{\left(-n {\left(p - 1\right)} p\right)^{\frac{3}{2}}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Kurtosis|\phantom{\verb!x!}\verb|=| -\frac{6 \, p^{2} - 6 \, p + 1}{n {\left(p - 1\right)} p} +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()

$\newcommand{\Bold}[1]{\mathbf{#1}}{\left(n - 1\right)} n p^{2} + n p$

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

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Mean|\phantom{\verb!x!}\verb|=| \frac{1}{p}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Variance|\phantom{\verb!x!}\verb|=| -\frac{p - 1}{p^{2}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Skewness|\phantom{\verb!x!}\verb|=| \frac{p^{2} - 3 \, p + 2}{p^{3} \left(-\frac{p - 1}{p^{2}}\right)^{\frac{3}{2}}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Kurtosis|\phantom{\verb!x!}\verb|=| -\frac{p^{2} - 6 \, p + 6}{p - 1} +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('Density Function: $f(x) =%s$'%str(latex(f(x)))) html('over the space $R = %s$'%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('Density Function: $f(x) =%s$'%str(latex(f(x)))) html('over the space $R = %s$'%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')

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Mean|\phantom{\verb!x!}\verb|=| \mu$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Variance|\phantom{\verb!x!}\verb|=| \mu$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Skewness|\phantom{\verb!x!}\verb|=| \frac{1}{\sqrt{\mu}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|Kurtosis|\phantom{\verb!x!}\verb|=| \frac{1}{\mu} +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.

353 - Topic 06 - Distributions

2739 days ago by Professor353

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

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

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

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

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

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

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

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

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