Autocovarianza en AM

Referencia: León García Ejemplo 9.9 p495, Gubner Ej10.35 p496

Sea X(t) = A cos(2πt), donde A es una variable aleatoria, con un comportamiento semejante a la figura.

Encontrar el valor esperado , la autocorrelación y autocovarianza de de X(t).


El valor esperado se calcula a continuación, note que la media varia respecto a t y que el valor es cero para valores de t donde cos(2πt) =0.

E[X(t)] = E[A cos(2\pi t)] E[X(t)] = E[A] cos(2\pi t)

La autocorrelación es:

R_X(t_1,t_2) = E[A cos(2\pi t_1) A cos(2\pi t_2)] = E[A^{2} cos(2\pi t_1) cos(2\pi t_2)] = E[A^{2}] cos(2\pi t_1) cos(2\pi t_2)

usando:

2 cos(x)cos(y) = cos(x-y) + cos(x+y) cos(x)cos(y) = \frac{ cos(x-y) + cos(x + y)}{2}

se reemplaza:

= E[A^{2}] \frac{1}{2}[cos(2\pi t_1 - 2\pi t_2) + cos(2\pi t_1 + 2\pi t_2)] R_X(t_1,t_2) = E[A^{2}] \frac{[cos(2\pi (t_1 - t_2)) + cos(2\pi (t_1 + t_2))]}{2}

se observa que el valor de autocorrelación depende de las diferencias de tiempo t1 y t2.


La autocovarianza es:

Cov_X(t_1,t_2) = R_X(t_1,t_2) - E[X(t_1)]E[X(t_2)] = E[A^{2}] cos(2\pi t_1) cos(2\pi t_2) - E[A] cos(2\pi t_1)E[A] cos(2\pi t_2) = E[A^{2}] cos(2\pi t_1) cos(2\pi t_2) - E[A]^2 cos(2\pi t_1)cos(2\pi t_2) = (E[A^{2}] - E[A]^2) cos(2\pi t_1)cos(2\pi t_2) = Var[A] cos(2\pi t_1)cos(2\pi t_2)

con el mismo procedimiento de cos(x)cos(y):

Cov_X(t_1,t_2) = Var[A] \frac{[cos(2\pi (t_1 - t_2)) + cos(2\pi (t_1 + t_2))]}{2}

Autocorrelación, Autocovarianza con variable tiempo

Referencia: León-García 9.2.2 p493

Se puede usar los momentos de muestras en el tiempo para parcialmente especificar un proceso aleatorio al resumir la información contenida en las cdf conjuntas.


Procesos aleatorios Contínuos en tiempo

Media:

m_X(t) = E[X(t)] = \int_{-\infty}^{\infty} x f_{X(t)}(x) dx

Varianza:

VAR[X(t)] = \int_{-\infty}^{\infty} ( x - m_X(t))^2 f_{X(t)}(x) dx

donde f_{X(t)}(x) es la pdf de X(t). Note que ambas son funciones determinísticas de tiempo.

Autocorrelación

R_X(t_1,t_2) = E[X(t_1,t_2)] = = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} xy f_{X(t_1),X(t_2)}(x,y) dx dy

Autocovarianza:

C_X(t_1,t_2) = E[\{X(t_1) - m_X(t_1) \} \{ X(t_2) - m_X(t_2)\} C_X(t_1,t_2) = R_X(t_1,t_2) - m_X(t_1) m_X(t_2)

coeficiente de correlación

\rho (t_1, t_2) = \frac{C_x(t_1,t_2)}{\sqrt{C_x(t_1,t_1)}\sqrt{C_x(t_2,t_2)}}

El coeficiente de correlación es la medida en la cual una variable aleatoria puede predecirse como una función lineal de otra.


Procesos aleatorios Discretos en tiempo

Media:

m_X(n) = E[X(n)]

Varianza:

VAR[X(n)] = E[(X(n) - m_X(n))^2]

Autocorrelación

R_X(n_1,n_2) = E[X(n_1,n_2)]

Autocovarianza:

C_X(n_1,n_2) = E[\{X(n_1) - m_X(n_1) \} \{ X(n_2) - m_X(n_2)\} C_X(n_1,n_2) = R_X(n_1,n_2) - m_X(n_1) m_X(n_2)

 

Funciones de una Variable aleatoria

Sea X una variable aleatoria y sea g(x) una función de valor real definida en el eje real.

Defina Y= g(X), esto es. Y está determinada por la evaluación de la función en g(x) en el valor que ha tomado la variable aleatoria X.  Entonces Y también es una variable aleatoria.

Las probabilidades de los valores para Y dependen de la función g(x) así como la función distribución acumulada de X.

Considere una función no lineal Y=g(X) como la que se muestra en la figura.

donde |dy| es la longitud del intervalo y < Y ≤ (y+dy).
De forma similar, la probabilidad que el evento en cada intervalo es aproximadamente

f_Y(y)=\left.\sum_{k} \frac{f_X(x)}{|dy/dx|} \right|_{x=x_k} =\left.\sum_{k} f_X(x) \left| \frac{dx}{dy} \right| \right|_{x=x_k}

Ejemplo: León García 4.30 p.176

Sea X el valor de las muestras de voltaje de una señal de voz, y suponga que X tiene una distribución uniforme en el intervalo [-4d,4d].

Sea Y = q(X), donde la función característica de entrada-salida de un cuantizador (convertidor analógico-digital) se muetra en la figura. Encuentre la función de probabilidad de masa para Y.

Solución: El evento {Y=q} para q en SY es equivalente al evento {X en Iq}, donde Iq es un intervalo de puntos equivalentes mapeados en representación al punto q. La pmf de Y se encuentra evaluando:
P[Y=q] = \int_{I_q} f_X(t) dt

Lo que permite ver fácilmente que la representación de un punto tiene un intervalo de longitud d mappeado en él. Entonces existirán ocho posibles salidas equiprobables, es decir, P[Y=q] = 1/8 para q en SY

Metodo de la Transformada

Referencia: Gubner 4.3 p159, León García 7.6.2 p 398

Los metodos de transformadas son muy útiles para cálculos que involucran derivadas e integrales de funciones. Muchos problemas involucran el uso de la «convolución» de dos funciones f1(x) * f2(x), cuyo cálculo se facilita si se tabaja con un método de transformadas.

Usar por ejemplo la transformada de Fourier, al cambiar de dominio convierte la convolución de funciones en una multiplicación, al resultado se le realiza la antitransformada y se obtiene el resultado buscado.

Función característica

Sea X una variable aleatoria contínua con función densidad de probabilidad f(x), entonces:

\Phi_{X} (\omega) = E\left[ e^{j \omega X} \right] = \int_{-\infty}^{\infty} f_{X}(x) e^{j\omega x} dx \text{donde } j=\sqrt{-1}

lo que también es la transformada de Fourier de f, la fórmula para invertir la transformada es:

f_{X}(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \Phi_{X}(\omega) e^{-j\omega x} d\omega

Es una transformación de la función desde el dominio del tiempo (t) al dominio de la frecuencia (ω).

Nota: En los libros de sistemas y señales, se define la transformada de Fourier de f por \int_{-\infty}^{\infty} e^{-j\omega x} dx . Para ser mas preciso, se debería decir φX(v) es la transformada de Fourier evaluada en -v.

Convolutions | Why X+Y in probability is a beautiful mess. 3Blue1Brown. 27 junio 2023.

cdf – Modulación AM

Referencia: Leon W Couch 4-2 p234, «El Aguacate» introducción

Modulación en Amplitud (AM)

La modulación es el proceso de codificación de la información fuente, sonido o moduladora, dentro de una señal pasabanda s(t), resultande o modulada. La señal modulada se obtiene de:

senal(t) = A_c[1+moduladora(t)] cos(\omega_c t) s(t) = A_c[1+m(t)] cos(\omega_c t)

donde:

\omega_c = 2\pi f_c

fc es la frecuencia de la portadora o «carrier».
Ac es la amplitud de la portadora.

Como ejemplo, si se quiere enviar una señal de sonido obtenida de un archivo.wav modulada en Amplitud, la señal m(t) será:

muestra_GuitarraCuerda.wav

Instrucciones en Python

# pmf de un sonido
# entrada es archivo01
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves
import scipy.stats as stats

# INGRESO 
# archivo01 = input('archivo de sonido 01: ' )
# k = int(input('muestras para ejemplo: '))
archivo01 = 'muestra_GuitarraCuerda.wav'
k = 500

# PROCEDIMIENTO
muestreo, sonido01 = waves.read(archivo01)

# Extrae un canal en caso de estéreo
canales = sonido01.shape
cuantos = len(canales)
canal = 0   
if (cuantos==1): # Monofónico
    uncanal = sonido01[:]  
if (cuantos>=2): # Estéreo
    uncanal = sonido01[:,canal]
    
moduladora = uncanal[0:k].astype(float)
dt = 1/muestreo
t  = np.arange(0,k*dt,dt)

# SALIDA GRAFICA
plt.plot(t,moduladora)
plt.title(' Moduladora m(t)')
plt.xlabel('t')
plt.ylabel('señal')
plt.plot()
plt.show()

Para el ejemplo, la señal de la portadora presentada tiene frecuencia de 5500 para que se pueda visualizar el efecto.
(Revisar frecuencias de portadoras AM estándares o ver el dial de un radio AM).

# Portadora:
fc = 5500
portadora = np.cos(2*np.pi*fc*t)

# SALIDA GRAFICA
plt.plot(t,portadora, color='orange')
plt.title(' Portadora')
plt.xlabel('t')
plt.ylabel('señal')
plt.plot()
plt.show()

Antes de aplicar la moduladora, se la normaliza para mantener la proporción en la gráfica

# normalizar y subir a positiva
moduladoranorm = moduladora/np.max(moduladora)
moduladora = (1+ moduladoranorm)

# Modular portadora
Ac = 1
modulada = Ac*moduladora*portadora

# SALIDA GRAFICA
plt.plot(t,moduladora,label='moduladora')
plt.plot(t,modulada,label='modulada')
plt.title(' Señal modulada S(t)')
plt.xlabel('t')
plt.ylabel('señal')
plt.legend()
plt.show()

cdf – una señal de sonido

Referencia: Leon W Couch apéndice B p675, «El Aguacate», pasillo introducción, Archivos de Audio.wav – Abrir, extraer una porción

De forma semejante al análisis para una señal triangular, se obtiene barriendo una ventana estrecha de Δx voltios de ancho, verticalmente a lo largo de las formas de onda y después midiendo la frecuencia relativa de la ocurrencia de voltajes en la ventana Δx.

El eje de tiempo se divide en n intervalos y la forma de onda aparece nΔx veces dentro de estos intervalos en la ventana Δx.

Se inicia abriendo el archivo 'muestra01_ElAguacateIntro.wav', seleccionando un canal en caso de ser estéreo. Para los cálculos se convierte a números reales ‘float’ de tipo estándar de python y para la gráfica se determina el rango de tiempo 't'.

# pmf de un sonido
# entrada es arcchivo01
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves
import scipy.stats as stats

# INGRESO 
# archivo01 = input('archivo de sonido 01:' )
archivo01 = 'muestra01_ElAguacateIntro.wav'

# PROCEDIMIENTO
muestreo, sonido01 = waves.read(archivo01)

# Extrae un canal en caso de estéreo
canales=np.shape(sonido01)
cuantos=len(canales)
canal = 0   
if (cuantos==1): # Monofónico
    uncanal=sonido01[:]  
if (cuantos>=2): # Estéreo
    uncanal=sonido01[:,canal]
    
senal=uncanal.astype(float)
dt=1/muestreo
n=len(senal)
t=np.arange(0,n*dt,dt)

Se grafica la señal en observación de uncanal:

# SALIDA GRAFICA
plt.plot(t,senal)
plt.title(' Sonido')
plt.xlabel('t')
plt.ylabel('señal')
plt.plot()
plt.show()


Aplicando luego el mismo algoritmo usado para señales triangulares y senoidales, manteniendo la forma discreta.

# Función de Probabilidad de Masa, PMF
m=40  # intervalos eje vertical

#PROCEDIMIENTO
relativa = stats.relfreq(senal, numbins = m )
deltax=relativa.binsize

# Eje de frecuencias, por cada deltax
senalmin=np.min(senal)
senalmax=np.max(senal)
senalrango=np.linspace(senalmin,senalmax,m)

# SALIDA
print('frecuencia relativa:')
print(relativa.frequency)
print('Rango de Señal')
print(senalrango)
frecuencia relativa:
[ 0.00023583  0.00073469  0.00089796  0.00053515  0.00084354  0.00179592  0.00245805  0.00306576  0.0039093   0.00504308  0.0066576   0.0099229  0.01693424  0.02671202  0.03568254  0.05160998  0.06758277  0.07323356  0.08575057  0.09217234  0.11898413  0.08494331  0.07535601  0.06454422  0.05382313  0.03821315  0.02429932  0.01573696  0.01075737  0.00767347  0.00576871  0.00365533  0.00245805  0.00176871  0.00145125  0.00103401  0.0008254   0.00106122  0.00137868  0.0004898 ]
Rango de Señal
[-32768.         -31087.61538462 -29407.23076923 -27726.84615385 -26046.46153846 -24366.07692308 -22685.69230769 -21005.30769231 -19324.92307692 -17644.53846154 -15964.15384615 -14283.76923077 -12603.38461538 -10923.          -9242.61538462  -7562.23076923  -5881.84615385  -4201.46153846  -2521.07692308   -840.69230769    839.69230769   2520.07692308   4200.46153846   5880.84615385   7561.23076923   9241.61538462  10922.          12602.38461538  14282.76923077  15963.15384615  17643.53846154  19323.92307692  21004.30769231  22684.69230769  24365.07692308  26045.46153846  27725.84615385  29406.23076923  31086.61538462  32767.        ]
# SALIDA Grafico de PMF
plt.bar(senalrango,relativa.frequency, width=deltax*0.8)
plt.xlabel('Amplitud de señal')
plt.ylabel('PMF')
plt.show()

# Función distribución acumulada
acumulada=np.cumsum(relativa.frequency)

# SALIDA CDF
plt.step(senalrango,relativa.frequency,label='pmf', where='post')
plt.step(senalrango,acumulada,label='cdf', where='post')
plt.xlabel('Amplitud de señal')
plt.title(' Función de distribuión acumulada , CDF')
plt.legend()
plt.show()

Para grabar los archivos puede usar np.savetxt(‘nombrearchivo.txt’, arreglodata )

Correlación

Referencia: León-García 5.6.2 p258, Gubner 2.4 p91

Correlación

La correlación entre dos variables aleatorias X y Y se define como E[XY].

La correlación determina cuando dos variables se encuentran linealmente relacionadas; es decir cuando una es función lineal de la otra.

R(X,Y) = E[XY]

Propiedades de la función de correlación

Simetría

R(X,Y) = R[Y,X] R(X,X) = E[X^2] \geq 0

Desigualdad de Cauchy-Schwarz

|R(X,Y)| = \sqrt{(E[X]^2 E[Y]^2)}

Covarianza

Retomando la función de covarianza de un proceso estocástico se muestra que:

Cov(X,Y) = E[XY] - E[X]E[Y] = = R(X,Y) - E[X]E[Y]

Coeficiente de correlación lineal

Al multiplicar una de las variables X o Y por un número se incrementa la covarianza, para una mejor medida se normaliza la covarianza y así tener los valores en una escala absoluta.

El coeficiente de correlación de X y Y se define por:

\rho_{X,Y} = \frac{Cov(X,Y)}{\sigma_X \sigma_Y} = \frac{E[XY]-E[X]E[Y]}{\sigma_X \sigma_Y} donde: \sigma_X=\sqrt{Var(X)} \sigma_Y=\sqrt{Var(Y)} -1 \leq \rho_{X,Y} \leq 1

El coeficiente de correlación es como máximo en magnitud 1.

Note que la correlación y el coeficiente de correlación no son lo mismo al resultar de la formula de covarianza.

Covarianza bivariada Ejemplo Ross 2_33

Referencia: Ejemplo Ross 2.33 p51

La función densidad conjunta de X,Y es

f(x,y) = \frac{1}{y} e^{-(y+x/y)} donde: 0<x , y<\infty

los intervalos también se expresan como: 0<x<∞  , 0<y<∞

a) Verifique que es una función de densidad conjunta

b) Encuentre la Cov(X,Y)


Solución propuesa

Para mostrar que f(x,y) es una función de densidad conjunta se debe mostar que es no negativa y que la integral doble es 1.

\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}f(x,y) dy dx = = \int_{0}^{\infty} \int_{0}^{\infty} \frac{1}{y} e^{-(y+x/y)} dy dx = = \int_{0}^{\infty} e^{-y} \int_{0}^{\infty} \frac{1}{y} e^{-x/y} dx dy

Resolviendo parcialmente el integral

\int_{0}^{\infty} \frac{1}{y} e^{-x/y} dx = \frac{1}{y} \int_{0}^{\infty} e^{-x/y} dx = = \frac{1}{y} \left. \frac{e^{-x/y}}{-\frac{1}{y}}\right|_{0}^{\infty} = -[e^{-\infty /y} - e^{0} ] = -(0-1) = 1

con el resultado se continua con el integral anterior:

= \int_{0}^{\infty} e^{-y} dy = \left. \frac{e^{-y}}{-1} \right|_{0}^{\infty} = - (e^{-\infty} - e^{0}) = = -(0-1) = 1

con lo que se comprueba que es una función densidad de probabilidad conjunta.


Covarianza (X,Y)

Se obtendrán primero las funciones marginales para determinar el valor esperado de cada una de ellas:

f_Y(y) = \int_{-\infty}^{\infty} f(x,y) dx = \int_{0}^{\infty} \frac{1}{y} e^{-(y+x/y)} dx = = \int_{0}^{\infty} \frac{1}{y} e^{-y} e^{-x/y} dx = e^{-y} \int_{0}^{\infty} \frac{1}{y}e^{-x/y} dx = e^{-y}

usando parte del integral anterior, se conoce que es igual a 1

