AJO_1617 Ejemplos 03.2017

hace 198 días por sevillad1617

Estimación de media y varianza en el laboratorio

mu = 84.963345 inte = [(85.39,85.62),(85.02,85.28),(85.12,85.35),(85.08,85.20),(84.82,85.19),(84.35,84.45),(84.65,84.95),(85.00,85.26)] p0 = sum([ line([(l[0],1+i),(l[1],1+i)],thickness=2) for i,l in enumerate(inte) ]) p1 = Graphics() p1 = line([(mu,0),(mu,len(inte)+1)],color="red",legend_label="$\mu$") (p0+p1).show(figsize=5) 
       
mu = 84.963345 inte = [(85.42,85.51),(84.36,85.23),(85.12,85.34),(84.58,85.14),(84.17,84.63),(84.92,85.33),(84.88,85.41),(84.84,85.29),(84.80,85.31)] p0 = sum([ line([(l[0],1+i),(l[1],1+i)],thickness=2) for i,l in enumerate(inte) ]) p1 = Graphics() p1 = line([(mu,0),(mu,len(inte)+1)],color="red",legend_label="$\mu$") (p0+p1).show(figsize=5) 
       
sigma2 = 0.497565^2 inte = [(0.02,0.55),(0.03,1.33),(0.03,0.88),(0.17,7.31),(0.11,5.03),(0.08,3.88),(0.15,6.68),(0.11,4.81)] p0 = sum([ line([(l[0],1+i),(l[1],1+i)],thickness=2) for i,l in enumerate(inte) ]) p1 = Graphics() p1 = line([(sigma2,0),(sigma2,len(inte)+1)],color="red",legend_label=u"¿ $\sigma^2$ ?") (p0+p1).show(figsize=5) 
       
 
       

Gráficas de distribuciones normales bidimensionales

import scipy.stats from mpl_toolkits.mplot3d import Axes3D from matplotlib.ticker import LinearLocator, FormatStrFormatter import matplotlib.pyplot as plt from matplotlib import cm import numpy as np #medias, varcovar = [11.5, -0.2], [[2., .6], [.6, 0.5]] #medias, varcovar = [11.5, -0.2], [[3., 0.3], [0.3, 0.5]] medias, varcovar = [20,5], matrix([[6.0,2], [2,1]]) #medias, varcovar = [20,5], diagonal_matrix([6.,1]) #medias, varcovar = [20,5], identity_matrix(2) rv = scipy.stats.multivariate_normal(medias, varcovar) 
       
medias, varcovar = [20,5], matrix([[6.0,0], [0,1]]) rv = scipy.stats.multivariate_normal(medias, varcovar) e = [sqrt(w).n() for w in matrix(varcovar).eigenvalues()] e = map(sqrt,matrix(varcovar).eigenvalues()) smax = max(e) smin = min(e) try: alfa = arctan(varcovar[0][1]/(smax^2-varcovar[0][0])) except ZeroDivisionError: alfa = 0 def eli(p,co="blue"): c = sqrt(r.qchisq(p,2))._sage_() return ellipse(medias, c*smax, c*smin, pi/2-alfa, color=co) verelipses = False n = 20000 puntos = rv.rvs([n]) nube_antigua = point(puntos,aspect_ratio=1,color="green",size=1) elipses = eli(0.39,"blue")+eli(0.95,"green") if verelipses else Graphics() origen = point((0,0),size=0) if True else Graphics() view(medias) view(varcovar) (nube_antigua+elipses+origen).show(figsize=7, aspect_ratio=1) 
       

                                
                            

                                
medias, varcovar = [20,5], matrix([[6.0,2], [2,1]]) rv = scipy.stats.multivariate_normal(medias, varcovar) e = [sqrt(w).n() for w in matrix(varcovar).eigenvalues()] e = map(sqrt,matrix(varcovar).eigenvalues()) smax = max(e) smin = min(e) try: alfa = arctan(varcovar[0][1]/(smax^2-varcovar[0][0])) except ZeroDivisionError: alfa = 0 def eli(p,co="blue"): c = sqrt(r.qchisq(p,2))._sage_() return ellipse(medias, c*smax, c*smin, pi/2-alfa, color=co) verelipses = False n = 20000 puntos = rv.rvs([n]) nube = point(puntos,aspect_ratio=1,color="red",size=1) elipses = eli(0.39,"blue")+eli(0.95,"green") if verelipses else Graphics() origen = point((0,0),size=0) if True else Graphics() view(medias) view(varcovar) (nube+nube_antigua+elipses+origen).show(figsize=7, aspect_ratio=1) 
       

                                
                            

                                
 
       
# see http://matplotlib.org/examples/axes_grid/scatter_hist.html for a 2D alternative fig = plt.figure(figsize=[8,6]) ax = fig.gca(projection='3d') cx, cy = 4, 4 xmin, xmax = medias[0]-cx*sqrt(varcovar[0][0]), medias[0]+cx*sqrt(varcovar[0][0]) xmin, xmax = xmin.n(), xmax.n() ymin, ymax = medias[1]-cy*sqrt(varcovar[1][1]), medias[1]+cy*sqrt(varcovar[1][1]) ymin, ymax = ymin.n(), ymax.n() X, Y = np.mgrid[xmin:xmax:0.1r, ymin:ymax:0.05r] pos = np.dstack((X,Y)) rv = scipy.stats.multivariate_normal(medias, varcovar) Z = rv.pdf(pos) margin = 2 ax.set_xlim(xmin-margin,xmax+margin) ax.set_ylim(ymin-margin,ymax+margin) n = 20000 X, Y = zip(*rv.rvs([n])) ax.scatter(X,Y,depthshade=False,s=0.01,color="red") hX = np.histogram(X,density=True,bins=20) hY = np.histogram(Y,density=True,bins=20) verhist = True if verhist: ax.bar(hX[1][0:-1],hX[0],ax.set_ylim()[1],zdir="y",linewidth=0) ax.bar(hY[1][0:-1],hY[0],ax.set_xlim()[0],zdir="x",linewidth=0) ax.zaxis.set_major_locator(LinearLocator(5)) ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f')) zmax = max(max(hX[0]),max(hY[0])) ax.set_zlim(0, zmax*1.1 if verhist else 1) # ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([1, 1, 0.25, 1])) # [scale_x, scale_y, scale_z, 1] ax.view_init(elev=20) ax.set_aspect(1,adjustable='box-forced') plt.savefig('') plt.close('all') 
       
 
       

Propagación lineal

# X,Y son N(0,1) # Z=X+Y # caso 1: independientes SXY = matrix([ [1,0] , [0,1] ]) A = matrix([ [1,1] ]) A*SXY*transpose(A) 
       
[2]
[2]
# X,Y son N(0,1) # Z=X+Y # caso 2: X=Y 
       
# X,Y son N(0,1) # Z=X+Y # caso 3: X=-Y 
       
 
       
# Z = A+B+C # mu_Z = mu_A + mu_B + mu_C = 200 sA = 0.0002 sB = 0.01 sC = 0.0002 S1 = diagonal_matrix([sA^2,sB^2,sC^2]) A = matrix([1,1,1]) S2 = A*S1*transpose(A) S2 
       
[0.000100080000000000]
[0.000100080000000000]
sZ = sqrt(S2[0,0]) sZ 
       
0.0100039992003198
0.0100039992003198
200-2*sZ,200+2*sZ 
       
(199.979992001599, 200.020007998401)
(199.979992001599, 200.020007998401)