Commit d6e9bca8 authored by Rafael Peretti Pezzi's avatar Rafael Peretti Pezzi

Estimativas de desvios das médias

parent 201e7a38
# -*- coding:utf-8 -*-
# Programa de visualização e estatísticas de dados meteorológicos
# Programa de visualização e estatísticas de dados meteorológicos
# coletados pela estação Meteorolog do Centro de Tecnologia Acadêmica
# Dados no formato DataHora Temperatura Umidade_do_Ar Solo Pressão Luminosidade
# Dados no formato DataHora Temperatura Umidade_do_Ar Solo Pressão Luminosidade
# DataHora no formato AAAAMMDDHHmmSS onde
# AAAA = Ano
# MM = Mês
......@@ -16,7 +16,7 @@
# Autor: Rafael Pezzi - Centro de Tecnologia Acadêmica - IF/UFRGS
#
# Este programa é softare livre, disponível sob os termos da
# GNU Affero General Public License V3
# GNU Affero General Public License V3
# http://www.gnu.org/licenses/agpl-3.0.html
import numpy as np
......@@ -27,9 +27,17 @@ import random as random
import datetime as datetime
import sys , os
# Realizar médias de n medidas aleatórias:
#n_medias = [1,5,10,50,100,1000,10000]
n_medias = [5,10,50,100,1000]
#n_medias = [5,10,50] # Possibilidades reduzidas para debug
#Número de medias a realizar por Monte Carlo
n_MC=5000
print "Estatísticas Meteorológicas Meteorolog - Centro de Tecnologia Acadêmica"
# Testa se o número de argumentos está correto
# Testa se o número de argumentos está correto
if len(sys.argv) < 2:
sys.stderr.write('Uso: python '+ sys.argv[0]+' arquivo_de_log\n' )
sys.exit(1)
......@@ -44,7 +52,7 @@ seleciona_intervalo = False
#seleciona_intervalo = True
# Defina aqui o intervalo de tempo a se analisado
# TODO: Incluir seleção pela linha de comando de execução
# TODO: Incluir seleção pela linha de comando de execução
sel_Inicio = datetime.datetime(2013,03,01,00,00,00)
sel_Fim = datetime.datetime(2013,03,31,23,59,59)
......@@ -52,7 +60,7 @@ sel_Fim = datetime.datetime(2013,03,31,23,59,59)
data, temp, u_ar, u_solo, pressao, lum = np.loadtxt(sys.argv[1], unpack=True,
converters={ 0: mdates.strpdate2num('%Y%m%d%H%M%S')})
## Seria útil incluir um algoritmo de filtragem para remover pontos expurios
## Seria útil incluir um algoritmo de filtragem para remover pontos expurios
## observados em algumas partes por mil medidas de temperatura
# Exibe informações dos dados
......@@ -65,8 +73,8 @@ if (seleciona_intervalo == True):
# Cria listas temporárias para seleção
# Acho que pode ser otimizado, usando iteradores direto em arrays
# crítico para series de dados gigantes
selec_data=[]; selec_temp=[]; selec_u_ar=[]
selec_u_solo=[]; selec_pressao=[]; selec_lum=[]
selec_data=[]; selec_temp=[]; selec_u_ar=[]
selec_u_solo=[]; selec_pressao=[]; selec_lum=[]
for j in enumerate(data):
if mdates.date2num(sel_Inicio) <= j[1] <= mdates.date2num(sel_Fim):
selec_data.append(data[j[0]])
......@@ -168,7 +176,7 @@ plt.figure()
plt.xticks(rotation=70)
plt.plot_date(x=data, y=temp, fmt="r.")
# Define os limites temporais do gráfico
# Define os limites temporais do gráfico
# útil para evitar escalas exageradas devido a dadas incorretas - ocorreu na estação sem relógio
plt.xlim(data[0], data[len(data)-1])
plt.ylabel(u"Temperatura (°C)")
......@@ -217,10 +225,15 @@ plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.2)
#'''
# Análise Monte Carlo para as médias
# Para fins de ilustração didática do teorema do limite central
# Para fins de ilustração didática do teorema do limite central e desvio da média
# Cria variaveis para calculo das médias
for tentativas in [5,10,50]:
temp_std_n=[]
umid_std_n=[]
pressao_std_n=[]
lum_std_n=[]
for tentativas in n_medias:
temp_med = 0.
temp_med_vec=[]
umid_med = 0.
......@@ -229,7 +242,7 @@ for tentativas in [5,10,50]:
pressao_med_vec=[]
lum_med = 0.
lum_med_vec=[]
for j in range(5000):
for j in range(n_MC):
for n in range(tentativas):
amostra=int(random.random()*len(data))
lum_med+=lum[amostra]
......@@ -245,28 +258,81 @@ for tentativas in [5,10,50]:
umid_med_vec.append(umid_med/tentativas)
pressao_med_vec.append(pressao_med/tentativas)
lum_med_vec.append(lum_med/tentativas)
temp_med = 0
umid_med = 0
umid_med =0
lum_med = 0
pressao_med = 0
# Convertentndo listas para matrizes para melhorar desempenho (pode - deve - ser feito anteriormente no programa)
temp_med_vec=np.asarray(temp_med_vec)
umid_med_vec=np.asarray(umid_med_vec)
pressao_med_vec=np.asarray(pressao_med_vec)
lum_med_vec=np.asarray(lum_med_vec)
# Pegando as médias dos resultados de monte carlo
temp_med = temp_med_vec.mean()
umid_med = umid_med_vec.mean()
lum_med = lum_med_vec.mean()
pressao_med = pressao_med_vec.mean()
lum_med_vec=np.asarray(lum_med_vec)
print " "
print str(tentativas)+" medidas:"
temp_med_std=temp_med_vec.std()
temp_std_n.append(temp_med_std)
print "Desvio Temp Média para MC com "+str(tentativas)+" medidas: "+str(temp_med_std)
# umid_med_vec.append(umid_med/tentativas)
umid_med_std=umid_med_vec.std()
umid_std_n.append(umid_med_std)
print "Desvio Umidade Média para MC com "+str(tentativas)+" medidas: "+str(umid_med_std)
# pressao_med_vec.append(pressao_med/tentativas)
pressao_med_std=pressao_med_vec.std()
pressao_std_n.append(pressao_med_std)
print "Desvio Pressão Média para MC com "+str(tentativas)+" medidas: "+str(pressao_med_std)
lum_med_std=lum_med_vec.std()
lum_std_n.append(lum_med_std)
print "Desvio lumin. Média para MC com "+str(tentativas)+" medidas: "+str(lum_med_std)
plt.figure()
plt.title(u"Temperatura - MC Média de: "+str(tentativas)+" medidas")
plt.hist(temp_med_vec, bins=50)
plt.title(u"Temp. MC "+str(tentativas)+" medidas: sigma = "+str("%0.3f" % temp_med_std)+u" | Média = "+str(temp_med))
n, bins, patches = plt.hist(temp_med_vec, bins=50,normed=1)
y = mlab.normpdf( bins,temp_med_vec.mean() , temp_med_std)
l = plt.plot(bins, y, 'k--', linewidth=1.5)
plt.xlim([T_min,T_max])
# plt.text(1.001*plt.xlim()[0],0.9*plt.ylim()[1],u"Média: "+str("%.2f" % temp_med_vec.mean())+u" °C\nDesvio Padrão: "+str("%.2f" % temp_med_vec.std())+" °C")
plt.figure()
plt.title(u"Umidade Relativa - MC Média de: "+str(tentativas)+" medidas")
plt.hist(umid_med_vec, bins=50)
# plt.title(u"Desvio Umidade Média para MC com "+str(tentativas)+" medidas: "+str(umid_med_std))
plt.title(u"Umid. MC "+str(tentativas)+" medidas: sigma = "+str("%0.3f" % umid_med_std)+u" | Média = "+str(umid_med))
n, bins, patches = plt.hist(umid_med_vec, bins=50,normed=1)
y = mlab.normpdf( bins,umid_med_vec.mean() , umid_med_std)
l = plt.plot(bins, y, 'k--', linewidth=1.5)
plt.xlim([U_ar_min,U_ar_max])
# plt.text(1.001*plt.xlim()[0],0.9*plt.ylim()[1],u"Média: "+str("%.2f" % umid_med_vec.mean())+u" %\nDesvio Padrão: "+str("%.2f" % umid_med_vec.std())+" %")
plt.figure()
plt.title(u"Pressão - MC Média de: "+str(tentativas)+" medidas")
plt.hist(pressao_med_vec, bins=50)
# plt.title(u"Desvio Pressão Média para MC com "+str(tentativas)+" medidas: "+str(pressao_med_std))
plt.title(u"Pressão MC "+str(tentativas)+" medidas: sigma = "+str("%0.3f" % pressao_med_std)+u" | Média = "+str(pressao_med))
n, bins, patches = plt.hist(pressao_med_vec, bins=50,normed=1)
y = mlab.normpdf( bins,pressao_med_vec.mean() , pressao_med_std)
l = plt.plot(bins, y, 'k--', linewidth=1.5)
plt.xlim([P_min,P_max])
# plt.text(1.001*plt.xlim()[0],0.9*plt.ylim()[1],u"Média: "+str("%.2f" % pressao_med_vec.mean())+u" Pa\nDesvio Padrão: "+str("%.2f" % pressao_med_vec.std())+" Pa")
plt.figure()
plt.title(u"Luminosidade - MC Média de: "+str(tentativas)+" medidas")
plt.hist(lum_med_vec, bins=50, normed=1)
# plt.title(u"Desvio lumin. Média para MC com "+str(tentativas)+" medidas: "+str(lum_med_std))
plt.title(u"Lum. MC "+str(tentativas)+" medidas: sigma = "+str("%0.3f" % lum_med_std)+u" | Média = "+str(lum_med))
n, bins, patches = plt.hist(lum_med_vec, bins=50, normed=1)
y = mlab.normpdf( bins,lum_med_vec.mean() , lum_med_std)
l = plt.plot(bins, y, 'k--', linewidth=1.5)
plt.xlim([L_min,L_max])
......@@ -279,7 +345,37 @@ for tentativas in [5,10,50]:
#plt.show()
#'''
#Apresenta os gráficos gerados
plt.show()
#print n_medias
#print temp_std_n
plt.figure()
plt.title(u"Desvio padrão para Temperatura média de N medidas")
plt.scatter(n_medias, temp_std_n)
for i,txt in enumerate(temp_std_n):
plt.annotate("%.3f" % txt, (n_medias[i]*1.01,temp_std_n[i]*1.01))
#for i, txt in enumerate(n):
# ax.annotate(txt, (z[i],y[i]))
plt.figure()
plt.title(u"Desvio padrão para Pressão média de N medidas")
plt.scatter(n_medias, pressao_std_n)
for i,txt in enumerate(pressao_std_n):
plt.annotate("%.3f" % txt, (n_medias[i]*1.01,pressao_std_n[i]*1.01))
plt.figure()
plt.title(u"Desvio padrão para Umidade Relativa média de N medidas")
plt.scatter(n_medias, umid_std_n)
for i,txt in enumerate(umid_std_n):
plt.annotate("%.3f" % txt, (n_medias[i]*1.01,umid_std_n[i]*1.01))
plt.figure()
plt.title(u"Desvio padrão para Luminosidade média de N medidas")
plt.scatter(n_medias, lum_std_n)
for i,txt in enumerate(lum_std_n):
plt.annotate("%.3f" % txt, (n_medias[i]*1.01,lum_std_n[i]*1.01))
#Apresenta os gráficos gerados
plt.show()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment