Funções de distribuição de probabilidade e o Python

Este artigo aborda a construção de um programa escrito em Python capaz de calcular funções de distribuição de probabilidade (PDF) a partir de arquivos *.txt ou *.csv contendo uma ou duas colunas de dados. Essas informações estatísticas são importantes para compreensão de diversos problemas complexos que usualmente são representados, de forma equivocada, via uma abordagem Gaussiana simples. Assim, esperamos que a ferramenta aqui apresentada seja tão útil para toda a comunidade, como é para o nosso grupo.

Introdução

As funções de distribuição de probabilidade (PDF) são usadas para descrever a dispersão de dados experimentais em torno de uma média. As distribuições de probabilidade são geralmente representadas em termos de integrais como uma função de densidade de probabilidade e podem ser interpretadas como um histograma construído com intervalos infinitesimais. A PDF mais comum é a distribuição Gaussiana ou normal, geralmente utilizada para o cálculo dos limites de confiança de dados experimentais [1].

A construção de PDFs deve satisfazer as seguintes condições:

  1. PDFs são sempre positivas;
  2. A área delimitada pela PDF f(x) é sempre igual a 1.0 [2];
  3. A probabilidade de um valor aleatório para estar dentro do intervalo (α, Ω) é igual à área delimitada pela curva dentro deste intervalo, como mostrado na Figura 1.
Probabilidades em Python: PDF bimodal hipotética

Figura 1: PDF bimodal hipotética

A maioria das investigações científicas assume que as flutuações de dados experimentais podem ser representadas por uma função de densidade normal. Bard [3] já provou que essa abordagem é muito útil, principalmente devido à sua simplicidade. No entanto, as flutuações experimentais podem seguir diferentes funções de distribuição, que podem ser úteis para uma interpretação mais apurada dos dados experimentais.

Assim, buscando acelerar a análise de dados, nosso laboratório decidiu construir uma ferramenta em Python, capaz de gerar PDFs de forma rápida e confiável. O programa desenvolvido, denominado PDF&Freq, é capaz de gerar melhores informações estatísticas, as quais foram usadas em diversos trabalhos do grupo onde foi usado para gerar PDFs de tamanho de poros [4], diâmetro de partículas [5] e mesmo a resistividade volumétrica [6] de nossos materiais.

Método

Desenvolvimento e descrição do programa

O programa foi escrito usando a linguagem de programação Python. O Python foi lançado em 1991 pelo programador Guido van Rossum e é uma linguagem de alto nível [7]. Uma característica em destaque desta linguagem é a priorização da legibilidade do código, tornando-a mais concisa e clara. O corpo do programa foi construído usando um ambiente de desenvolvimento integrado denominado PYTHON-IDLE, o qual permitiu a rápida importação de módulos, além do fácil uso de algumas funções básicas. O programa foi dividido em vários trechos, buscando facilitar a sua compreensão.

Identificação do programa e listagem de atualizações