se requiere el valor esperado E[y]:

E[Y] = \int_{-\infty}^{\infty} y f(y) dy = \int_{0}^{\infty} y e^{-y} dy =

integración por partes dv=e^{-y}dy , v = - e^{-y} , u=y, du=dy

\int u dv = uv - \int v du E[Y] = \left. - ye^{-y} \right|_{0}^{\infty} - \int_{0}^{\infty} (-e^{-y}) dy = -(\infty e^{-\infty} - 0 e^{-0}) - \left. e^{-y} \right|_{0}^{\infty} = = 0 -[e^{-\infty}- e^{0}] = -[0 -1] E[Y] =1

El valor esperado E[x] se obtiene como:

E[x] = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x f(x,y) dy dx = = \int_{0}^{\infty} \int_{0}^{\infty} x \frac{1}{y} e^{-(y+x/y)} dy dx = = \int_{0}^{\infty} e^{-y} \int_{0}^{\infty} \frac{x}{y} e^{-x/y} dx dy =

para el integral respecto a x:

\int_{0}^{\infty} \frac{x}{y} e^{-x/y} dx =

integración por partes dv=\frac{1}{y}e^{-x/y} dx , v = \frac{\frac{1}{y}}{-(1/y)} e^{-x/y}, u=x, du=dx

