Referencia: León García Ejercicio 5.16 pag 252
Ejemplo
Encuentre la constante c de normalización y las pdf marginales de la siguiente función:
Solución
la función es válida en la región mostrada:
La constante c se encuentra cumpliendo la condición de normalización:
Determinar las marginales, conociendo que c=2:
Tarea: Verificar que las funciones marginales cumplen que el integral es 1
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # Función evaluada def fxydensidad(X,Y): n,m =np.shape(X) Z=np.zeros(shape=(n,m),dtype=float) c=2 for i in range(0,n,1): for j in range(0,m,1): x=X[i,j] y=Y[i,j] if (y>=0 and y<=x): z=c*np.exp(-x)*np.exp(y) Z[i,j]=z return(Z) # PROGRAMA --------- # INGRESO # Rango de evaluación xa = 0 xb = 4 ya = 0 yb = 4 # muestras por eje nx = 500 ny = 500 # PROCEDIMIENTO # Matriz de evaluación y = np.linspace(ya,yb,ny) x = np.linspace(xa,xb,nx) X,Y = np.meshgrid(x,y) # Evalúa la función Z = fxydensidad(X,Y) # SALIDA figura = plt.figure(1) grafica = figura.add_subplot(1, 1, 1, projection='3d') grafica.plot_wireframe(X, Y, Z, rstride=10, cstride=10) plt.show()
# Zona de integración def arealimite(X): n = len(X) yinferior = np.zeros(n,dtype=int) ysuperior = X return(yinferior, ysuperior) # PROCEDIMIENTO yinferior, ysuperior = arealimite(x) # SALIDA plt.plot(x, yinferior) plt.plot(x, ysuperior) plt.fill_between(x, yinferior, ysuperior, where=(ysuperior>=yinferior)) plt.show()
# Forma de las marginales def marginalx(x): fx = 2*np.exp(-x)*(1-np.exp(-x)) return(fx) # PROCEDIMIENTO deltax = (xb-xa)/nx fx = marginalx(x) Fx = np.cumsum(fx*deltax) integrax = np.sum(fx*deltax) # SALIDA plt.plot(x,fx, label='f(x)') plt.plot(x,Fx, label='F(x)') plt.xlabel('x') plt.legend() plt.show()
# Forma de las marginales def marginaly(y): fy = 2*np.exp(-2*y) return(fy) # PROCEDIMIENTO deltay = (yb-ya)/ny fy = marginaly(y) Fy = np.cumsum(fy*deltay) integray = np.sum(fy*deltay) # SALIDA print(' El integral sobre el area es: ', integray) plt.plot(y,fy, label='f(y)') plt.plot(y,Fy, label='F(y)') plt.xlabel('y') plt.legend() plt.show()