5  Algoritmos para Médias, Variâncias e Covariâncias

Muitos algoritmos para o cálculo de médias, variâncias e covariâncias são imprecisos, podendo gerar resultados finais contendo grandes erros. A necessidade de utilizarmos algoritmos eficientes para realizarmos estas tarefas simples são evidentes e serão descritos neste capítulo.

5.1 Introdução

Felizmente o Python utiliza algoritmos precisos para cálculo da média, da variância e de covariância. Vamos buscar esclarecer como a utilização de algoritmo ineficientes podem levar a resultados inconsistentes e imprecisos. Nosso objetivo neste capítulo é apresentar os algoritmos eficientes para estimarmos estes parâmetros quando as fórmulas convencionais podem falhar se utilizadas nos algoritmos diretamente.

Estes algoritmos eficientes são particularmente úteis quando os dados possuem grande magnitude ou estão muito próximos de zero. Neste caso particular, algumas planilhas eletrônicas, como o Excel nas suas versões mais antigas, podiam falhar (McCullough and Wilson 1999). O conhecimento de algoritmos que conduzirão a maiores precisões numéricas pode levar o pesquisador a não cometer os mesmos erros encontrados em alguns softwares.

5.2 Algoritmos Univariados

A ideia básica de utilizarmos algoritmos para média, para a soma de potências dos desvios em relação a média ou para soma de produtos de desvios é aplicarmos recursivamente as fórmulas existentes. Se desejamos, por exemplo, obter a média de \(n\) observações, podemos inicialmente calcular a média da primeira observação, que é a própria. Em um segundo estágio, propor uma expressão para atualizarmos a média da primeira observação, contemplando a segunda e assim sucessivamente. O mesmo procedimento de atualização é aplicado às variâncias ou às covariâncias, no caso de termos mais de uma variável.

Para uma amostra de tamanho \(n\) dada por \(X_1\), \(X_2\), \(\cdots\), \(X_n\), a média e a variância amostral convencionais são obtidas, respectivamente, por:

\[ \bar{X}_. = \dfrac{\displaystyle \sum_{i=1}^n X_i}{n}, \qquad(5.1)\] e \[ S^2 = \dfrac{1}{n-1}\left[\sum_{i=1}^n X_i^2 - \dfrac{\left(\displaystyle \sum_{i=1}^n X_i\right)^2}{n}\right]. \qquad(5.2)\]

Alguns algoritmos existentes procuraram melhorar os algoritmos dos livros textos que são as fórmulas das equações (Equation 5.1 e Equation 5.2) procurando fazer adaptações, repassando a amostra duas vezes. Estes algoritmos podem ser eficientes do ponto de vista da precisão, mas não são rápidos, justamente por repassarem duas vezes os dados. West (1979) propôs utilizar um algoritmo que faz uma única passada, atualizando a média e a variância em cada nova observação. Podemos facilmente mostrar que para uma amostra de tamanho \(n\), que a média é igual a \(X_1\) se \(n=1\), \((X_1+X_2)/2\) se \(n=2\) e assim por diante. No \((k-1)\)-ésimo passo podemos especificar o estimador da média por: \[\begin{align*} \bar{X}_{k-1} = & \frac{\displaystyle \sum_{i=1}^{k-1} X_i}{k-1}. \end{align*}\]

No k-ésimo passo teremos observado \(X_k\) e a média atualizada é:

\[ \bar{X}_{k} = \frac{\displaystyle \sum_{i=1}^{k} X_i}{k}. \qquad(5.3)\]

