Estatísticas simples em hidrologia
Nelson Luís Dias
[email protected]
9 de setembro de 2009
Este exercício baseia-se em um arquivo de dados razoavelmente realista. O arquivo
patos.dat (clique com o botão direito do mouse no ícone para salvá-lo em disco:
)
contém dados de vazão diária no Rio dos Patos (afluente do Rio Iguaçu), medidos na estação hidrométrica com código ANA 64620000. A listagem 1 mostra um trecho do início do
arquivo de dados.
Listagem 1: Trecho do inicio do arquivo patos.dat
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
1930
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
1.00
0.60
0.60
0.60
0.60
0.60
0.60
0.60
0.60
1.00
2.88
O arquivo começa na linha 1, em 01/01/1930, mas o trecho mostra o início real das
medições, que foi em 20/05/1930. Note que até esta data há um código de erro (-1) para a
vazão, ou seja: não há dados. A organização do arqivo é muito simples (intencionalmente):
cada linha contém: ano, mês, dia do mês e vazão (em m3 s−1 ). Os dados vão até 31/12/1999.
1
Os dados serão processados em Python. Vamos começar abrindo o arquivo de dados e
armazenando as vazões diárias na variável qdado:
hcalpato.pyi≡
#!/usr/bin/python
# -*- coding: iso-8859-1 -*# –––––––––––––––––––––––––––––––––––––––
# –> calpato.py
#
# Nelson Luís Dias
# 2009-09-08T15:11:56
# 2009-09-08T15:12:01
# –––––––––––––––––––––––––––––––––––––––
# –––––––––––––––––––––––––––––––––––––––
# abre o arquivo de entrada
# –––––––––––––––––––––––––––––––––––––––
inp = open(’patos.dat’,’rt’)
# –––––––––––––––––––––––––––––––––––––––
# aloca um array com dados por ano e por dia
# –––––––––––––––––––––––––––––––––––––––
from numpy import zeros
nanos = 1999 - 1931 + 1
qdado = zeros((nanos,366),float)
from datetime import date
# –––––––––––––––––––––––––––––––––––––––
# lê a primeira linha
# –––––––––––––––––––––––––––––––––––––––
umal = inp.readline()
while umal:
# –––––––––––––––––––––––––––––––––––––––
# separa os campos
# –––––––––––––––––––––––––––––––––––––––
campo = umal.split()
ano = int(campo[0])
mes = int(campo[1])
dia = int(campo[2])
qq = float(campo[3])
# –––––––––––––––––––––––––––––––––––––––
# lê a próxima linha
# –––––––––––––––––––––––––––––––––––––––
umal = inp.readline()
# –––––––––––––––––––––––––––––––––––––––
# pula 1930
# –––––––––––––––––––––––––––––––––––––––
if ano > 1930:
# –––––––––––––––––––––––––––––––––––––––
# calcula o dia corrido do ano (dia juliano), e o índice do ano
# –––––––––––––––––––––––––––––––––––––––
iano = ano - 1931
dj = (date(ano,mes,dia) - date(ano,1,1)).days
qdado[iano,dj] = qq
2
Agora varremos o array de dados diários, e calculamos a vazão média. Para anos bissextos,
a média é sobre 366 dias; para anos normais, é sobre 365. Isto é feito automaticamente com
a variável ndias, e com a função isleap do módulo padrão calendar de Python:
hcalpato.pyi+≡
from calendar import isleap
ndias = (365,366)
# –––––––––––––––––––––––––––––––––––––––
# inicializa as vazões médias anuais com zeros
# –––––––––––––––––––––––––––––––––––––––
qano = zeros(nanos,float)
# –––––––––––––––––––––––––––––––––––––––
# loop para calcular as vazões médias anuais
# –––––––––––––––––––––––––––––––––––––––
for iano in range(nanos):
nd = ndias[isleap(iano+1931)]
for dj in range(nd):
qano[iano] += qdado[iano,dj]
qano[iano] /= nd
Informações úteis são a média e o desvio padrão das vazões médias anuais:
hcalpato.pyi+≡
# –––––––––––––––––––––––––––––––––––––––
# média e desvio-padrão não-enviesado
# –––––––––––––––––––––––––––––––––––––––
from math import sqrt
soma = 0.0
for i in range(nanos):
soma += qano[i]
qmed = soma/nanos
soma = 0.0
for i in range(nanos):
soma += (qano[i] - qmed)**2
qdvp = sqrt(soma/(nanos-1))
print ’med = ’, ’%5.2f’ % qmed
print ’dvp = ’, ’%5.2f’ % qdvp
3
Média e desvio-padrão são
q = 23,79 m3 s−1 ,
s∗q = 10,10 m3 s−1 .
(1)
(2)
Como vimos em sala de aula, se Q é normalmente distribuída com média µ e desviopadrão σ, então a variável “reduzida”
Z=
Q−µ
σ
(3)
possui fda
!
1 1
z
FZ (z) = + erf √ .
2 2
2
(4)
Portanto, pelo método dos momentos,
!
1 1
q − 23,79
√
FQ (q) = + erf
.
2 2
10,10 2
(5)
O próximo passo é ordenar a variável qano, e calcular posições de plotagem; com isto,
temos calculada a função distribuição acumulada empírica:
hcalpato.pyi+≡
qano.sort()
pp = zeros(nanos,float)
# –––––––––––––––––––––––––––––––––––––––
# calcula posições de plotagem Weibull
# –––––––––––––––––––––––––––––––––––––––
for i in range(nanos):
pp[i] = float(i+1)/(nanos+1)
# –––––––––––––––––––––––––––––––––––––––
# abre o arquivo com as posições de plotagem
# –––––––––––––––––––––––––––––––––––––––
fpp = open(’pp.dat’,’wt’)
for i in range(nanos):
fpp.write("%8.4f %8.4f\n" % (qano[i],pp[i]))
fpp.close()
4
A vazão média anual é o resultado da soma de um grande número de valores de vazão
média diária. Portanto, à primeira vista é razoável supor que a vazão média anual possua
distribuição normal. Vamos primeiro verificar quão “boa” é esta hipótese plotando (5) contra
a fda empírica. Para plotar as fda empírica e ajustada, usamos Gnuplot. O truque aqui é
gerar o script de Gnuplot de dentro do programa Python:
hcalpato.pyi+≡
# –––––––––––––––––––––––––––––––––––––––
# abre o arquivo gnuplot
# –––––––––––––––––––––––––––––––––––––––
plt = open(’calpato.plt’,’wt’)
pltstr = ’’’
set encoding iso_8859_1
set terminal postscript enhanced eps 18
Fq(x) = 0.5 + 0.5*erf( (x-23.79)/(10.10*sqrt(2.0)))
set xlabel ’vazão média anual q (m^3&{,}s^{/Symbol -1})’
set ylabel ’P\{ Q {/Symbol \243} q \}’
set nokey
set grid
set output ’pp.eps’
set format y ’%3.1f’
plot ’pp.dat’ with lines lt 1 lw 2, Fq(x) with lines lt 2 lw 2
!epstopdf pp.eps
exit
’’’
plt.write(pltstr)
plt.close()
from os import system
system(’gnuplot calpato.plt’)
5
1.0
0.9
0.8
P{ Q ≤ q }
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
0
10
20
30
40
50
60
3 −1
vazão média anual q (m s )
Figura 1: Fda empírica e fda da distribuição normal ajustada às vazões do Rio dos Patos
O resultado é mostrado na figura 1. Claramente, a cauda direita (vazões altas) da
fda empírica é mais “pesada” que a da normal ajustada: a distribuição normal subestima
as vazões com tempos de retorno (de excedência) altos. Na “vida real”, seria aconselhável
buscar uma distribuição que ajustasse melhor a cauda direita, se a preocupação fosse estimar
vazões médias anuais altas.
A tentação (na verdade, a indicação!), portanto, é tentar um teste de bondade de ajuste
aos dados. Usaremos Kolmogorov-Smirnov, assim como descrito por Benjamin e Cornell
(1970, p. 466 – 468). Atenção, o teste não vale quando a distribuição teórica é ajustada com
os próprios dados da distribuição empírica (conforme advertem Benjamin e Cornell (1970),
e também o artigo sobre o teste na Wikipedia
(http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test)!
Ao aplicá-lo para testar o ajuste de uma distribuição cujos parâmetros foram calculados a
partir dos próprios dados da amostra, portanto, nós estaremos fazendo uma concessão de
engenharia ao rigor matemático. O teste é o seguinte: define-se a distribuição empírica com
posição de plotagem
i
(6)
Fn (Xi ) = ,
n
e calcula-se a estatística
Dn = max |FX (Xi ) − Fn (Xi )|
(7)
i
onde por hipótese a amostra está ordenada: X1 ≤ X2 ≤ . . . ≤ Xn . Então, a variável (veja
(Press et al., 1992, Eq. 14.3.9))
√
√
K = ( n + 0.12 + 0.11/ n)Dn
(8)
6
possui fda de Kolmogorov:
P {K ≤ k} = 1 − 2
∞
X
(−1)i−1 e−2i
2 k2
.
(9)
i=1
Uma implementação despretensiosa de (9) é
hcalpato.pyi+≡
# –––––––––––––––––––––––––––––––––––––––
# fda de Kolmogorov; teste de Kolmogorov-Smirnov
# –––––––––––––––––––––––––––––––––––––––
from math import exp
def kol(n,dn):
# –––––––––––––––––––––––––––––––––––––––
# ajuste usando 14.3.9 de Numerical Recipes
# –––––––––––––––––––––––––––––––––––––––
adj = sqrt(n) + 0.12 + 0.11/sqrt(n)
# –––––––––––––––––––––––––––––––––––––––
# agora usa a fórmula assintótica
# –––––––––––––––––––––––––––––––––––––––
k = adj*dn
eps = 1.0e-6
pkold = 0.0
pknew = 1.0
i = 1
while abs(pknew - pkold) > eps:
pkold = pknew
pknew = pknew -2*(-1)**(i-1)*exp(-2*i**2*k**2)
print ’pknew = ’, pknew
i += 1
return pknew
Para aplicar o teste precisamos, também, recalcular as posições de plotagem:
hcalpato.pyi+≡
# –––––––––––––––––––––––––––––––––––––––
# recalcula a posição de plotagem para Kolmogorov-Smirnov
# –––––––––––––––––––––––––––––––––––––––
for i in range(nanos):
pp[i] = float(i+1)/nanos
O cálculo de Dn é simples, em Python. Primeiro, precisamos da mesma FQ (q) que já
tínhamos antes em Gnuplot, definida por (5):
hcalpato.pyi+≡
# –––––––––––––––––––––––––––––––––––––––
# fda normal
# –––––––––––––––––––––––––––––––––––––––
from scipy.special import erf
def Fq(q,med,dvp):
return (1.0 + erf((q - med)/(sqrt(2)*dvp)))/2.0
Agora que dispomos de FQ (q), calculamos as diferenças entre a distribuição acumulada
empírica e a teórica para todos os pontos amostrais; tiramos o seu máximo:
hcalpato.pyi+≡
# –––––––––––––––––––––––––––––––––––––––
# dn de Kolmogorov
# –––––––––––––––––––––––––––––––––––––––
dn = max([ abs(Fq(qano[i],qmed,qdvp) - pp[i]) for i in range(nanos)])
7
e verificamos então qual é a probabilidade de isto acontecer:
hcalpato.pyi+≡
# –––––––––––––––––––––––––––––––––––––––
# P { Dn > dn } = signif
# –––––––––––––––––––––––––––––––––––––––
signif = 1.0 - kol(nanos,dn)
print ’n = ’,nanos
print ’dn = ’, dn
print ’signif = ’, signif
8
O resultado do teste é desapontador: dn = 0,0839, e P {Dn > dn } = 0,698. O teste
de Kolmogorov-Smirnov não consegue diferenciar nossos dados anuais de uma distribuição
normal com os parâmetros dados, apesar da evidência visual na cauda direita. O programa
calpato.py inteiro está na listagem 2. O seu código está no anexo calpato.py (
).
Listagem 2: Listagem completa de calpato.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#!/ usr / bin / python
# -* - coding : iso -8859 -1 -* # -----------------------------------------------------------------------------# --> calpato . py
#
# Nelson Luís Dias
# 2009 -09 -08 T15 :11:56
# 2009 -09 -08 T15 :12:01
# -----------------------------------------------------------------------------# -----------------------------------------------------------------------------# abre o arquivo de entrada
# -----------------------------------------------------------------------------inp = open ( ' patos . dat ' , ' rt ')
# -----------------------------------------------------------------------------# aloca um array com dados por ano e por dia
# -----------------------------------------------------------------------------from numpy import zeros
nanos = 1999 - 1931 + 1
qdado = zeros (( nanos ,366) , float )
from datetime import date
# -----------------------------------------------------------------------------# lê a primeira linha
# -----------------------------------------------------------------------------umal = inp . readline ()
while umal :
# -----------------------------------------------------------------------------# separa os campos
# -----------------------------------------------------------------------------campo = umal . split ()
ano = int ( campo [0])
mes = int ( campo [1])
dia = int ( campo [2])
qq = float ( campo [3])
# -----------------------------------------------------------------------------# lê a próxima linha
# -----------------------------------------------------------------------------umal = inp . readline ()
# -----------------------------------------------------------------------------# pula 1930
# -----------------------------------------------------------------------------if ano > 1930:
# -----------------------------------------------------------------------------# calcula o dia corrido do ano ( dia juliano ) , e o índice do ano
# -----------------------------------------------------------------------------iano = ano - 1931
dj = ( date ( ano , mes , dia ) - date ( ano ,1 ,1)). days
qdado [ iano , dj ] = qq
from calendar import isleap
ndias = (365 ,366)
# -----------------------------------------------------------------------------# inicializa as vazões médias anuais com zeros
# -----------------------------------------------------------------------------qano = zeros ( nanos , float )
# -----------------------------------------------------------------------------# loop para calcular as vazões médias anuais
# -----------------------------------------------------------------------------for iano in range ( nanos ):
nd = ndias [ isleap ( iano +1931)]
for dj in range ( nd ):
9
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
qano [ iano ] += qdado [ iano , dj ]
qano [ iano ] /= nd
# -----------------------------------------------------------------------------# média e desvio - padrão não - enviesado
# -----------------------------------------------------------------------------from math import sqrt
soma = 0.0
for i in range ( nanos ):
soma += qano [ i ]
qmed = soma / nanos
soma = 0.0
for i in range ( nanos ):
soma += ( qano [ i ] - qmed )**2
qdvp = sqrt ( soma /( nanos -1))
print ' med = ', '%5.2 f ' % qmed
print ' dvp = ', '%5.2 f ' % qdvp
qano . sort ()
pp = zeros ( nanos , float )
# -----------------------------------------------------------------------------# calcula posições de plotagem Weibull
# -----------------------------------------------------------------------------for i in range ( nanos ):
pp [ i ] = float ( i +1)/( nanos +1)
# -----------------------------------------------------------------------------# abre o arquivo com as posições de plotagem
# -----------------------------------------------------------------------------fpp = open ( ' pp . dat ' , ' wt ')
for i in range ( nanos ):
fpp . write ("%8.4 f %8.4 f \ n " % ( qano [ i ] , pp [ i ]))
fpp . close ()
# -----------------------------------------------------------------------------# abre o arquivo gnuplot
# -----------------------------------------------------------------------------plt = open ( ' calpato . plt ' , ' wt ')
pltstr = '''
set encoding iso_8859_1
set terminal postscript enhanced eps 18
Fq ( x ) = 0.5 + 0.5* erf ( (x - 2 3. 79 ) /( 10 .1 0 * sqrt (2.0)))
set xlabel ' vazão média anual q ( m ^3&{ ,} s ^{/ Symbol -1}) '
set ylabel 'P \{ Q {/ Symbol \243} q \} '
set nokey
set grid
set output ' pp . eps '
set format y '%3.1 f '
plot ' pp . dat ' with lines lt 1 lw 2 , Fq ( x ) with lines lt 2 lw 2
! epstopdf pp . eps
exit
'''
plt . write ( pltstr )
plt . close ()
from os import system
system ( ' gnuplot calpato . plt ')
# -----------------------------------------------------------------------------# fda de Kolmogorov ; teste de Kolmogorov - Smirnov
# -----------------------------------------------------------------------------from math import exp
def kol (n , dn ):
# -----------------------------------------------------------------------------# ajuste usando 14.3.9 de Numerical Recipes
# -----------------------------------------------------------------------------adj = sqrt ( n ) + 0.12 + 0.11/ sqrt ( n )
# -----------------------------------------------------------------------------# agora usa a fórmula assintótica
# -----------------------------------------------------------------------------k = adj * dn
eps = 1.0 e -6
pkold = 0.0
pknew = 1.0
10
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
i = 1
while abs ( pknew - pkold ) > eps :
pkold = pknew
pknew = pknew -2*( -1)**( i -1)* exp ( -2* i **2* k **2)
print ' pknew = ', pknew
i += 1
return pknew
# -----------------------------------------------------------------------------# recalcula a posição de plotagem para Kolmogorov - Smirnov
# -----------------------------------------------------------------------------for i in range ( nanos ):
pp [ i ] = float ( i +1)/ nanos
# -----------------------------------------------------------------------------# fda normal
# -----------------------------------------------------------------------------from scipy . special import erf
def Fq (q , med , dvp ):
return (1.0 + erf (( q - med )/( sqrt (2)* dvp )))/2.0
# -----------------------------------------------------------------------------# dn de Kolmogorov
# -----------------------------------------------------------------------------dn = max ([ abs ( Fq ( qano [ i ] , qmed , qdvp ) - pp [ i ]) for i in range ( nanos )])
# -----------------------------------------------------------------------------# P { Dn > dn } = signif
# -----------------------------------------------------------------------------signif = 1.0 - kol ( nanos , dn )
print 'n = ', nanos
print ' dn = ', dn
print ' signif = ', signif
Referências
Benjamin, J. R. e Cornell, C. A. (1970). Probability, Statistics and Decision for Civil Engineers. McGraw-Hill, New York.
Press, W. H., Teukolsky, S. A., Vetterling, W. T., e Flannery, B. P. (1992). Numerical
Recipes in C. Cambridge University Press, Cambridge, UK.
11