| Hosted by CoCalc | Download
# 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()
3D rendering not yet implemented
# 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)
zx=z_x = xye(xy)+e(xy)\displaystyle x y e^{\left(x y\right)} + e^{\left(x y\right)}
zy=z_y = x2e(xy)\displaystyle x^{2} e^{\left(x y\right)}
zxx=z_{xx} = xy2e(xy)+2 ye(xy)\displaystyle x y^{2} e^{\left(x y\right)} + 2 \, y e^{\left(x y\right)}
zyy=z_{yy} = x3e(xy)\displaystyle x^{3} e^{\left(x y\right)}
zxy=z_{xy} = x2ye(xy)+2 xe(xy)\displaystyle x^{2} y e^{\left(x y\right)} + 2 \, x e^{\left(x y\right)}
zyx=z_{yx} = x2ye(xy)+2 xe(xy)\displaystyle x^{2} y e^{\left(x y\right)} + 2 \, x e^{\left(x y\right)}
# 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())
zxx=z_{xx} = −2 x2(x2+y2)2+1x2+y2\displaystyle -\frac{2 \, x^{2}}{{\left(x^{2} + y^{2}\right)}^{2}} + \frac{1}{x^{2} + y^{2}}
zyy=z_{yy} = −2 y2(x2+y2)2+1x2+y2\displaystyle -\frac{2 \, y^{2}}{{\left(x^{2} + y^{2}\right)}^{2}} + \frac{1}{x^{2} + y^{2}}
zxx+zyy=z_{xx} + z_{yy} = 0\displaystyle 0
# 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)
1.041.98=1.04^{1.98}= 1.08075191020342\displaystyle 1.08075191020342
estimate = 1.08000000000000\displaystyle 1.08000000000000
# 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())
zxyz_{xy}= ∂∂xu(x,y)∂∂yu(x,y)D0,0(z)(u(x,y),v(x,y))+∂∂yu(x,y)∂∂xv(x,y)D0,1(z)(u(x,y),v(x,y))+∂∂xu(x,y)∂∂yv(x,y)D0,1(z)(u(x,y),v(x,y))+∂∂xv(x,y)∂∂yv(x,y)D1,1(z)(u(x,y),v(x,y))+∂2∂x∂yu(x,y)D0(z)(u(x,y),v(x,y))+∂2∂x∂yv(x,y)D1(z)(u(x,y),v(x,y))\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)
# 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()



# 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()
3D rendering not yet implemented
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()

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()
3D rendering not yet implemented
# 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])
[[x=0\displaystyle x = 0, y=0\displaystyle y = 0]]
Function value: [0\displaystyle 0]
[[x=2\displaystyle x = \sqrt{2}, y=12 2\displaystyle y = \frac{1}{2} \, \sqrt{2}, λ=(−72)\displaystyle \lambda = \left(-\frac{7}{2}\right)], [x=−2\displaystyle x = -\sqrt{2}, y=−12 2\displaystyle y = -\frac{1}{2} \, \sqrt{2}, λ=(−72)\displaystyle \lambda = \left(-\frac{7}{2}\right)], [x=1\displaystyle x = 1, y=(−1)\displaystyle y = \left(-1\right), λ=1\displaystyle \lambda = 1], [x=(−1)\displaystyle x = \left(-1\right), y=1\displaystyle y = 1, λ=1\displaystyle \lambda = 1]]
Function value: [212\displaystyle \frac{21}{2}, 212\displaystyle \frac{21}{2}, −3\displaystyle -3, −3\displaystyle -3]
# 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)
23\displaystyle \frac{2}{3} = 23\displaystyle \frac{2}{3}
# 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))
525\displaystyle \frac{52}{5}
458\displaystyle \frac{45}{8}
−sin⁡(1)+1\displaystyle -\sin\left(1\right) + 1
12\displaystyle \frac{1}{2}
# 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))
32 π\displaystyle \frac{3}{2} \, \pi
# 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))
19 (3 π−4)R3\displaystyle \frac{1}{9} \, {\left(3 \, \pi - 4\right)} R^{3}
# 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)`
2 a2\displaystyle 2 \, a^{2}
#################
### Chapter 6 ###
#################