A pergunta que fazemos é “podemos expressar a média do \(k\)-ésimo passo em função da média do \((k-1)\)-ésimo passo?” A resposta a esta pergunta é sim e o resultado nos fornece o algoritmo desejado. A partir da equação (Equation 5.3) obtemos: \[\begin{align*} \bar{X}_{k} = & \frac{\displaystyle \sum_{i=1}^{k-1} X_i + X_k}{k}\\ = & \frac{\displaystyle (k-1)\sum_{i=1}^{k-1} X_i}{(k-1)k}+\dfrac{X_k}{k}\\ = & \frac{\displaystyle (k-1)\bar{X}_{k-1}}{k}+\dfrac{X_k}{k}\\ = & \bar{X}_{k-1} -\frac{\bar{X}_{k-1}}{k}+\dfrac{X_k}{k} \end{align*}\] resultando na equação recursiva final \[ \bar{X}_{k} = \bar{X}_{k-1} +\dfrac{X_k-\bar{X}_{k-1}}{k}, \qquad(5.4)\] para \(2\le k\le n\), sendo que \(\bar{X}_1=X_1\).

Da mesma forma se definirmos a soma de quadrados corrigidas das \(k\) primeiras observações amostrais \(1<k\le n\) por: \[\begin{align*} W2_k = & \sum_{i=1}^k X_i^2 - \dfrac{\left(\displaystyle \sum_{i=1}^k X_i\right)^2}{k} \end{align*}\] veremos que a variância correspondente é dada por \(S^2_k\) \(=\) \(W2_k/(k-1)\). Se expandirmos esta expressão isolando o \(k\)-ésimo termo e simplificarmos a expressão resultante teremos:

\[ \begin{aligned} W2_k = \sum_{i=1}^{k-1} X_i^2 +X_{k}^2- \dfrac{\left(\displaystyle\sum_{i=1}^{k-1} X_i +X_k\right)^2}{k}\\ = W2_{k-1} + k\left(X_k-\bar{X}_{k}\right)^2/(k-1) \end{aligned} \qquad(5.5)\]

A expressão que desenvolvemos (Equation 5.5) é equivalente a apresentada por West (1979). Isso pode ser demonstrado facilmente se substituirmos \(\bar{X}_k\) obtida na equação (Equation 5.4) na equação (Equation 5.5), de onde obtivemos:

\[ W2_k = W2_{k-1} + (k-1)\left(X_k-\bar{X}_{k-1}\right)^2/k \qquad(5.6)\] para \(2\le k\le n\), sendo que \(W2_{1}=0\). A variância é obtida por \(S^2 = W2_n/(n-1)\).

Para demonstrarmos diretamente a expressão (Equation 5.6) temos: \[\begin{align*} W2_k = & \sum_{i=1}^{k-1} X_i^2 +X_{k}^2- \dfrac{\left(\displaystyle\sum_{i=1}^{k-1} X_i +X_k\right)^2}{k}\\ = & \sum_{i=1}^{k-1} X_i^2 +X_k^2 - \dfrac{\left(\displaystyle\sum_{i=1}^{k-1} X_i\right)^2 + 2X_k\displaystyle\sum_{i=1}^{k-1} X_i+X_k^2}{k}\\ = & \sum_{i=1}^{k-1} X_i^2 +X_k^2 - \dfrac{(k-1)\left(\displaystyle\sum_{i=1}^{k-1} X_i\right)^2 + 2(k-1)X_k\displaystyle\sum_{i=1}^{k-1} X_i+ (k-1) X_k^2}{k(k-1)}\\ = & \sum_{i=1}^{k-1} X_i^2 -\dfrac{\left(\displaystyle\sum_{i=1}^{k-1} X_i\right)^2}{k-1} + X_k^2 + \dfrac{\left(\displaystyle\sum_{i=1}^{k-1} X_i\right)^2}{k(k-1)} -\dfrac{2(k-1) X_k \displaystyle\sum_{i=1}^{k-1} X_i}{k(k-1)} -\dfrac{X_k^2}{k}\\ = & W2_{k-1} + \dfrac{\left(\displaystyle\sum_{i=1}^{k-1} X_i\right)^2}{k(k-1)} -\dfrac{2(k-1) X_k \displaystyle\sum_{i=1}^{k-1} X_i}{k(k-1)} +\dfrac{(k-1) X_k^2}{k}\\ = & W2_{k-1} + \dfrac{(k-1)\bar{X}_{k-1}^2}{k} -\dfrac{2(k-1) X_k \bar{X}_{k-1}}{k} +\dfrac{(k-1) X_k^2}{k}, \end{align*}\] resultando em \[\begin{align*} W2_k = W2_{k-1} + (k-1)\left(X_k-\bar{X}_{k-1}\right)^2/k. \end{align*}\]