A primeira etapa do programa visa a identificação do mesmo para o usuário, fornecendo dados sobre a versão usada e a data de criação. Assim, foi utilizada a função print, que apresenta os seus objetos na tela do computador. Vários comentários foram acrescidos ao longo das linhas do programa usando o símbolo hash (#), o que permitiu manter um registro das atualizações.

#!/usr/bin/env python
# -*- coding: cp850 -*-

print "Programa para calculo de PDF's"
print "e para silulacao de valores mais provaveis"
print "Fernando Gomes - UFRJ"
print "V0.4 - 19/05/16"

#lista de atualizacoes:
#19/05/16 - Programa versao 0.4
#16/03/15 - Modifiquei a forma de verificar se ha uma ou duas linhas
#18/03/12 - Inseri rotina para detectar se o arquivo de entrada possui uma ou duas colunas
#16/06/10 - Programa versao 0.3
#21/05/10 - Programa versao 0.2
#18/03/10 - Programa versao 0.1

Importação dos módulos necessários

Os módulos numpy, matplotlib, os, pylab, random, scipy, statistics e string, necessários para o funcionamento do programa, foram invocados conforme descrito abaixo.

#Importacao dos modulos
from numpy import *
from numpy import log
from statistics import pdf
from matplotlib.patches import Polygon
from os import getcwd, listdir
from pylab import *
from random import random as rnd
from scipy import integrate
from statistics import pdf
from statistics import mean
from statistics import variance
import string
from string import upper

Escolha do kernel

Os histogramas são muito úteis para o tratamento de dados discretos e remetem intuitivamente, a alguma forma de estimativa de densidade. A heurística diz que é possível construir distribuições contínuas, desde que estejam disponíveis mais de 20 experiências [8]. Para isso, a escolha de um kernel capaz de descrever as nuances amostrais da amostra é fundamental. Existem várias opções disponíveis para o cálculo de estimativas de densidade usando o Python. Cada uma delas depende de quais são os objetivos particulares da análise [9]. Vários tipos de funções do kernel são disponíveis. Os mais comuns, como o Epanechnikov, Uniform, Triangle, Gaussian, Quartic/biweight, Triweight e Cosine estão disponíveis no programa.

print
print
print "Escolha o kernel:"
print
print "-'E' or 'Epanechnikov' : Epanechnikov kernel (padrao)"
print "-'U' or 'Uniform' : Uniform kernel"
print "-'T' or 'Triangle' : Triangle kernel"
print "-'G' or 'Gaussian' : Gaussian kernel"
print "-'B' or 'Biweight' : Quartic/biweight kernel"
print "-'3' or 'Triweight' : Triweight kernel"
print "-'C' or 'Cosine' : Cosine kernel"
print
k = raw_input("Escolha o kernel (E, U, T, G, B, 3 ou C)")
print
if k == "":
  k = 'G'
  print "Vc escolheu a opcao padrao (G)"

Recolha de informações gerais

Escolhido o kernel, o utilizador deve inserir a legenda desejada na abscissa do gráfico PDF, conforme descrito abaixo.

print
print
print
print "Voce deve escolher a legenda da abscissa"
print "Ex.: Diameter (nm)"
abscissa=raw_input("Qual eh a legenda da abscissa? ")
print
print
print
print "A legenda da abscissa eh:",abscissa
print
print
correto=raw_input("A legenda esta correta? (S/N) ")
while correto=="n" or correto=="N":
  abscissa=raw_input("Qual eh a legenda correta?")
  print
  print abscissa
  correto=raw_input("A legenda esta correta? (S/N) ")
  print

Além disso, o utilizador pode optar por gerar diferentes quantidades de valores simulados mais prováveis. A opção padrão é de 5 números.

print
print "Simulacao dos valores mais provaveis"
print "Padrao = 5"
print
qp = raw_input("Quantos pontos? ")
print
if qp=="":
  qp=5
else:
  qp=int(qp)

gravar_log=open('log.dat','w') #Cria arquivo de log
gravar_log.write("kernel= "+k+'\n')
string_log = 'Arquivo ; LI(95%) ; MAXPROB; LS(95%) ; Media ; DesvPad ; Area'
gravar_log.write(string_log+'\n')
gravar_log.close()

Como observação, destacamos que o número de pontos para a construção da PDF, por padrão, é igual a 500.

npdf = 500 #Numero de pontos para a contrucao da PDF

Criação do arquivo de log

Um arquivo de log é criado para armazenar todas as informações relevantes como, o nome da amostra (Arquivo), o limite de confiança inferior com 95% de probabilidade (LI95%), a probabilidade máxima (MAXPROB), o limite de confiança superior com 95% de probabilidade (LS95%), o valor médio (Media), o desvio padrão (DesvPad) e a área sob a curva PDF (Area).

gravar_log=open('log.dat','w') #Cria arquivo de log
gravar_log.write("kernel= "+k+'\n')
string_log = 'Arquivo ; LI(95%) ; MAXPROB; LS(95%) ; Media ; DesvPad ; Area'
gravar_log.write(string_log+'\n')
gravar_log.close()

Identificação dos arquivos para análise

Em seguida, são selecionados todos os arquivos *.txt ou *.csv presentes na pasta. Esses arquivos são listados e depois abertos, um a um, para a análise PDF, conforme descrito abaixo.

Conj_arquivos=listdir(getcwd())
arquivos=[]
for i in range(0,len(Conj_arquivos)):
  if Conj_arquivos[i][-3]=='t' and Conj_arquivos[i][-2]=='x' and Conj_arquivos[i][-3]=='t':
    arquivos.append(Conj_arquivos[i][0:-4])
  if Conj_arquivos[i][-3]=='c' and Conj_arquivos[i][-2]=='s' and Conj_arquivos[i][-3]=='v':
    arquivos.append(Conj_arquivos[i][0:-4])

print type(arquivos)
print arquivos
for i in range(0, len(arquivos)):
  print str(arquivos[i])
  nomearquivodados = string.rstrip(str(arquivos[i]),'\n')
  nomearquivodados = string.rstrip(str(arquivos[i]),'.txt')
  nomearquivodados = string.rstrip(str(arquivos[i]),'.csv')
  nomefig = nomearquivodados
  print nomearquivodados
  try:
    arquivo=open(nomearquivodados+".csv", "r") #Abre arquivo para leitura
  except:
    arquivo=open(nomearquivodados+".txt", "r") #Abre arquivo para leitura
  gravar=open(nomearquivodados+".rep", "w") #Report
  gravar_pdf=open(nomearquivodados+".fdp", "w")  #Distribui‡Æo de probabilidade
  gravar_freq=open(nomearquivodados+".frq", "w") #Arquivo com a frequencia das repeticoes
  gravar_log=open('log.dat','a') #Arquivo de log
  dados=[]
  for line in arquivo:
    linha_extra=line.split(',')
    if len(linha_extra)>1: #Verifica se ha uma coluna ou duas
      linha_extra[1].rstrip('\r\n')
      dados.append(float(linha_extra[1])) #usar dados da segunda coluna
    else:
      dados.append(float(line)) #para o caso do arquivo ter apenas uma coluna