\int u dv = uv - \int v du = = \left. x(-e^{-x/y})\right|_{0}^{\infty} - \int_{0}^{\infty} - e^{-x/y} dx = = [- \infty ye^{-\infty/y} - (-0e^{0/y})] + \int_{0}^{\infty} e^{-x/y} dx = = [-0+0] + \left. \frac{e^{-x/y}}{-(1/y)}\right|_{0}^{\infty} = =-y[e^{-\infty/y} - e^{-0/y}] = y

reemplazando en E[x] y usando el resultado del integral de E[y]:

E[x] = \int_{0}^{\infty} y e^{-y} dy = 1

Para calcular el valor de E[XY]:

E[XY] = =\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} xy f(x,y) dy dx = =\int_{0}^{\infty} \int_{0}^{\infty} xy \frac{1}{y} e^{-(y+x/y)} dy dx = = \int_{0}^{\infty} y e^{-y} \int_{0}^{\infty} \frac{x}{y} e^{-x/y} dx dy = = \int_{0}^{\infty} y e^{-y} (y) dy = \int_{0}^{\infty} y^2 e^{-y} dy =

nuevamente por partes: dv = e^{-y} dy, u=y^2 se obtiene

E[XY] = \int_{0}^{\infty} y^2e^{-y} dy = =\left. -y^2 e^{-y} \right|_{0}^{\infty} + \int_{0}^{\infty} 2ye^{-y}dy = 2E[Y] =2