Podemos generalizar essa expressão para computarmos a covariância \(S_{(x,y)}\) entre uma variável \(X\) e outra \(Y\). A expressão para a soma de produtos é dada por:

\[ W2_{k,(x,y)} = W2_{k-1,(x,y)} + (k-1)\left(X_k-\bar{X}_{k-1}\right)\left(Y_k-\bar{Y}_{k-1}\right)/k \qquad(5.7)\] para \(2\le k\le n\), sendo \(W2_{1,(x,y)}=0\). O estimador da covariância é obtido por \(S_{(x,y)} = W2_{n,(x,y)}/(n-1)\).

Podemos observar, analisando todas estas expressões, que para efetuarmos os cálculos da soma de quadrados e da soma de produtos corrigida necessitamos do cálculo das médias das variáveis no \(k\)-ésimo passo ou no passo anterior. A vantagem é que em uma única passagem pelos dados, obtemos todas as estimativas. Podemos ainda estender estes resultados para obtermos somas das terceira e quarta potências dos desvios em relação a média. As expressões que derivamos para isso são:

\[ W3_k = W3_{k-1}+\dfrac{(k^2-3k+2)(X_k-\bar{X}_{k-1})^3}{k^2}-\dfrac{3(X_k-\bar{X}_{k-1})W2_{k-1}}{k} \qquad(5.8)\] e \[ \begin{aligned} W4_k =& W4_{k-1}+\dfrac{(k^3-4k^2+6k-3)(X_k-\bar{X}_{k-1})^4}{k^3}+\\ & + \dfrac{6(X_k-\bar{X}_{k-1})^2W2_{k-1}}{k^2} -\dfrac{4(X_k-\bar{X}_{k-1})W3_{k-1}}{k} \end{aligned} \qquad(5.9)\] para \(2\le k\le n\), sendo \(W3_{1}=0\) e \(W4_{1}=0\).

Desta forma podemos implementar a função medsq() que retorna a média, a soma de quadrado, cubo e quarta potência dos desvios em relação a média e variância a partir de um vetor de dados \(\mathbf{X}\) de dimensão \(n\).

# função para retornar a média, somas de desvios em relação
#  a média # ao quadrado, ao cubo e quarta potência e variância
def medsq(x):
  n = len(x)
  if (n <= 1):
    print('Vetor deve ter mais de 1 elemento!')
    return
  xb = x[0]  
  w2 = 0
  w3 = 0
  w4 = 0
  for j in range(1, n):
    d = x[j] - xb
    i = j + 1.0
    w4  = w4 + (i**3 - 4 * i**2 + 6 * i - 3) * d**4 / i**3 + \
      6 * w2 * d**2 / i**2 - 4 * w3 * d / i
    w3 = w3 + (i**2 - 3 * i + 2) * d**3 / i**2 - 3 * w2 * d / i
    w2 = w2 + (i - 1) * d**2 / i  
    xb = xb + d / i
  s2 = w2 / (n - 1)
  res = {'media': xb, 'variância': s2, 'SQ2': w2, 'W3': w3, 'W4': w4}
  return res
# Exemplo
x = [1, 2, 3, 4, 5, 7, 8]
res = medsq(x)
for key in res.keys():
    print(key, ":", res[key])
