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

Cálculos estatísticos

Em seguida é calculada a média, o desvio padrão e a PDF da amostra. Nesta etapa também é determinado onde ocorre o máximo de probabilidade na PDF.

  dados_med = 0
  pdf1 = 0
  pdf2 = 0
  pdf3 = 0
  dados_med = mean(dados) #Calcula a média
  dados_desvpad = (variance(dados))**.5 #Calcula o desvio padrao

  y0, x0 = pdf(dados, kernel = k, n = npdf)
  y_t0, x_t0 = pdf(dados, kernel = k, n = npdf)
  y_int0, x_int0 = pdf(dados, kernel = k, n = npdf)
  y = array(y0)
  x = array(x0)
  y_t = array(y_t0)
  x_t = array(x_t0)
  y_int = array(y_int0)
  x_int = array(x_int0)
  gravar.write(nomearquivodados+'\n')
  testey = sort(y_t)
  probMax = float(testey[len(y_t)-1:len(y_t)])
  for i in range(0,len(x)):
    if y_t[i] == probMax:
      ponto_x = i #Posição da lista onde a probabilidade é máxima
      centro_x = x_t[i] #Valor da propriedade em analise onde a probabilidade é máxima
      string_MaxProb = "MaxProb = "+ str(probMax) + " @ " + str(x_t[i])
      print string_MaxProb
      gravar.write(string_MaxProb+'\n')
  integral = 0
  cont = 0
  somatorio_Prob=0.0
  lim_cont = 500 - ponto_x
  for i in range(0,len(x)):
    somatorio_Prob=somatorio_Prob+float(y[i]) #Calculo do somatorio para transformar a escala de probabilidade

Posteriormente é calculada a probabilidade acumulada, necessária para a delimitação da região correspondente a 95% de probabilidade a partir da máxima probabilidade calculada.

  PAC=[] #probabilidade acumulada
  xs=[] #Valor da propriedade para ser usado na simulacao
  for i in range(0,len(x_int)):
    xs.append(x[i])
    string_PDF = str(x[i]) + ',' +str(y[i])
    gravar_pdf.write(string_PDF+'\n')
    integral=integrate.trapz(y_int[0:i],x_int[0:i])
    PAC.append(integral)
    if integral/.025 < 1.0:# > 0.023 and integral < 0.025:
      ponto025 = i
    if integral/.975 <1.0:# > 0.975 and integral < 0.976:
      ponto975 = i

A penúltima etapa dos cálculos estatísticos consiste em determinar, a partir da PDF, cinco valores mais prováveis.

  dadosimulado=0.0
  prob=0.0
  nps=0 #número pontos simulados
  PAC.sort(reverse=True)
  xs.sort(reverse=True)
  #qp=5 #quantos pontos
  while nps<qp:
    for i in range(1,len(xs)/10):
      prob=rnd()
      if float(PAC[i])>prob:
        dadosimulado=float(xs[i-1])+(float(xs[i])-float(xs[i-1]))*(prob-PAC[i-1])/(PAC[i]-PAC[i-1])
        if dadosimulado>0:
          gravar_freq.write(str(dadosimulado)+'\n')
          nps=nps+1
        if nps>=qp:
          break

Então, os principais resultados são salvos.

  string_LInf = "Limite inferior = " + str(x_t[ponto025])
  print string_LInf
  gravar.write(string_LInf+'\n')
  string_LSup = "Limite superior = " + str(x_t[ponto975])
  gravar.write(string_LSup+'\n')
  print string_LSup
  string_Resist = "Propriedade (95%LC) = (-" + str(float(centro_x)-float(x_t[ponto025])) + ";" + str(float(centro_x)) + "; +" + str(float(x_t[ponto975])-float(centro_x)) +") "
  gravar.write(string_Resist+'\n')
  print string_Resist
  string_Area = "área integrada: "+str(integral)
  gravar.write(string_Area+'\n')
  print string_Area
  string_log = nomearquivodados+';'+str(x_t[ponto025])+';'+str(float(centro_x))+';'+str(x_t[ponto975])+';'+ str(dados_med) +';'+ str(dados_desvpad)+';'+str(integral)
  gravar_log.write(string_log+'\n')
  print somatorio_Prob

E, por último, o gráfico contendo o PDF é gerado.

  ax = subplot(111)
  quem = str(nomearquivodados)
  plot(x,y)
  grid()
  dados_cientifico= "%.5g" %(dados_med)
  titulo = "Sample: "+upper(quem) +"; Average value. = "+ dados_cientifico
 
  title(titulo)
  nome_graf1 = nomefig
  ylabel('Probability density',size=18)
  xlabel(abscissa, size = 18)
  savefig(nome_graf1)
  a = x_int[ponto025]
  b = x_int[ponto975]

  # Gera regiao sombreada
  ix = x_int[ponto025:ponto975]
  iy = y_int[ponto025:ponto975]
  verts = [(a,0)] + zip(ix,iy) + [(b,0)]
  poly = Polygon(verts, facecolor='0.8', edgecolor='k')
  ax.add_patch(poly)
  a=float('%.2f'%a)
  centro_x=float('%.2f'%centro_x)
  b=float('%.2f'%b)
  ax.set_xticks((a,centro_x,b))
  nome_graf2 = nomefig+"-P95"
  savefig(nome_graf2)
  close('all')

  gravar.close
  gravar_pdf.close
  gravar_log.close
  print"Processo concluído!"
  #show()

