︠7ad8647f-da30-4412-92cc-78fccf40db91︠
# Example 5.1.5
x, y,t,r = var('x y t r')
def cf(r,t): return 2*sin(t)*cos(t)*r
P = parametric_plot3d((r*sin(t),r*cos(t),sin(t)*cos(t)),(r,0,.7),(t,-pi,pi),color=(cf,colormaps.ocean),opacity=0.7,frame=false)
PL=plot3d(0, (-1,1), (-1,1),opacity=0.2,mash=false,frame=false)
from sage.plot.plot3d.plot3d import axes
T = axes(1, .5);
(P+PL+T).show()
︡17d5be20-d3f1-48db-867d-c03caa107627︡{"file":{"filename":"baa4b3fb-c53b-487a-bd8b-9a8d45341185.sage3d","uuid":"baa4b3fb-c53b-487a-bd8b-9a8d45341185"}}︡{"done":true}︡
︠51780d2b-c5b7-4746-8d30-95a627d8c1c9︠
# Example 5.1.10
x,y = var('x y')
z = x*exp(x*y);
z_x = diff(z,x);show('$z_x = $',z_x)
z_y = diff(z,y);show('$z_y = $',z_y)
z_xx = diff(z_x,x);show('$z_{xx} = $',z_xx)
z_yy = diff(z_y,y);show('$z_{yy} = $',z_yy)
z_xy = diff(z_y,x);show('$z_{xy} = $',z_xy)
z_yx = diff(z_x,y);show('$z_{yx} = $',z_yx)
︡f2ce1488-960e-4bea-99db-b2ba9ab5fbd1︡{"html":"
$z_x = $ $\\displaystyle x y e^{\\left(x y\\right)} + e^{\\left(x y\\right)}$
"}︡{"html":"$z_y = $ $\\displaystyle x^{2} e^{\\left(x y\\right)}$
"}︡{"html":"$z_{xx} = $ $\\displaystyle x y^{2} e^{\\left(x y\\right)} + 2 \\, y e^{\\left(x y\\right)}$
"}︡{"html":"$z_{yy} = $ $\\displaystyle x^{3} e^{\\left(x y\\right)}$
"}︡{"html":"$z_{xy} = $ $\\displaystyle x^{2} y e^{\\left(x y\\right)} + 2 \\, x e^{\\left(x y\\right)}$
"}︡{"html":"$z_{yx} = $ $\\displaystyle x^{2} y e^{\\left(x y\\right)} + 2 \\, x e^{\\left(x y\\right)}$
"}︡{"done":true}︡
︠5d34470c-38c5-4dff-a999-cf76efed58d5︠
# Example 5.1.11
x,y = var('x y')
z = ln(sqrt(x^2+y^2));
z_xx = diff(z,x,2);show('$z_{xx} = $',z_xx)
z_yy = diff(z,y,2);show('$z_{yy} = $',z_yy)
show('$z_{xx} + z_{yy} = $', (z_xx+z_yy).simplify_full())
︡1f241663-0081-43f8-8be8-c01e27fe6f55︡{"html":"$z_{xx} = $ $\\displaystyle -\\frac{2 \\, x^{2}}{{\\left(x^{2} + y^{2}\\right)}^{2}} + \\frac{1}{x^{2} + y^{2}}$
"}︡{"html":"$z_{yy} = $ $\\displaystyle -\\frac{2 \\, y^{2}}{{\\left(x^{2} + y^{2}\\right)}^{2}} + \\frac{1}{x^{2} + y^{2}}$
"}︡{"html":"$z_{xx} + z_{yy} = $ $\\displaystyle 0$
"}︡{"done":true}︡
︠a06fd64a-5a15-419a-a021-e92ce02a5585︠
# Example estimate 1.04^{1.98}
show("$1.04^{1.98}=$",exp(ln(1.04)*1.98))
z(x,y) = x^y;
x0=1;y0=2;
dx = 0.04; dy=-0.02;
show("estimate = ", z(x0,y0)+ diff(z,x)(x0,y0)*dx+diff(z,y)(x0,y0)*dy)
︡d234c7e5-e740-4862-bc1b-758442373df5︡{"html":"$1.04^{1.98}=$ $\\displaystyle 1.08075191020342$
"}︡{"html":"estimate = $\\displaystyle 1.08000000000000$
"}︡{"done":true}︡
︠5650e9b6-c8ea-493c-83c9-c2ea84c4452a︠
# Section 5.1.3.3
x, y = var('x y')
u=function('u')(x,y)
v=function('v')(x,y)
z=function('z')(u,v)
show("$z_{xy}$=",diff(diff(z,x),y).expand())
︡d469e73d-868c-444e-96b0-b47917fdfff1︡{"html":"$z_{xy}$= $\\displaystyle \\frac{\\partial}{\\partial x}u\\left(x, y\\right) \\frac{\\partial}{\\partial y}u\\left(x, y\\right) \\mathrm{D}_{0, 0}\\left(z\\right)\\left(u\\left(x, y\\right), v\\left(x, y\\right)\\right) + \\frac{\\partial}{\\partial y}u\\left(x, y\\right) \\frac{\\partial}{\\partial x}v\\left(x, y\\right) \\mathrm{D}_{0, 1}\\left(z\\right)\\left(u\\left(x, y\\right), v\\left(x, y\\right)\\right) + \\frac{\\partial}{\\partial x}u\\left(x, y\\right) \\frac{\\partial}{\\partial y}v\\left(x, y\\right) \\mathrm{D}_{0, 1}\\left(z\\right)\\left(u\\left(x, y\\right), v\\left(x, y\\right)\\right) + \\frac{\\partial}{\\partial x}v\\left(x, y\\right) \\frac{\\partial}{\\partial y}v\\left(x, y\\right) \\mathrm{D}_{1, 1}\\left(z\\right)\\left(u\\left(x, y\\right), v\\left(x, y\\right)\\right) + \\frac{\\partial^{2}}{\\partial x\\partial y}u\\left(x, y\\right) \\mathrm{D}_{0}\\left(z\\right)\\left(u\\left(x, y\\right), v\\left(x, y\\right)\\right) + \\frac{\\partial^{2}}{\\partial x\\partial y}v\\left(x, y\\right) \\mathrm{D}_{1}\\left(z\\right)\\left(u\\left(x, y\\right), v\\left(x, y\\right)\\right)$
"}︡{"done":true}︡
︠a1949be9-b198-45a1-9251-060056d481d6︠
# Example 5.1.5
x, y = var('x y')
def cf(x,y): return abs(x*y)
P = plot3d(x*y,(x,-1,1),(y,-1,1),color=(cf,colormaps.rainbow),opacity=0.7,frame=false)
from sage.plot.plot3d.plot3d import axes
T = axes(1, .5,color='black');
(P+T).show()
︡92bf3a65-f2f0-453d-9920-a9836392e064︡
︠969744e3-60f8-4c71-b8f4-21e425c410d6︠
︡70b943a3-f0dc-4acb-b585-c71e2df9c0fe︡
︠3f94e6af-bf82-4b5e-9c06-5f11180747c1︠
# Example 5.1.5-1
def cf(r,t): return abs(3*r^2*cos(t)*sin(t))
T = Cylindrical('height', ['radius', 'azimuth'])
r, t, z = var('r theta z')
#P=plot3d(r^2, (r, 0, 1), (t, -pi, pi), color=(cf,colormaps.rainbow),opacity=0.7,frame=false)
PP = plot3d(r^2*cos(t)*sin(t), (r, 0, 1), (t, 0, 2*pi), transformation=T,color=(cf,colormaps.rainbow),axes=true, adaptive=True,frame=false,opacity=0.7,aspect_ratio=[1,1,1.5])
from sage.plot.plot3d.plot3d import axes
CD = axes(1, .5,color='black');
(PP+CD).show()
︡6adcb2ab-f3a0-41d9-9b64-5b599853649c︡{"file":{"filename":"d211b952-3259-47cb-bbc9-ff0cb71ac13d.sage3d","uuid":"d211b952-3259-47cb-bbc9-ff0cb71ac13d"}}︡{"done":true}
︠1ac6f426-e6be-4666-ba1a-1cacbe02f32d︠
def cf(r,theta): return abs(r^2)
T = Cylindrical('height', ['radius', 'azimuth'])
r, t, z = var('r theta z')
#P=plot3d(r^2, (r, 0, 1), (t, -pi, pi), color=(cf,colormaps.rainbow),opacity=0.7,frame=false)
PP = plot3d(r^2, (r, 0, 1), (t, 0, 2*pi), transformation=T,color=(cf,colormaps.rainbow),axes=true, frame=false,opacity=0.7,aspect_ratio=[1,1,1.5])
from sage.plot.plot3d.plot3d import axes
CD = axes(1, .5,color='black');
(PP+CD).show()
︡a177d8ce-8501-4a34-bb90-7cd2bd052b2f︡
︠fb31991a-8598-4105-a268-109cde961d86︠
def cf(r,theta): return abs(r^2)
T = Cylindrical('height', ['radius', 'azimuth'])
r, t, z = var('r theta z')
#P=plot3d(r^2, (r, 0, 1), (t, -pi, pi), color=(cf,colormaps.rainbow),opacity=0.7,frame=false)
PP = plot3d(2*r^2, (r, 0, 0.5), (t, 0, 2*pi), transformation=T,color=(cf,colormaps.rainbow),axes=true,
frame=false,opacity=0.7,aspect_ratio=1)
def cff(r,t): return cos(4*r)
PPP = plot3d((1-cos(4*r))/4, (r, 0, 1), (t, 0, 1.5*pi), transformation=T,
#color=(cff,colormaps.rainbow),mesh=true,
frame=false,opacity=0.7,aspect_ratio=1)
from sage.plot.plot3d.plot3d import axes
CD = axes(1.5, .5,color='black');
(PP+PPP+CD).show()
︡6af0f03b-1ea3-47d4-a8bf-d135028651bb︡{"file":{"filename":"52929655-2aab-4935-8c44-3f676e81bbc1.sage3d","uuid":"52929655-2aab-4935-8c44-3f676e81bbc1"}}︡{"done":true}︡{"done":true}︡
︠438597ec-baa1-49be-abd2-0a2124e27104︠
# Example 5.1.30
x,y,lam= var('x y lambda_');
F = 2*x^2+6*x*y+y^2
S0=solve([diff(F,x)==0,diff(F,y)==0],x,y)
show(S0)
show("Function value:",[F.subs(s) for s in S0])
G = x^2+2*y^2-3
L = F + lam*G
Lx = diff(L,x)
Ly = diff(L,y)
Llam = diff(L,lam)
S = solve([Lx==0,Ly==0,Llam==0],x,y,lam)
show(S)
show("Function value:",[F.subs(s) for s in S])
︡e3899ca9-23c3-4e0f-9fa4-711b6cfae93a︡{"html":"[[$\\displaystyle x = 0$, $\\displaystyle y = 0$]]
"}︡{"html":"Function value: [$\\displaystyle 0$]
"}︡{"html":"[[$\\displaystyle x = \\sqrt{2}$, $\\displaystyle y = \\frac{1}{2} \\, \\sqrt{2}$, $\\displaystyle \\lambda = \\left(-\\frac{7}{2}\\right)$], [$\\displaystyle x = -\\sqrt{2}$, $\\displaystyle y = -\\frac{1}{2} \\, \\sqrt{2}$, $\\displaystyle \\lambda = \\left(-\\frac{7}{2}\\right)$], [$\\displaystyle x = 1$, $\\displaystyle y = \\left(-1\\right)$, $\\displaystyle \\lambda = 1$], [$\\displaystyle x = \\left(-1\\right)$, $\\displaystyle y = 1$, $\\displaystyle \\lambda = 1$]]
"}︡{"html":"Function value: [$\\displaystyle \\frac{21}{2}$, $\\displaystyle \\frac{21}{2}$, $\\displaystyle -3$, $\\displaystyle -3$]
"}︡{"done":true}︡
︠24d420ba-fbdf-4668-9e0b-b1805bb36315︠
# Example on ppt page 17
x,y = var('x y')
I1 = integral(integral(y,y,0,sqrt(1-x^2)),x,-1,1)
I2 = integral(integral(y+y^2*sin(x),y,0,sqrt(1-x^2)),x,-1,1)
show(I1,"=",I2)
︡cd7d0b5b-4e6b-498c-8850-618d66e4dedd︡{"html":"$\\displaystyle \\frac{2}{3}$ = $\\displaystyle \\frac{2}{3}$
"}︡{"done":true}︡
︠cc7af456-9394-4195-aba0-246545a8d22d︠
# Example on ppt page 15
x, y = var('x y')
show(integral(integral(x+2*y,y,0,x^2),x,0,2))
show(integral(integral(x*y,x,y^2,y+2),y,-1,2))
assume(x>0);assume(y>0);
# show(integral(integral(sin(y)/y,y,x,sqrt(x)),x,0,1)) # cause error!
show(integral(integral(sin(y)/y,x,y^2,y),0,1))
#Example on ppt page 16
# show(integral(integral(cos(y)/y,y,x,pi/6),x,0,pi/6)) # cause error!
show(integral(integral(cos(y)/y,x,0,y),y,0,pi/6))
︡bb20c203-dd94-43e9-8b92-d0276ae99352︡{"html":"$\\displaystyle \\frac{52}{5}$
"}︡{"html":"$\\displaystyle \\frac{45}{8}$
"}︡{"html":"$\\displaystyle -\\sin\\left(1\\right) + 1$
"}︡{"html":"$\\displaystyle \\frac{1}{2}$
"}︡{"done":true}︡
︠257ea159-fddb-40f3-a5e2-2dc076d140bf︠
# Example on ppt page 21
r,t=var('r t')
show(integral(integral(r^3,r,0,2*cos(t)),t,-pi/2,pi/2))
︡a83fc108-8f1c-4256-ad38-d71be37eee78︡{"html":"$\\displaystyle \\frac{3}{2} \\, \\pi$
"}︡{"done":true}︡
︠8700790c-57fa-442a-ac41-7a3329d76960︠
# Example on ppt page 22
r,t,R=var('r t R')
assume(R>0)
#show(integral(r*sqrt(R^2-r^2),r,0,R*cos(t)).expand())
show(2*integral(1/3*R^3*(1-(sin(t))^3),t,0,pi/2))
︡c9268c12-4a2a-4bf7-901e-b7a3a4f4c151︡{"html":"$\\displaystyle \\frac{1}{9} \\, {\\left(3 \\, \\pi - 4\\right)} R^{3}$
"}︡{"done":true}︡
︠5515afd2-f44f-4c3f-9067-20a15c44125c︠
# Example on ppt page 23
r, t,a = var('r t a')
assume(a>0)
A = 4*integral(integral(r,r,0,sqrt(2*a^2*((cos(t))^2-(sin(t))^2))),t,0,pi/4)
show(A)`
︡0a447f9c-0ed9-4980-b3b1-2d3748fa0140︡{"html":"$\\displaystyle 2 \\, a^{2}$
"}︡{"done":true}︡
︠b70420ab-c2a0-4bb3-a7dd-f54bf59d6483︠
#################
### Chapter 6 ###
#################
︡edf08eb1-39c1-4eb4-8076-d2462de0ca59︡
︠fcc2730e-46f1-4fd1-b5d8-857bef8e5fde︠
# Example 6.2.1
P = binomial(100-3,4)*binomial(3,1)/binomial(100,5);
show(P,'=',float(P))
︡5362f37f-80f0-4ae1-894e-c10cf45b1b6d︡{"html":"$\\displaystyle \\frac{893}{6468}$ = $\\displaystyle 0.138064316636$
"}︡{"done":true}︡
︠f557910d-6511-404a-841d-f3b85d8332a1︠
# Example
P = binomial(50-3,2)*binomial(3,2)/binomial(50,4)+binomial(50-3,1)*binomial(3,3)/binomial(50,4)
show(float(P))
︡0e880a9b-6786-4d0d-884c-13245786bfe8︡{"html":"$\\displaystyle 0.0142857142857$
"}︡{"done":true}︡
︠5c22f7d8-c9d8-4d5b-a169-353931fb7c41︠
# Example 6.2.3
def Pc(n):
return factorial(n)*binomial(365,n)/365^n
for i in xrange(2,51):
print i, 'students; ', float(1-Pc(i))
︡789896cd-4f60-44e1-9126-5f51db0280de︡{"stdout":"2 students; 0.0027397260274\n3 students; 0.00820416588478\n4 students; 0.0163559124666\n5 students; 0.0271355736998\n6 students; 0.0404624836491\n7 students; 0.056235703096\n8 students; 0.0743352923517\n9 students; 0.0946238338892\n10 students; 0.116948177711\n11 students; 0.141141378322\n12 students; 0.167024788838\n13 students; 0.194410275232\n14 students; 0.223102512005\n15 students; 0.252901319764\n16 students; 0.283604005253\n17 students; 0.315007665297\n18 students; 0.346911417872\n19 students; 0.379118526032\n20 students; 0.411438383581\n21 students; 0.443688335165\n22 students; 0.475695307663\n23 students; 0.507297234324\n24 students; 0.538344257915\n25 students; 0.568699703969\n26 students; 0.598240820136\n27 students; 0.626859282263\n28 students; 0.654461472342\n29 students; 0.680968537478\n30 students; 0.706316242719\n31 students; 0.730454633729\n32 students; 0.75334752785\n33 students; 0.774971854176\n34 students; 0.79531686462\n35 students; 0.814383238875\n36 students; 0.83218210638\n37 students; 0.848734008216\n38 students; 0.864067821082\n39 students; 0.878219664367\n40 students; 0.891231809818\n41 students; 0.903151611482\n42 students; 0.914030471562\n43 students; 0.923922855656\n44 students; 0.932885368551\n45 students; 0.940975899466\n46 students; 0.948252843367\n47 students; 0.954774402833\n48 students; 0.960597972879\n49 students; 0.965779609323\n50 students; 0.970373579578\n"}︡{"done":true}︡
︠b06f685d-e7d6-4d17-8f32-2b5559c8db66︠
# Example on ppt page 40
show(1-(1-0.004)^100)
︡a0254bbc-8a60-483a-a3f1-6ec2fd1c732b︡{"html":"$\\displaystyle 0.330217428727354$
"}︡{"done":true}︡
︠2742e332-20d5-4ff8-96af-149853461f1d︠
# Example on ppt page 43
def b(k,n,p):
return binomial(n,k)*p^k*(1-p)^(n-k)
P = sum([b(i,10,12/60) for i in xrange(0,6)])
show(float(P))
# Example 6.3.8
P = 1-b(0,6,1/20)-b(1,6,1/20);
show(float(P))
︡89cab4a7-e51b-4fe6-87d3-9ab7eb27b64c︡{"html":"$\\displaystyle 0.9936306176$
"}︡{"html":"$\\displaystyle 0.032773828125$
"}︡{"done":true}︡
︠3b9b0acd-b177-4e6c-8cf5-1c7291b57f97︠
# Example on ppt page 47
PA = 0.0004
PB_A=0.94
PB_AC=0.04
show(PA*PB_A/(PA*PB_A+(1-PA)*PB_AC))
︡965339db-ec5f-46eb-b5c0-c18288d38408︡{"html":"$\\displaystyle 0.00931615460852329$
"}︡{"done":true}︡
︠666670b3-daae-4d85-93fb-caa0bc220b52︠
# Example on 6.5 ppt page 8
from scipy.stats import binom
binom_dist = binom(5,0.1)
# bar_chart([binom_dist.pmf(x) for x in range(6)])
for x in range(6):
print "P(X=%d) = %f"%(x,binom_dist.pmf(x))
︡14ade4de-e0c2-4a0d-9dd0-74cbad042c78︡{"stdout":"P(X=0) = 0.590490\nP(X=1) = 0.328050\nP(X=2) = 0.072900\nP(X=3) = 0.008100\nP(X=4) = 0.000450\nP(X=5) = 0.000010\n"}︡{"done":true}︡
︠a8a63eb3-a142-4390-8045-c8933bea81e3︠
# Example on 6.5 ppt page 9
from scipy.stats import binom
binom_dist = binom(8,0.3)
# bar_chart([binom_dist.pmf(x) for x in range(6)])
for x in range(9):
print "P(X=%d) = %f"%(x,binom_dist.pmf(x))
print "P(X>=2) = %f"%(1-binom_dist.pmf(0)-binom_dist.pmf(1)) # or use: (1-binom_dist.cdf(1))
︡9070ff71-ede2-456b-86fa-2c18f90c8012︡{"stdout":"P(X=0) = 0.057648\nP(X=1) = 0.197650\nP(X=2) = 0.296475\nP(X=3) = 0.254122\nP(X=4) = 0.136137\nP(X=5) = 0.046675\nP(X=6) = 0.010002\nP(X=7) = 0.001225\nP(X=8) = 0.000066\n"}︡{"stdout":"P(X>=2) = 0.744702\n"}︡{"done":true}︡
︠8633e150-f9f6-44f8-9ab3-499426969b6b︠
# Example 6.5.3
from scipy.stats import poisson
for k in xrange(0,20):
show("$P(X\leq %i)=%f$"%(k, poisson.cdf(k, mu=10)))
︡06ba65d0-c0b3-469f-ad41-c1ab81d85d5d︡{"html":"$P(X\\leq 0)=0.000045$
"}︡{"html":"$P(X\\leq 1)=0.000499$
"}︡{"html":"$P(X\\leq 2)=0.002769$
"}︡{"html":"$P(X\\leq 3)=0.010336$
"}︡{"html":"$P(X\\leq 4)=0.029253$
"}︡{"html":"$P(X\\leq 5)=0.067086$
"}︡{"html":"$P(X\\leq 6)=0.130141$
"}︡{"html":"$P(X\\leq 7)=0.220221$
"}︡{"html":"$P(X\\leq 8)=0.332820$
"}︡{"html":"$P(X\\leq 9)=0.457930$
"}︡{"html":"$P(X\\leq 10)=0.583040$
"}︡{"html":"$P(X\\leq 11)=0.696776$
"}︡{"html":"$P(X\\leq 12)=0.791556$
"}︡{"html":"$P(X\\leq 13)=0.864464$
"}︡{"html":"$P(X\\leq 14)=0.916542$
"}︡{"html":"$P(X\\leq 15)=0.951260$
"}︡{"html":"$P(X\\leq 16)=0.972958$
"}︡{"html":"$P(X\\leq 17)=0.985722$
"}︡{"html":"$P(X\\leq 18)=0.992813$
"}︡{"html":"$P(X\\leq 19)=0.996546$
"}︡{"done":true}︡
︠e8441736-3a79-4cf9-bf57-cdb708a215b9︠
# Example 2 on ppt page 13
from math import factorial,exp
from scipy.stats import poisson
mu=var("mu")
S = solve([mu/1==mu^2/2],mu)
show(S)
show(2.0^4/(factorial(4))*exp(-2.0))
show(poisson.pmf(k=4,mu=2))
︡e7464967-24be-4e2a-aa88-35a9518c7eb4︡{"html":"[$\\displaystyle \\mu = 0$, $\\displaystyle \\mu = 2$]
"}︡{"html":"$\\displaystyle 0.0902235221577418$
"}︡{"html":"$\\displaystyle 0.0902235221577$
"}︡{"done":true}︡
︠a57117e8-8946-4ff3-9266-385c16aefa3a︠
# Example 3 on ppt page 13
from scipy.stats import poisson
show(poisson.pmf(k=8, mu=4))
show(poisson.cdf(k=9,mu=4))
︡d31d4d51-95ed-442a-98ae-52f293102b15︡{"html":"$\\displaystyle 0.0297701813049$
"}︡{"html":"$\\displaystyle 0.991867757203$
"}︡{"done":true}︡
︠4a3630ce-421e-4cd0-9307-923b70fbb2af︠
# Example on ppt page 14-2
from scipy.stats import poisson,binom
show("Poisson: \t\t",poisson.pmf(5, mu=1))
show("Binomial:\t\t", binom.pmf(5,10000,0.0001))
︡ed6f9eb6-166f-440c-9b48-ec31980933f7︡{"html":"Poisson: \t\t $\\displaystyle 0.00306566200976$
"}︡{"html":"Binomial:\t\t $\\displaystyle 0.00306397596597$
"}︡{"done":true}︡
︠0f84e12c-5a46-4e63-a522-87fbbb8a177e︠
# Example on ppt page 14-3
from scipy.stats import poisson,binom
show("Poisson:",1-poisson.cdf(1, mu=150.0/500))
show("Binomial:", 1-binom.cdf(1,150,1.0/500))
︡7b0db564-5f84-4bc2-80da-391954e616fb︡{"html":"Poisson: $\\displaystyle 0.0369363131138$
"}︡{"html":"Binomial: $\\displaystyle 0.0367803266518$
"}︡{"done":true}︡
︠1dd94900-96cb-4088-85b8-1e4f5ac9274b︠
from scipy.stats import binom
from sage.plot.point import point
from sage.plot.line import line
n=5
p=0.1
t = text("cdf of B(%i,%g)"%(n,p), (n,1.1), rgbcolor=(1,0,0))
P = [(x,binom.cdf(x,n,p)) for x in range(-1,n+1)]
t = t+ point(P[1:],size=50)
for x, y in P:
t = t+ line([(x,y),(x+1,y)],thickness=2)
t.show()
︡c21a1103-9920-44ae-b6fb-150d3dabc400︡{"file":{"filename":"/home/user/.sage/temp/project-89fc128f-9cc1-4f1d-9784-c596ccec7fa3/628/tmp_tXdq0F.svg","show":true,"text":null,"uuid":"33bd8b64-16f9-45b4-b15c-dd264955c6fa"},"once":false}︡{"done":true}︡
︠ebea0f99-c86d-4f25-8549-7b296f11bd0d︠
from scipy.stats import norm
p=plot(norm.cdf,xmin=-3.5,xmax=3.5)
show(p+text("$\Phi$", (0.2,1)))
q=plot(norm.pdf,xmin=-3.5,xmax=3.5,rgbcolor=(0,.5,1))
show(q)
show(p+q)
︡c3f7754c-9591-4a9b-9f89-96648524b267︡{"file":{"filename":"/home/user/.sage/temp/project-89fc128f-9cc1-4f1d-9784-c596ccec7fa3/148/tmp_35mM9P.svg","show":true,"text":null,"uuid":"295f45cc-dd8c-41b2-99a6-ae375a2f83d0"},"once":false}︡{"file":{"filename":"/home/user/.sage/temp/project-89fc128f-9cc1-4f1d-9784-c596ccec7fa3/148/tmp_gYy9H9.svg","show":true,"text":null,"uuid":"0b817251-e9f4-4eca-86be-7569f722a59f"},"once":false}︡{"file":{"filename":"/home/user/.sage/temp/project-89fc128f-9cc1-4f1d-9784-c596ccec7fa3/148/tmp_Isj4j5.svg","show":true,"text":null,"uuid":"5d53c4ed-e419-4c64-95eb-c672f496c755"},"once":false}︡{"done":true}︡
︠114fd58f-3763-41fa-9d7b-c7ba3377d828︠
# Example on Section 6.5 ppt p33
from scipy.stats import norm
show(400*norm.ppf(0.9))
show(1-2*norm.cdf(-1))
show(1-2*norm.cdf(-2))
show(1-2*norm.cdf(-3))
︡958cb73f-829d-427d-9966-2695077c25ca︡{"html":"$\\displaystyle 512.620626218$
"}︡{"html":"$\\displaystyle 0.682689492137$
"}︡{"html":"$\\displaystyle 0.954499736104$
"}︡{"html":"$\\displaystyle 0.997300203937$
"}︡{"done":true}{"html":"$\\displaystyle 512.620626218$
"}︡{"html":"$\\displaystyle 0.682689492137$
"}︡{"html":"$\\displaystyle 0.954499736104$
"}︡{"html":"$\\displaystyle 0.997300203937$
"}︡{"done":true}︡
︠0d1f5a4c-4043-45ee-b388-e872b4fc5142s︠
# Example on Section 6.5 ppt 6
from scipy.stats import norm
show("Route 1",norm.cdf(70,loc=50,scale=10),", Route 2",norm.cdf(70,loc=60,scale=4))
show("Route 1",norm.cdf(60,loc=50,scale=10),", Route 2",norm.cdf(60,loc=60,scale=4))
︡545a3290-a64b-43ff-be65-ab0dd975c5a2︡{"html":"Route 1 $\\displaystyle 0.9772498680518208$ , Route 2 $\\displaystyle 0.9937903346742238$
"}︡{"html":"Route 1 $\\displaystyle 0.8413447460685429$ , Route 2 $\\displaystyle 0.5$
"}︡{"done":true}
︠ed81a69f-d62b-4973-af4c-8110d9061ca8s︠
# 2d normal distribution
from sage.plot.contour_plot import ContourPlot
︡e1e36c88-4078-4d86-ae95-036545241cf7︡{"done":true}
︠ef8b3734-0aea-4dfa-ae28-a597bfd8d56bs︠
x,y,u1, u2, s1,s2,rho=var("x y u1 u2 s1 s2 rho")
#def mvpdf(x,y,u1,u2,s1,s2,rho):
# assert(abs(rho)<1)
# return
multpdf = 1/(2*pi*sqrt(1-rho^2))*exp(-0.5/(1-rho^2)*((x-u1)^2/s1^2+(y-u2)^2/s2^2-2*rho*(x-u1)*(y-u2)/s1/s2))
multpdf1 = 30*multcdf.subs(u1=0,u2=0,s1=1,s2=1,rho=0)
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
w = plot3d(multpdf1,(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(60, 'rgbtuple'),optical=0.6,frame=False)
show(xaxis+yaxis+w)
multpdf1 = 30*multcdf.subs(u1=0,u2=0,s1=1,s2=1,rho=0.8)
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
w = plot3d(multpdf1,(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(60, 'rgbtuple'),optical=0.6,frame=False)
show(xaxis+yaxis+w)
multpdf1 = 30*multcdf.subs(u1=0,u2=0,s1=1,s2=1,rho=-0.8)
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
w = plot3d(multpdf1,(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(100, 'rgbtuple'),optical=0.6,frame=False)
show(xaxis+yaxis+w)
︡e641ad96-bea2-4798-b2ce-ea1bf69a238d︡{"stderr":"Error in lines 3-3\nTraceback (most recent call last):\n File \"/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py\", line 1188, in execute\n flags=compile_flags) in namespace, locals\n File \"\", line 1, in \nNameError: name 'multcdf' is not defined\n"}︡{"done":true}
︠8467981c-1c08-4016-bd0d-e4271d53f41c︠
k = var('k')
g = plot(exp(-ln(k)/k),(k,2,30))
show(g)
︡928cf82d-4aa3-471e-b637-a338ce547aa7︡{"file":{"filename":"/home/user/.sage/temp/project-89fc128f-9cc1-4f1d-9784-c596ccec7fa3/357/tmp_pL8khq.svg","show":true,"text":null,"uuid":"e70edd52-2bdd-4047-bf42-8a4756091085"},"once":false}︡{"done":true}︡
︠80aea1db-c263-40a2-8854-96a67955d4ff︠
show((2*120+2*130+3*125+110+135+140.0)/10)
D1 = [0.028, 0.121,0.234,0.267,0.200,0.103,0.037,0.009,0.001]
D2 = [0.001, 0.010,0.044, 0.117, 0.205, 0.247, 0.205, 0.117, 0.044, 0.010]
show(sum([i*p for i,p in enumerate(D1)]))
show(sum([i*p for i,p in enumerate(D2)]))
︡04518336-c42b-4a47-82b6-d9815c0dea2c︡
︠7d3a10c8-2b08-47a9-ac27-627b9858e9ff︠
D = [(-1,1/8),(0,1/4),(2,3/8),(3,1/4)]
show(sum([x*p for x,p in D]))
show(sum([x^2 *p for x,p in D]))
show(sum([(-2*x+1) *p for x,p in D]))
︡62ad95b7-1075-40be-a880-c3c4ec362421︡{"html":"$\\displaystyle \\frac{11}{8}$
"}︡{"html":"$\\displaystyle \\frac{31}{8}$
"}︡{"html":"$\\displaystyle -\\frac{7}{4}$
"}︡{"done":true}︡
︠f9cd07ea-2764-4234-9413-bcf9eca00a54︠
︡3d058483-4a2b-43ed-8d1c-8266e5579529︡
︠483425ba-2679-4c9e-87b7-ad474009020ds︠
x, y = var('x y')
def f(s,t,r):
return 1/(2*pi*s*t*sqrt(1-r^2))*exp(-(x^2/(s^2)-2*r*x*y/s/t+y^2/(t^2))/(2*(1-r^2)))
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
w = plot3d(30*f(1,1,0),(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(60, 'rgbtuple'),optical=0.6,frame=False)
show(xaxis+yaxis+w)
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
w = plot3d(30*f(1,1,0.5),(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(120, 'rgbtuple'),optical=0.6,frame=False)
show(xaxis+yaxis+w)
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
w = plot3d(30*f(1,1,-0.5),(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(120, 'rgbtuple'),optical=0.6,frame=False)
show(xaxis+yaxis+w)
︡00ee71e1-d8a4-43ea-a77a-c189d7fcc4b2︡{"file":{"filename":"6200beca-97fa-48de-91c9-8475cff669d1.sage3d","uuid":"6200beca-97fa-48de-91c9-8475cff669d1"}}︡{"file":{"filename":"1625e04c-fce9-457d-8554-15f0c3515f98.sage3d","uuid":"1625e04c-fce9-457d-8554-15f0c3515f98"}}︡{"file":{"filename":"2e1a263e-f47e-422e-b64b-2919ce71202e.sage3d","uuid":"2e1a263e-f47e-422e-b64b-2919ce71202e"}}︡{"done":true}
︠c39ed7a2-e529-458c-90da-6e49f3c95456s︠
from scipy.stats import binom
from sage.plot.point import point
from sage.plot.line import line
n=2
p=0.6
t = text("cdf of B(%i,%g)"%(n,p), (n,1.1), rgbcolor=(1,0,0))
P = [(x,binom.cdf(x,n,p)) for x in range(-1,n+1)]
for x, y in P:
t = t+ line([(x,y),(x+1,y)],thickness=2,rgbcolor=(0,0,1))
t = t+ point(P[1:],size=50)
for x,y in P[:-1]:
t=t+circle((x+1,y),0.03,fill=True,facecolor=(1,1,1))
t.show()
︡9ca3d3e7-e023-4fd0-8246-911cbc6c621c︡{"file":{"filename":"/home/user/.sage/temp/project-89fc128f-9cc1-4f1d-9784-c596ccec7fa3/468/tmp_hg12SB.svg","show":true,"text":null,"uuid":"9f24c713-3dee-408a-843f-c75311e2586b"},"once":false}︡{"done":true}
︠bd89cd3f-cca8-42bb-8419-210560c7e629︠