media : 4.285714285714286
variância : 6.571428571428572
SQ2 : 39.42857142857143
W3 : 22.04081632653064
W4 : 391.4518950437319

5.3 Algoritmos para Vetores Médias e Matrizes de Covariâncias

Vamos apresentar nesta seção a extensão multivariada para obtermos o vetor de médias e as matrizes de somas de quadrados e produtos e de covariâncias. Por essa razão não implementamos uma função específica para obtermos a covariância entre duas variáveis \(X\) e \(Y\). Seja uma amostra aleatória no espaço \(\mathbb{R}^p\) dada por \(\mathbf{X}_1\), \(\mathbf{X}_2\), \(\cdots\), \(\mathbf{X}_j\), \(\cdots\), \(\mathbf{X}_n\), sendo que estes vetores serão dispostos em uma matriz \(\mathbf{X}\) de dimensões \((n \times p)\). Para estendermos os resultados da seção anterior, utilizaremos as mesmas expressões, tomando o devido cuidado para adaptá-las para lidar com operações matriciais e vetoriais. Assim, para estimarmos o vetor de médias populacionais, devemos em vez de utilizar o estimador clássico dos livros textos dado por \[\begin{align*} \mathbf{\bar{X}}= & \dfrac{\displaystyle\sum_{i=1}^n \mathbf{X}_i}{n}, \end{align*}\] utilizaremos a expressão recursiva dada por \[ \mathbf{\bar{X}}_{k} = \mathbf{\bar{X}}_{k-1} +\dfrac{\mathbf{X}_k-\mathbf{\bar{X}}_{k-1}}{k}, \qquad(5.10)\] para \(2\le k\le n\), sendo que \(\mathbf{\bar{X}}_1=\mathbf{X}_1\).

Da mesma forma, a adaptação para \(p\) dimensões da expressão (Equation 5.6) é direta e o resultado obtido é:

\[ \mathbf{W}_{k} = \mathbf{W}_{k-1} + (k-1)\left(\mathbf{X}_k-\mathbf{\bar{X}}_{k-1}\right)\left(\mathbf{X}_k-\mathbf{\bar{X}}_{k-1}\right)^\top/k \qquad(5.11)\] para \(2\le k\le n\), sendo \(\mathbf{W}_{1}=\mathbf{0}\), uma matriz de zeros de dimensões \((p \times p)\). O estimador da matriz de covariâncias é obtido por \(\mathbf{S} = \mathbf{W}_{n}/(n-1)\).

Implementamos a função medcov() apresentada a seguir para obtermos o vetor de médias, a matriz de somas de quadrados e produtos e a matriz de covariâncias. O argumento desta função deve ser uma matriz de dados multivariados com \(n\) linhas (observações) e \(p\) colunas (variáveis). O programa resultante e um exemplo são apresentados na sequência. Escolhemos um exemplo onde geramos uma matriz de dados de uma normal multivariada.

# função para retornar o vetor de médias, a matriz de 
# somas de # quadrados e produtos e a matriz de covariâncias
import numpy as np
def medcov(x):
  n = x.shape[0]
  p = x.shape[1]
  if (n <= 1):
    print('Matriz deve ter mais de 1 linha!')
    return
  xb = x[0,:]  
  w  = np.zeros((p, p))
  for j in range(1, n):
    d = x[j,:] - xb
    i = j + 1.0
    w = w + (i - 1) * np.outer(d, d)  / i
    xb = xb + d / i
  s = w / (n - 1)
  res = {'media': xb, 'covariância': s, 'SQP': w}
  return res