Resultados

Foram selecionados 20 valores entre 11 e 19, mostrados na Tabela 1, capazes de produzir uma distribuição bimodal.

Tabela 1: Vinte valores entre 11 e 19 selecionados para o teste
ContagemValor
111
212
313
414
513
612
711
813
918
1019
1117
1218
1317
1416
1511
1612
1713
1811
1912
2013

Os valores mostrados na Tabela 1 foram usados para construir um histograma, com o auxílio do software QtiPlot®, mostrado na Figura 2.

Probabilidades em Python: histograma
Figura 2: Histograma gerado a partir dos dados da Tabela 1.

Em seguida, os dados foram salvos num arquivo denominado teste.txt e fornecidos ao programa PDF&Freq. Como primeira prova, foram geradas PDFs usando diferentes kernels. Os resultados são mostrados na Figura 3.

Probabilidades em Python: PDFs obtidos usando os kernels Epanechnikov (a), Triangle (b), Gaussian (c) e Cosine (d)
Figura 3: PDFs obtidos usando os kernels Epanechnikov (a), Triangle (b), Gaussian (c) e Cosine (d)

Comprovada a influência da escolha do kernel sobre os resultados, os cálculos seguintes foram feitos usando o kernel Gaussiano. Os arquivos gerados pelo programa PDF&Freq são mostrados na Figura 4.

Probabilidades em Python: Arquivos gerados pelo programa PDF&Freq
Figura 4: Arquivos gerados pelo programa PDF&Freq

Entre estes arquivos, os que contêm as principais respostas geradas pelo programa são três:

  1. Arquivo log.dat, onde são listados o kernel escolhido e as demais informações previamente descritas.
    kernel= g
    Arquivo ; LI(95%) ; MAXPROB; LS(95%) ; Media ; DesvPad ; Area
    teste;8.980685;12.2418703;20.0478976;13.8;2.66754371;0.999612524
  2. Arquivo teste.frq, onde estão listados os cinco valores mais prováveis tomados para os valores aleatório da variável prob quando esta mais se aproxima de 1.
    9.40974719
    9.74840852
    9.27565556
    12.9108177
    19.9949159
  3. O arquivo teste.rep, onde são listadas a ordenada e a abscissa referente à máxima probabilidade, bem como os limites de confiança, com 95% de probabilidade.
    MaxProb = 0.152862861298 @ 12.2418703658
    Propriedade (95%LC) = (-3.26118472476;12.2418703658; +7.80602726672)

Assim, as principais informações estatísticas referentes à construção das PDFs são armazenadas numa única pasta, de maneira muito prática e rápida.

Conclusões

O maior impacto deste trabalho consiste na construção de uma ferramenta computacional que permite a rápida construção de distribuições de probabilidade a partir de diferentes kernels. Isso só foi possível com o uso da linguagem de programação Python, a qual permitiu desenvolver rapidamente um programa capaz de lidar com diversos arquivos de dados, produzindo informações estatísticas valiosas na forma de tabelas e gráficos, os quais são fundamentais para o nosso grupo de pesquisa e para os demais grupos lidem com esse tipo de informação ou problema.

Agradecimentos

Os autores agradecem ao Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), à Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), à Financiadora de Estudos e Projetos (FINEP PRESAL Ref.1889/10) e à Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) por toda a ajuda financeira e pelas bolsas concedidas.

Bibliografia

  1. Souza Jr., F. G., Pinto, J. C. & Soares, B. G. SBS/Pani – DBSA mixture plasticized with DOP and NCLS – Effect of the plasticizers on the probability density of volume resistivity measurements. Eur. Polym. J. 43, 2007–2016 (2007).
  2. Luce, B. & Anthony, O. A Primer on Bayesian Statistics in Health Economics and Outcomes Research. (MEDTAP International, Incorporated, 2003).
  3. Bard, Y. Nonlinear Parameter Estimation. (Academic Press, 1974).
  4. Grance, E. G. O. et al. New petroleum absorbers based on lignin-CNSL-formol magnetic nanocomposites. J. Appl. Polym. Sci. (2012).
  5. Araújo, A., Botelho, G., Oliveira, M. & Machado, A. V. Influence of clay organic modifier on the thermal-stability of PLA based nanocomposites. Appl. Clay Sci. 88–89, 144–150 (2014).
  6. de Souza Junior, F. G. et al. Conducting and magnetic mango fibers. Ind. Crops Prod. 68, 97–104 (2015).
  7. Souza Jr., F. G. & Varela, A. Construção de ferramenta de aquisição e inspeção de dados eletromecânicos usando Python. Revista PROGRAMAR 34, 32–38 (2012).
  8. Schwaab, M. Análise de Dados Experimentais: I. Fundamentos de Estatística e Estimação de Parâmetros. (Editora E-papers).
  9. Kernel Density Estimation in Python. Disponível em: http://jakevdp.github.io/blog/2013/12/01/kernel-density-estimation/. (Acedido: 19 Maio 2016)