%matplotlib inline
import random
import numpy as np
import scipy as sc
import scipy.stats as st
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
Datos de oleaje y viento y normal 2D¶
Vamos a examinar datos de viento y oleaje de una boya cerca del puerto de Cádiz.
- Hm0 : Altura significante Espectral (m)
- DirM : Direccion Media de PROCEDENCIA del Oleaje (0=N,90=E)
- VelV : Velocidad Media del Viento (m/s)
- DirV : Direccion Media de PROCEDENCIA del Viento (0=N,90=E)
También tiene otros datos que vamos a ignorar porque no están presentes o no nos interesan.
simar=pd.read_csv('/home/OYE/SIMAR_6006052.csv')
simar.head()
simar.info()
simar.describe()
Histogramas de altura de ola y velocidad del viento¶
plt.figure(figsize=(9,9))
plt.hist(simar['Hm0'], bins=200)
plt.xlabel('Hm0')
plt.title('altura media de ola Hm0')
plt.show()
#Tenemos que quitar primero los datos NA (Not Available)
notnas = simar['VelV'].notna()
velV_sin_na = simar['VelV'][notnas]
#en vez del histograma de la librería matplotlib,
#usamos el distplot de la librería "seaborn"
sns.distplot(velV_sin_na)
No parecen normales:
- no son simétricas
- sólo toman valores positivos
Distribución de la velocidad del viento¶
Dibujamos el histograma de velocidad del viento junto con la función de densidad de la Normal con la misma media y desviación típica.
Recordamos que la normal que mejor ajusta a los datos en el sentido de "máxima verosimilitud" es la $\operatorname{Normal}(\mu, \sigma^2)$ que tiene su parámetro $\mu$ igual a la media muestral: $$ \bar{x} = \frac{1}{n}\sum_{j=1}^n x_j $$ y su parámetro $\sigma^2$ igual a la varianza muestral: $$ \hat{\sigma}^2 = \frac{1}{n}\sum_{j=1}^n (x_j - \bar{x})^2 $$
VelV = simar['VelV'][simar['VelV'].notna()]
meanV, stdV = VelV.mean(), VelV.std()
plt.figure(figsize=(16,12))
plt.hist(VelV, bins=30, density=1,alpha=0.9)
xs = np.arange(meanV-3*stdV,meanV+3*stdV,0.01)
N = st.norm(loc=meanV, scale=stdV)
ys = N.pdf(xs)
lines = plt.plot(xs,ys, lw=3)
plt.xlabel('VelV')
plt.show()
La distribución de Rayleigh ajusta mucho mejor, si escogemos el parámetro adecuadamente. Para la distribución de Rayleigh, el valor de máxima verosimilitud para el parámetro de escala es:
$${\displaystyle {\widehat {\sigma }}\approx {\sqrt {{\frac {1}{2N}}\sum _{i=1}^{N}x_{i}^{2}}}}$$plt.figure(figsize=(16,12))
plt.hist(VelV, bins=50, density=1,alpha=0.9, label='Velocidad del viento')
scale_data = np.sqrt((VelV**2).sum()/(2*len(VelV)))
R = st.rayleigh(scale=scale_data)
xs = np.arange(meanV-3*stdV,meanV+5*stdV,0.01)
ys = R.pdf(xs)
lines = plt.plot(xs,ys, lw=3, label='Rayleigh(%.4f)'%scale_data)
plt.legend()
plt.xlabel('VelV')
plt.show()
El resultado es igual de revelador usando la función de distribución empírica
def dibuja_distribucion_empirica(muestra, dist=None):
mean, std = dist.mean(), dist.std()
minx = min(min(muestra)-1, mean-3*std)
maxx = max(max(muestra)+1, mean+3*std)
xs = np.concatenate([[minx], sorted(muestra), [maxx]])
pasoy = 1/len(muestra)
ys = np.concatenate( [np.arange(0,1,pasoy), [1,1]])
plt.step(xs, ys, where='post', label='cdf empírica')
if dist:
xs = np.arange(minx,maxx,std*0.01)
ys = dist.cdf(xs)
plt.plot(xs,ys, label = 'cdf')
plt.legend()
Comparamos con la normal
data = simar['VelV'][simar['VelV'].notna()]
mean, std = data.mean(), data.std()
N = st.norm(loc=mean, scale=std)
dibuja_distribucion_empirica(data,N)
plt.xlabel('VelV')
plt.show()
Comparamos con la distribución de Rayleigh
data = simar['VelV'][simar['VelV'].notna()]
mean, std = data.mean(), data.std()
R = st.rayleigh(scale=scale_data)
dibuja_distribucion_empirica(data,R)
plt.xlabel('VelV')
plt.show()
Maximum likelihood con st.rayleigh.fit(data, floc=0)¶
En lugar de usar las fórmulas de la tabla (o de la wikipedia, o de cualquier libro de texto), podemos encontrar los valores de los parámetros que maximizan la verosimilitud numéricamente, usando la sintaxis:
loc0, scale0 = st.rayleigh.fit(VelV, floc=0)
Donde:
locyscaleson los dos parámetros quescipy.statsreconoce parascipy.stats.rayleighfloc=0indica que queremos fijarloca0, porque no queremos que la desplace.VelVson los datos
#Obtenemos un valor muy similar al anterior
st.rayleigh.fit(VelV, floc=0)
scale_data
Distribución de la altura del oleaje¶
HM = simar['Hm0'][simar['Hm0'].notna()]
meanH, stdH = HM.mean(), HM.std()
plt.figure(figsize=(8,6))
plt.hist(HM, bins=50, density=1,alpha=0.9)
xs = np.arange(meanH-3*stdH,meanH+5*stdH,0.01)
loc0, scale0 = st.rayleigh.fit(HM, floc=0)
R = st.rayleigh(scale=scale0)
ys = R.pdf(xs)
lines = plt.plot(xs,ys, lw=3)
plt.xlabel('Hm0')
plt.show()
R = st.rayleigh(scale=meanH/np.sqrt(np.pi/2))
dibuja_distribucion_empirica(HM,R)
plt.xlabel('Hm0')
plt.show()
La distribución de Rayleigh ajusta mejor que una normal, pero no ajusta tan bien como lo hacía para la velocidad del viento. Se suele modelizar como una distribución de Weibull, que es algo más difícil de ajustar:
- En la wikipedia puedes leer que no se conoce una expresión cerrada para el estimador de máxima verosimilitud.
- El método
st.weibull_min.fit(HM, floc=0)devuelve unos valores muy malos de los parámetros, porque el métodofitnecesita un punto de partida decente. Cualquier númeroc_inimayor que uno es un buen punto de partida parast.weibull_min.fit(HM, c_ini, floc=0)
c_ini = 2
c0, loc0, scale0 = st.weibull_min.fit(HM, c_ini, floc=0)
print(c0, loc0, scale0)
W = st.weibull_min(c=c0, loc=loc0, scale=scale0)
dibuja_distribucion_empirica(HM,W)
plt.xlabel('Hm0')
plt.show()
En este caso concreto tampoco ajusta mucho mejor que una Rayleigh :-/
Relaciones entre variables¶
Una correlación interesante: ¿a mayor velocidad del viento mayor altura del oleaje?
simar['VelV'].mean(), simar['Hm0'].mean()
simar.plot.scatter(y = 'Hm0', x = 'VelV', s=10, alpha=0.1)
Observamos correlación positiva entre las variables.
simar[['Hm0','VelV']].corr()
Intentamos ajustar el vector (Altura significativa, velocidad del viento) por una normal multivariante...
percentiles_canonicos = [
1-2*st.norm(loc=0,scale=1).cdf(-k) for k in (1,2,3)
]
def plotMN(MN, sample_size=1000, npoints=100):
#Muestra aleatoria de la normal multivariable
means, Sigma = MN.mean, MN.cov
# Create grid coordinates for plotting
std1, std2 = np.sqrt([Sigma[0,0], Sigma[1,1]])
X1 = np.linspace(means[0] - 4*std1, means[0] + 4*std1, npoints)
X2 = np.linspace(means[1] - 4*std2, means[1] + 4*std2, npoints)
xx, yy = np.meshgrid(X1,X2, indexing='xy')
Z = np.zeros((X1.size,X2.size))
for i,x1i in enumerate(X1):
for j, x2j in enumerate(X2):
Z[j,i] = MN.pdf([x1i,x2j])
#Queremos dibujar exactamente los contornos que capturan
# los percentiles habituales (~68%,95%,99%):
#https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Interval
factor = 1/(2*np.pi*np.sqrt(np.linalg.det(Sigma)))
E = st.expon(scale=2).ppf
levels = [
factor*np.exp(-0.5*E(p))
for p in reversed(percentiles_canonicos)
]
plt.contour(xx, yy, Z, cmap=plt.cm.Set2_r, levels=levels)
plt.contourf(xx, yy, Z, cmap=plt.cm.Set2_r, alpha=0.4, levels=levels)
La normal multivariante que mejor ajusta a los datos es la que tiene parámetro $\mathbf{\mu}$ igual a la media muestral $\bar{\mathbf{x}}=\frac{1}{n}\sum_{j=1}^n \mathbf{x}_j$ y parámetro $\Sigma$ igual a la matriz de covarianzas de la muestra $\frac{1}{n}\sum_{j=1}^n (\mathbf{x}_j-\bar{\mathbf{x}})^t\cdot (\mathbf{x}_j-\bar{\mathbf{x}})$.
El conjunto de puntos (Altura significativa, velocidad del viento) no parece sacado de una normal multivariante.
data = simar[['Hm0','VelV']]
means = data.mean()
Sigma = data.cov()
print(means)
print(Sigma)
MN = st.multivariate_normal(mean=means, cov=Sigma)
fig = plt.figure(figsize=(9,9))
fig.suptitle('Altura significativa vs velocidad del viento', fontsize=20)
plotMN(MN, npoints=200)
plt.scatter(data[data.columns[0]], data[data.columns[1]], s=6, alpha=0.3)
plt.xlabel('Hm0')
plt.ylabel('VelV')
plt.show()
Descomponemos el viento en componentes¶
Según la wikipedia:
One example where the Rayleigh distribution naturally arises is when wind velocity is analyzed into its orthogonal 2-dimensional vector components. Assuming that each component is uncorrelated, normally distributed with equal variance, and zero mean, then the overall wind speed (vector magnitude) will be characterized by a Rayleigh distribution.
Descomponemos el viento en sus componentes NS y EW, y vemos que la distribución se parece razonablemente a una normal multivariante, cuyos parámetros hemos encontrado igual que antes.
VNS_EW = pd.DataFrame({
'VEW': np.sin(simar['DirV']*np.pi/180)*simar['VelV'],
'VNS': np.cos(simar['DirV']*np.pi/180)*simar['VelV'],
})
data = VNS_EW
means = data.mean()
Sigma = data.cov()
print(means)
print(Sigma)
MN = st.multivariate_normal(mean=means, cov=Sigma)
fig = plt.figure(figsize=(9,9))
fig.suptitle('Componentes NS-EW del viento', fontsize=20)
plotMN(MN, npoints=200)
plt.scatter(data['VEW'], data['VNS'], s=4, alpha=0.1)
plt.xlabel('VelV (E-W)')
plt.ylabel('VelV (N-S)')
plt.show()
Sin embargo, para los datos de oleaje, para los que los datos de altura significativa de la ola no sigue una distribución de Rayleigh, la normal multivariante no es una buena aproximación para las componentes direccionales.
El mecanismo que genera este fenómeno es otro...
HNS_EW = pd.DataFrame({
'HNS': np.cos(simar['DirM']*np.pi/180)*simar['Hm0'],
'HEW': np.sin(simar['DirM']*np.pi/180)*simar['Hm0'],
})
data = HNS_EW
means = data.mean()
Sigma = data.cov()
print(means)
print(Sigma)
MN = st.multivariate_normal(mean=means, cov=Sigma)
fig = plt.figure(figsize=(9,9))
fig.suptitle('Componentes NS-EW del oleaje', fontsize=20)
plotMN(MN, npoints=200)
plt.scatter(data['HEW'], data['HNS'], s=4, alpha=0.1)
plt.xlabel('Hm0 (E-W)')
plt.ylabel('Hm0 (N-S)')
plt.show()
Estudiamos las correlaciones entre las componentes del viento y el oleaje. Aunque la normal multivariante no es un buen modelo para el oleaje, hay correlaciones importantes entre las componentes NS del viento y el oleaje, y no las hay entre las componentes NS de una y EW de la otra.
data = pd.DataFrame({
'HNS': np.cos(simar['DirM']*np.pi/180)*simar['Hm0'],
'VNS': np.cos(simar['DirV']*np.pi/180)*simar['VelV'],
'HEW': np.sin(simar['DirM']*np.pi/180)*simar['Hm0'],
'VEW': np.sin(simar['DirV']*np.pi/180)*simar['VelV'],
})
print(data.corr())
Es una pena que los datos de esta boya no traigan el oleaje desglosado en "mar de viento" y "mar de fondo". En fin, podríamos seguir afinando, pero lo vamos a dejar aquí...
Ejercicio¶
Ajusta una normal multivariante a los datos:
- Tm02 : Periodo Medio Espectral Momentos 0 y 2 (s)
- Tp : Periodo de pico espectral (s)
Preguntas
- ¿Hay correlación, positiva o negativa, entre ambas variables?
- ¿Es la normal multivariante un buen modelo?
data = simar[['Tp', 'Tm02']]
data.mean()
Ejercicio¶
Ajusta una normal multivariante a los datos:
- log(Tm02) : logaritmo del Periodo Medio Espectral Momentos 0 y 2 (s)
- log(Tp) : logaritmo del Periodo de pico espectral (s)
Preguntas
- ¿Hay correlación, positiva o negativa, entre ambas variables?
- ¿Es la normal multivariante un buen modelo?
logTpTm = pd.DataFrame({
'logTp': np.log(simar['Tp']),
'logTm': np.log(simar['Tm02']),
})
logTpTm.mean()