# Example 6.2.1
P = binomial(100-3,4)*binomial(3,1)/binomial(100,5);
show(P,'=',float(P))
8936468\displaystyle \frac{893}{6468} = 0.138064316636\displaystyle 0.138064316636
# 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))
0.0142857142857\displaystyle 0.0142857142857
# 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))
2 students; 0.0027397260274 3 students; 0.00820416588478 4 students; 0.0163559124666 5 students; 0.0271355736998 6 students; 0.0404624836491 7 students; 0.056235703096 8 students; 0.0743352923517 9 students; 0.0946238338892 10 students; 0.116948177711 11 students; 0.141141378322 12 students; 0.167024788838 13 students; 0.194410275232 14 students; 0.223102512005 15 students; 0.252901319764 16 students; 0.283604005253 17 students; 0.315007665297 18 students; 0.346911417872 19 students; 0.379118526032 20 students; 0.411438383581 21 students; 0.443688335165 22 students; 0.475695307663 23 students; 0.507297234324 24 students; 0.538344257915 25 students; 0.568699703969 26 students; 0.598240820136 27 students; 0.626859282263 28 students; 0.654461472342 29 students; 0.680968537478 30 students; 0.706316242719 31 students; 0.730454633729 32 students; 0.75334752785 33 students; 0.774971854176 34 students; 0.79531686462 35 students; 0.814383238875 36 students; 0.83218210638 37 students; 0.848734008216 38 students; 0.864067821082 39 students; 0.878219664367 40 students; 0.891231809818 41 students; 0.903151611482 42 students; 0.914030471562 43 students; 0.923922855656 44 students; 0.932885368551 45 students; 0.940975899466 46 students; 0.948252843367 47 students; 0.954774402833 48 students; 0.960597972879 49 students; 0.965779609323 50 students; 0.970373579578 
# Example on ppt page 40
show(1-(1-0.004)^100)
0.330217428727354\displaystyle 0.330217428727354
# 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))
0.9936306176\displaystyle 0.9936306176
0.032773828125\displaystyle 0.032773828125
# 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))
0.00931615460852329\displaystyle 0.00931615460852329
# 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))
P(X=0) = 0.590490 P(X=1) = 0.328050 P(X=2) = 0.072900 P(X=3) = 0.008100 P(X=4) = 0.000450 P(X=5) = 0.000010 
# 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))
P(X=0) = 0.057648 P(X=1) = 0.197650 P(X=2) = 0.296475 P(X=3) = 0.254122 P(X=4) = 0.136137 P(X=5) = 0.046675 P(X=6) = 0.010002 P(X=7) = 0.001225 P(X=8) = 0.000066 P(X>=2) = 0.744702 
# 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)))
P(X≤0)=0.000045P(X\leq 0)=0.000045
P(X≤1)=0.000499P(X\leq 1)=0.000499
P(X≤2)=0.002769P(X\leq 2)=0.002769
P(X≤3)=0.010336P(X\leq 3)=0.010336
P(X≤4)=0.029253P(X\leq 4)=0.029253
P(X≤5)=0.067086P(X\leq 5)=0.067086
P(X≤6)=0.130141P(X\leq 6)=0.130141
P(X≤7)=0.220221P(X\leq 7)=0.220221
P(X≤8)=0.332820P(X\leq 8)=0.332820
P(X≤9)=0.457930P(X\leq 9)=0.457930
P(X≤10)=0.583040P(X\leq 10)=0.583040
P(X≤11)=0.696776P(X\leq 11)=0.696776
P(X≤12)=0.791556P(X\leq 12)=0.791556
P(X≤13)=0.864464P(X\leq 13)=0.864464
P(X≤14)=0.916542P(X\leq 14)=0.916542
P(X≤15)=0.951260P(X\leq 15)=0.951260
P(X≤16)=0.972958P(X\leq 16)=0.972958
P(X≤17)=0.985722P(X\leq 17)=0.985722
P(X≤18)=0.992813P(X\leq 18)=0.992813
P(X≤19)=0.996546P(X\leq 19)=0.996546
# 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))
[μ=0\displaystyle \mu = 0, μ=2\displaystyle \mu = 2]
0.0902235221577418\displaystyle 0.0902235221577418
0.0902235221577\displaystyle 0.0902235221577
# 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))
0.0297701813049\displaystyle 0.0297701813049
0.991867757203\displaystyle 0.991867757203
# 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))
Poisson: 0.00306566200976\displaystyle 0.00306566200976
Binomial: 0.00306397596597\displaystyle 0.00306397596597
# 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))
Poisson: 0.0369363131138\displaystyle 0.0369363131138
Binomial: 0.0367803266518\displaystyle 0.0367803266518
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()

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)

# 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))
512.620626218\displaystyle 512.620626218
0.682689492137\displaystyle 0.682689492137
0.954499736104\displaystyle 0.954499736104
0.997300203937\displaystyle 0.997300203937
0.682689492137\displaystyle 0.682689492137
0.954499736104\displaystyle 0.954499736104
0.997300203937\displaystyle 0.997300203937
# 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))
Route 1 0.9772498680518208\displaystyle 0.9772498680518208 , Route 2 0.9937903346742238\displaystyle 0.9937903346742238
Route 1 0.8413447460685429\displaystyle 0.8413447460685429 , Route 2 0.5\displaystyle 0.5
# 2d normal distribution
from sage.plot.contour_plot import ContourPlot

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)
Error in lines 3-3 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute flags=compile_flags) in namespace, locals File "", line 1, in <module> NameError: name 'multcdf' is not defined 

k = var('k')
g = plot(exp(-ln(k)/k),(k,2,30))
show(g)


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

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]))
118\displaystyle \frac{11}{8}
318\displaystyle \frac{31}{8}
−74\displaystyle -\frac{7}{4}



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)
3D rendering not yet implemented
3D rendering not yet implemented
3D rendering not yet implemented
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()