Distribuciones Estadísticas en Hidrología: Distribución Normal Ing. Iván Arturo Arturo Ayala Bizarro,
[email protected] Bizarro,
[email protected] Modificado 16 de abril de 2016 Hidrología General
Distribución Normal También denominada campana de Gauss (la distribución de probabilidad es similar a una campana) •
Función de densidad de probabilidad f (x) f (x) =
Estandarizando con z =
1 √ e σ 2π
−
2
(
x−µ σ
)
2
(1)
x−µ σ
√ 12π e
f (z ) = •
1
−z 2
(2)
2
Función de distribución acumulada F (x) x
1 F (x) = σ 2π
√ √
F (z ) =
e
−
1 2
(
x−µ σ
2
) dx
(3)
−∞
z
1 2π
e
−z2 2
(4)
dz
−∞
Solución aproximada según Abramowitz y Stegun (1965) F (z ) = 1
− f (z)(0.4361836w − 0.1201676w
2
+ 0.937280w3 )
(5)
donde: w =
1 1 + 0 .33267 Z
si Z < 0 la 0 la función de distribución acumulada es 1 •
(6)
| |
− F (Z )
Parámetros Para ambos método (Máxima Verosimilitud y Momentos) son iguales
σ =
n
1
− µ =
n
n
i=1
n
1
−1 i
(7)
xi
(xi
µ)
2
1/2
(8)
=1
Donde: x es la varaible independiente, µ el parámetro de localización igual a la media, σ parámetro de escala igual a la desviación estandard y n tamaño de la muestra. 1
Recepción de datos Importar los datos de la variable. # Definir la ubicación de la carpeta de trabajo setwd( C:/Users/abia/OneDrive/IVAN/IVAN_PROGRAMACION_R/Rpubs_IvanAyala/02.-Distribuciones Estadísticas
Datos<- read.table ( Data_Pmax.txt , header=TRUE) x<-sort (Datos$Pmax) #summary(x)
Distribución Empírica Existen muchas relaciones para representar una Distribución Empírica (DE). Para precipitaciones y caudales en hidrología, recomiendan el suso de la distribución propuesta por Weilbull, cuya relación es. P e =
m
(9)
n+1
donde m es el número de orden correspondiente a cada valor de la serie ordenados de menor a mayor, y n es el tamaño o número total de la muestra. n<-length (x) P_E <- c() for(i in 1:n) { P_E[i] <- i/(n+1) } P_E ## ## ## ## ## ## ##
[1] [7] [13] [19] [25] [31] [37]
0.02380952 0.16666667 0.30952381 0.45238095 0.59523810 0.73809524 0.88095238
0.04761905 0.19047619 0.33333333 0.47619048 0.61904762 0.76190476 0.90476190
0.07142857 0.21428571 0.35714286 0.50000000 0.64285714 0.78571429 0.92857143
0.09523810 0.23809524 0.38095238 0.52380952 0.66666667 0.80952381 0.95238095
0.11904762 0.26190476 0.40476190 0.54761905 0.69047619 0.83333333 0.97619048
plot(x,P_E, pch=18, col="blue", type="b", lty=1 ,ylab = "P",
main = "Distribución Empírica") grid()
2
0.14285714 0.28571429 0.42857143 0.57142857 0.71428571 0.85714286
Distribución Empírica 0 . 1 8 . 0 6 . 0 P 4 . 0 2 . 0 0 . 0
15
20
25
30
35 x
Cálculo de los parámetros # Cálculo de parámetros Método de Momento y Máxima verosimilitud # Parámetro de Posición. mu <- mean(x);mu
## [1] 29.10244 sigma <- sd(x);sigma
#
Parámetro de Escala.
## [1] 8.527353
Cálculo de la función de distribución acumulada •
Mediante la ecuación de Abramowitz y Stegun (1965) F_AS <- c() for(i in 1:n) { z<-(x[i]-mu)/sigma fz <- 1/(sqrt(2*pi))*exp(-0.5*z^2) w <- 1./(1+0.33267*abs(z))
3
F x
40
45
F_AS[i]<-1-fz*(0.4361836*w-0.1201676*w**2+0.937280*w**3) if(z<0){ F_AS[i]<-1-F_AS[i] } } F_AS ## ## ## ## ## ## ## •
[1] [7] [13] [19] [25] [31] [37]
0.04908434 0.12275747 0.32361276 0.43928759 0.59715931 0.75540283 0.92862637
0.05276376 0.12759032 0.34927917 0.44855050 0.60620816 0.75540283 0.93915405
0.08680805 0.15935489 0.35799717 0.45784211 0.63299495 0.76271016 0.96885482
0.09645090 0.20577245 0.36678971 0.48584033 0.65050912 0.77344284 0.97687728
0.10685901 0.22287071 0.39357183 0.50925349 0.68456401 0.85962628 0.98450532
0.10685901 0.27090858 0.40716308 0.52328022 0.70510218 0.86222498
0.10687041 0.22286919 0.39357770 0.50924183 0.68457105 0.85961587 0.98451571
0.10687041 0.27090162 0.40717213 0.52326461 0.70510989 0.86221436
Mediante comandos R # Función de Distribución Normal Acumulada F_N <- pnorm(x, mu, sigma);F_N
## ## ## ## ## ## ##
[1] [7] [13] [19] [25] [31] [37]
0.04908523 0.12276891 0.32360634 0.43930299 0.59715127 0.75540733 0.92861932
0.05276584 0.12760158 0.34927595 0.44856699 0.60620224 0.75540733 0.93914958
0.08681779 0.15936323 0.35799544 0.45785912 0.63299495 0.76271373 0.96886035
0.09646168 0.20577363 0.36678966 0.48585371 0.65051230 0.77344491 0.97688563
plot (x, F_N, type="l", pch=19, col="red",lty=2, xlab="",
ylab="Probabilidad Acumulada" , main="") lines(x, P_E, pch=18, col="blue", type="b", lty=1) # Agregar una línea legend ( topleft , legend= c( Distribución Normal , Distribución Empírica ), lty=1:2, col= c( red , blue ), cex=0.8, text.font=4) grid ()
4
0 . 1 a d a l u m u c A d a d i l i b a b o r P
Distribución Normal Distribución Empírica
8 . 0 6 . 0 4 . 0 2 . 0
15
20
25
30
35
40
45
Cálculo de probabilidades Las probabilidas se asignan a partir del periodo de retorno T r al cual se desea someter. Mayor información
en: Dr. Victor Miguel Ponce Tr<-c(2,5,10,15,20,25,50,100,200,250,500,1000) # Periodos de retorno p<-1-1/Tr;p
## ##
[1] 0.5000000 0.8000000 0.9000000 0.9333333 0.9500000 0.9600000 0.9800000 [8] 0.9900000 0.9950000 0.9960000 0.9980000 0.9990000
Cálculo de Cuantiles Q<-qnorm(p, mu, sigma) Rsult<-data.frame(Tr,p,Q);Rsult
## ## ## ## ## ## ##
1 2 3 4 5 6
Tr 2 5 10 15 20 25
p 0.5000000 0.8000000 0.9000000 0.9333333 0.9500000 0.9600000
Q 29.10244 36.27924 40.03068 41.90273 43.12869 44.03116
5
## ## ## ## ## ##
7 50 8 100 9 200 10 250 11 500 12 1000
0.9800000 0.9900000 0.9950000 0.9960000 0.9980000 0.9990000
46.61548 48.94003 51.06744 51.71757 53.64554 55.45394
plot(Tr,Q, pch=23, col="#000099", bg="#FF6666", type="b", lty=1 ,
ylab = "Cuantil", xlab = "Periodo de retorno (años)" , main = "Cuantiles Vs Periodo de retorno" ) grid()
Cuantiles Vs Periodo de retorno 5 5 0 5 l i t n a u C
5 4 0 4 5 3 0 3
0
200
400
600
800
1000
Periodo de retorno (años)
Test para la distribución teórica El test de prueba, se realiza mediante el test de Smirnov-Kolmogorov . Una forma de validar, es crear una serie de datos con una distribución conocida y realizar la prueba de Smirnov Kolmogorov. x2<-rbeta(1000, shape1=50, shape2=2, ncp = 0) tsk<-ks.test(x2,"pnorm",mean(x2),sd(x2));tsk ## ## One-sample Kolmogorov-Smirnov test ## ## data: x2
6
# Distribución Beta
## D = 0.091837, p-value = 9.448e-08 ## alternative hypothesis: two-sided # Ahora el test para beta tsk<-ks.test(x2,"pbeta",shape1=50, shape2=2, ncp = 0);tsk
## ## One-sample Kolmogorov-Smirnov test ## ## data: x2 ## D = 0.025208, p-value = 0.5488 ## alternative hypothesis: two-sided
Por lo tanto, la pruebas para la serie x es: tsk<-ks.test(x,"pnorm",mu,sigma);tsk
# Smirnov Kolmogorov
## Warning in ks.test(x, "pnorm", mu, sigma): ties should not be present for ## the Kolmogorov-Smirnov test ## ## One-sample Kolmogorov-Smirnov test ## ## data: x ## D = 0.06752, p-value = 0.9921 ## alternative hypothesis: two-sided
. . . comentarios?
7