# Exemplo
n = 1000
p = 5
x = np.random.multivariate_normal(np.zeros(p),np.identity(p),n)
medcov(x)
# comparar com resultado da nossa função
np.cov(x, rowvar = False)
np.mean(x, axis = 0)
{'media': array([-0.06713584, -0.01274922, -0.07107321, -0.02224068,  0.00835804]),
 'covariância': array([[ 1.05888790e+00,  1.94848414e-02, -4.15627698e-02,
         -6.51670629e-02, -1.22762016e-02],
        [ 1.94848414e-02,  9.91799418e-01, -8.29786121e-04,
         -1.58557832e-02,  1.55196004e-02],
        [-4.15627698e-02, -8.29786121e-04,  1.00532576e+00,
          1.16397452e-02,  3.06002146e-02],
        [-6.51670629e-02, -1.58557832e-02,  1.16397452e-02,
          1.01894586e+00,  1.80459064e-02],
        [-1.22762016e-02,  1.55196004e-02,  3.06002146e-02,
          1.80459064e-02,  1.10124964e+00]]),
 'SQP': array([[ 1.05782901e+03,  1.94653566e+01, -4.15212071e+01,
         -6.51018958e+01, -1.22639254e+01],
        [ 1.94653566e+01,  9.90807618e+02, -8.28956335e-01,
         -1.58399274e+01,  1.55040808e+01],
        [-4.15212071e+01, -8.28956335e-01,  1.00432043e+03,
          1.16281054e+01,  3.05696144e+01],
        [-6.51018958e+01, -1.58399274e+01,  1.16281054e+01,
          1.01792692e+03,  1.80278605e+01],
        [-1.22639254e+01,  1.55040808e+01,  3.05696144e+01,
          1.80278605e+01,  1.10014839e+03]])}
array([[ 1.05888790e+00,  1.94848414e-02, -4.15627698e-02,
        -6.51670629e-02, -1.22762016e-02],
       [ 1.94848414e-02,  9.91799418e-01, -8.29786121e-04,
        -1.58557832e-02,  1.55196004e-02],
       [-4.15627698e-02, -8.29786121e-04,  1.00532576e+00,
         1.16397452e-02,  3.06002146e-02],
       [-6.51670629e-02, -1.58557832e-02,  1.16397452e-02,
         1.01894586e+00,  1.80459064e-02],
       [-1.22762016e-02,  1.55196004e-02,  3.06002146e-02,
         1.80459064e-02,  1.10124964e+00]])
array([-0.06713584, -0.01274922, -0.07107321, -0.02224068,  0.00835804])

Felizmente, devido ao Python (numpy) usar precisão dupla e utilizar algoritmos de ótima qualidade para estas tarefas, não precisaremos nos preocupar com a implementação de funções como as apresentadas neste capítulo. Mas se formos utilizar um compilador da linguagem Pascal, Fortran ou C e C\(^{++}\), deveremos fazer uso destes algoritmos, pois somente assim alcançaremos elevada precisão, principalmente se tivermos lidando com dados de grande magnitude ou muito próximos de zero.

5.4 Exercícios

  1. Mostrar a equivalência existente entre as expressões (Equation 5.5) e (Equation 5.6).

  2. Implementar em Python uma função para obtermos a soma de produtos, equação (Equation 5.7), e covariância entre as \(n\) observações de duas variáveis e cujos valores estão dispostos em dois vetores \(\mathbf{X}\) e \(\mathbf{Y}\). Criar uma matriz com \(n\) linhas e \(p=2\) colunas de algum exemplo e utilizar a função medcov() para comparar os resultados obtidos.

  3. Os coeficientes de assimetria e curtose amostrais são funções das somas de potências dos desvios em relação a média. O coeficiente de assimetria é dado por \(\sqrt{b_1}\) \(=\) \((W3_n/n)/(W2_n/n)^{3/2}\) e o coeficiente de curtose é \(b_2 = (W4_n/n)/(W2_n/n)^2\). Implementar uma função Python, utilizando a função medsq() para estimar o coeficiente de assimetria e curtose univariados.

  4. Utilizar uma amostra em uma área de sua especialidade e determinar a média, variância, soma de quadrados, soma de desvios ao cubo e na quarta potência. Determinar também os coeficientes de assimetria e curtose.