Distribuciones Estadsticas en Hidrologa: Distribucin
Normal
Ing. Ivn Arturo Ayala Bizarro, ayalabizarro@[Link]
Modificado 16 de abril de 2016
Hidrologa General
Distribucin Normal
Tambin denominada campana de Gauss (la distribucin de probabilidad es similar a una campana)
Funcin de densidad de probabilidad f (x)
f (x) =
Estandarizando con z =
1 x 2
1
e 2 ( )
2
(1)
1 z2
f (z) = e 2
2
(2)
Funcin de distribucin acumulada F (x)
F (x) =
e 2 (
1
) dx
(3)
1
F (z) =
2
z 2
2
dz
(4)
Solucin aproximada segn Abramowitz y Stegun (1965)
F (z) = 1 f (z)(0.4361836w 0.1201676w2 + 0.937280w3 )
(5)
donde:
w=
1
1 + 0.33267|Z|
(6)
si Z < 0 la funcin de distribucin acumulada es 1 F (Z)
Parmetros
Para ambos mtodo (Mxima Verosimilitud y Momentos) son iguales
n
=
"
1X
xi
n i=1
n
1 X
2
=
(xi )
n 1 i=1
(7)
#1/2
(8)
Donde: x es la varaible independiente, el parmetro de localizacin igual a la media, parmetro de escala
igual a la desviacin estandard y n tamao de la muestra.
1
Recepcin de datos
Importar los datos de la variable.
# Definir la ubicacin de la carpeta de trabajo
setwd('C:/Users/abia/OneDrive/IVAN/IVAN_PROGRAMACION_R/Rpubs_IvanAyala/02.-Distribuciones Estadsticas')
Datos<- [Link]('Data_Pmax.txt', header=TRUE)
x<-sort(Datos$Pmax)
#summary(x)
Distribucin Emprica
Existen muchas relaciones para representar una Distribucin Emprica (DE). Para precipitaciones y caudales
en hidrologa, recomiendan el suso de la distribucin propuesta por Weilbull, cuya relacin es.
Pe =
m
n+1
(9)
donde m es el nmero de orden correspondiente a cada valor de la serie ordenados de menor a mayor, y n es
el tamao o nmero 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 = "Distribucin Emprica")
grid()
0.14285714
0.28571429
0.42857143
0.57142857
0.71428571
0.85714286
0.0
0.2
0.4
0.6
0.8
1.0
Distribucin Emprica
15
20
25
30
35
x
Clculo de los parmetros
# Clculo de parmetros Mtodo de Momento y Mxima verosimilitud
mu <- mean(x);mu
# Parmetro de Posicin.
## [1] 29.10244
sigma <- sd(x);sigma
Parmetro de Escala.
## [1] 8.527353
Clculo de la funcin de distribucin acumulada Fx
Mediante la ecuacin 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
40
45
F_AS[i]<-1-fz*(0.4361836*w-0.1201676*w**2+0.937280*w**3)
}
F_AS
##
##
##
##
##
##
##
[1]
[7]
[13]
[19]
[25]
[31]
[37]
if(z<0){
F_AS[i]<-1-F_AS[i]
}
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
# Funcin de Distribucin 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 lnea
legend('topleft', legend=c('Distribucin Normal','Distribucin Emprica'),
lty=1:2, col=c('red', 'blue'), cex=0.8, [Link]=4)
grid()
1.0
0.2
0.4
0.6
0.8
Probabilidad Acumulada
Distribucin Normal
Distribucin Emprica
15
20
25
30
35
40
45
Clculo de probabilidades
Las probabilidas se asignan a partir del periodo de retorno Tr al cual se desea someter. Mayor informacin
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
Clculo de Cuantiles
Q<-qnorm(p, mu, sigma)
Rsult<-[Link](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 (aos)",
main = "Cuantiles Vs Periodo de retorno" )
grid()
45
40
30
35
Cuantil
50
55
Cuantiles Vs Periodo de retorno
200
400
600
800
1000
Periodo de retorno (aos)
Test para la distribucin terica
El test de prueba, se realiza mediante el test de Smirnov-Kolmogorov. Una forma de validar, es crear una
serie de datos con una distribucin conocida y realizar la prueba de Smirnov Kolmogorov.
x2<-rbeta(1000, shape1=50, shape2=2, ncp = 0)
tsk<-[Link](x2,"pnorm",mean(x2),sd(x2));tsk
##
## One-sample Kolmogorov-Smirnov test
##
## data: x2
6
# Distribucin Beta
## D = 0.091837, p-value = 9.448e-08
## alternative hypothesis: two-sided
# Ahora el test para beta
tsk<-[Link](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<-[Link](x,"pnorm",mu,sigma);tsk
# Smirnov Kolmogorov
## Warning in [Link](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?