En consecuencia:

Cov(X,Y) = E[XY] - E[X]E[Y] = 2-(1)(1) = 1

Covarianza bivariada

Referencia: Ross 2.5.3 p50

Covarianza

La covarianza y varianza de dos variables aleatorias X y Y se definen como:

Cov(X,Y) = E[(X-E[X])(Y- E[Y])] = E[XY] - E[X]E[Y]

Si X y Y son independientes, entonces Cov(X,Y) = 0

Propiedades

  1. Cov(X,X) = Var(X)
  2. Cov(X,Y) = Cov(Y,X)
  3. Cov(cX,Y) = c Cov(X,Y)
  4. Cov(X,Y+Z) = Cov(X,Y) + Cov(X,Z)

 

Pdf Bivariada Ejercicio02

Referencia: León-García 5.17 p 253

Ejemplo

Encuentre P[ X + Y ≤ 1] de la funcion en el ejemplo 5.16 mostrada a continuación :

f_{X,Y}(x,y)= \begin{cases} c e^{-x} e^{-y} & , 0\leq y \leq x < \infty \\ 0 & ,\text{otro caso} \end{cases}

Solución

La la región para integración es [ X + Y ≤ 1] donde la pdf no es cero. Se obtine la probabilidad del evento al añadir(integrar) rectángulos infinitesimales de ancho dy como se indica en la figura:

X + Y \leq 1 Y \leq 1 - X P[ X + Y \leq 1] = \int_{0}^{1/2} \int_{y}^{1-y} 2 e^{-x} e^{-y} dx dy = \int_{0}^{1/2} 2 e^{-y} \int_{y}^{1-y} e^{-x} dx dy = \int_{0}^{1/2} 2 e^{-y} \left. [-e^{-x}] \right|_{y}^{1-y} dy = = \int_{0}^{1/2} 2 e^{-y} [-e^{-(1-y)}-(-e^{-y})] dy = \int_{0}^{1/2} [2 e^{-2y}- 2 e^{-y-(1-y)}] dy = = \int_{0}^{1/2} [2 e^{-2y}- 2 e^{-1}] dy = \left. \left[ 2\frac{e^{-2y}}{-2} - 2 e^{-1}y\right] \right|t_{0}^{1/2} = = [ - e^{-2 (1/2)} - 2 e^{-1} (1/2) ] - [ -e^{0}-0] = -e^{-1}-e^{-1} +1 P[ X + Y \leq 1] = 1- 2e^{-1} = 0.26424111765711533

que limita la figura que genera la función a:


Instrucciones en Python

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

# 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 and (x + y)<=1):
                z = c*np.exp(-x)*np.exp(y)
                Z[i,j] = z
    return(Z)

# PROGRAMA
# INGRESO
# Rango de evaluación
xa = 0
xb = 1.5
ya = 0
yb = 1.5
#muestras por eje
nx = 200
ny = 200

# PROCEDIMIENTO
# Matriz de evaluación
y = np.linspace(ya,yb,ny)
x = np.linspace(xa,xb,nx)
X,Y = np.meshgrid(x,y)
# Evalua la función
Z = fxydensidad(X,Y)

# Zona de integración
def arealimite(X):
    n = len(X)
    yinferior = np.zeros(n,dtype=float)
    ysuperior = np.zeros(n,dtype=float)
    for i in range(0,n,1):
        x = X[i]
        if (x<0.5):
            y = x
        if (x>=0.5):
            y = 1 - x
        ysuperior[i] = y
    
    return(yinferior, ysuperior)

# PROCEDIMIENTO
yinferior , ysuperior = arealimite(x) 

# SALIDA GRAFICAS
figura1 = plt.figure(1)
plt.plot(x,yinferior)
plt.plot(x,ysuperior)
plt.fill_between(x, yinferior, ysuperior,
                 where= (ysuperior>=yinferior))
plt.xlabel('x')

# SALIDA
figura2 = plt.figure(2)
grafica2 = figura2.add_subplot(1, 1, 1, projection='3d')
grafica2.plot_wireframe(X, Y, Z, rstride=10, cstride=10)
plt.show()