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}, \tag{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]. \tag{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 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}. \tag{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}, \tag{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} \tag{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 \tag{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 \tag{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} \tag{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} \tag{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, 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}, \tag{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 \tag{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)
res = medcov(x)
for keys,values in res.items():
    print(keys)
    print(values)
# comparar com resultado da nossa função
np.cov(x, rowvar = False, ddof = 1)
np.mean(x, axis = 0)
media
[ 0.03925061 -0.00917821  0.0590073   0.03687259  0.03115346]
covariância
[[ 0.99349167  0.01967448  0.03356703 -0.00719789  0.01570097]
 [ 0.01967448  1.02298624  0.0260564   0.00435293  0.02836326]
 [ 0.03356703  0.0260564   1.08074563 -0.00361641  0.007306  ]
 [-0.00719789  0.00435293 -0.00361641  0.94235813 -0.0237797 ]
 [ 0.01570097  0.02836326  0.007306   -0.0237797   1.09604901]]
SQP
[[ 992.49817754   19.65480478   33.53346475   -7.19068914   15.68526864]
 [  19.65480478 1021.96325757   26.03034316    4.34857569   28.33489577]
 [  33.53346475   26.03034316 1079.66488285   -3.61279364    7.29869703]
 [  -7.19068914    4.34857569   -3.61279364  941.41577529  -23.75592142]
 [  15.68526864   28.33489577    7.29869703  -23.75592142 1094.95296269]]
array([[ 0.99349167,  0.01967448,  0.03356703, -0.00719789,  0.01570097],
       [ 0.01967448,  1.02298624,  0.0260564 ,  0.00435293,  0.02836326],
       [ 0.03356703,  0.0260564 ,  1.08074563, -0.00361641,  0.007306  ],
       [-0.00719789,  0.00435293, -0.00361641,  0.94235813, -0.0237797 ],
       [ 0.01570097,  0.02836326,  0.007306  , -0.0237797 ,  1.09604901]])
array([ 0.03925061, -0.00917821,  0.0590073 ,  0.03687259,  0.03115346])

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 Estudos de Casos

Vamos considerar dois casos de instabilidades em computação numéricas com soluções clássicas, que são conhecidos por Log-Sum-Exp e Soma de Kahan. O truque Log-Sum-Exp (LSE) é usado principalmente para evitar overflow e underflow (de forma passiva) no cálculo de exponenciais, enquanto a soma de Kahan é usada para reduzir o erro de arredondamento acumulado na adição de números de ponto flutuante.

O truque LSE resolve a instabilidade causada por “faixa dinâmica excessiva”. Ao lidar com funções exponenciais, a diferença de magnitude entre as entradas pode ser enorme, fazendo com que os resultados intermediários excedam o intervalo representável dos números de ponto flutuante.

Já a soma de Kahan resolve a instabilidade causada pela “acumulação de erros”. Operações de ponto flutuante têm erros de arredondamento, e erros de cálculo únicos podem ser negligenciáveis, mas ao acumular através de milhares de operações, esses pequenos erros se amplificam, fazendo com que o resultado final perca completamente o significado. Algoritmos como o de Kahan mostram que, através de um projeto algorítmico cuidadoso, é possível controlar o limite de erro para que não cresça ilimitadamente com o tamanho dos dados.

5.4.1 LOg-Sum-Exp (LSE)

Sua função principal é calcular com segurança \(\log \left( \sum_{i=1}^{n} e^{x_i}\right)\). Se calculado diretamente, quando \(x_i\) é muito grande, \(e^{x_i}\) causará overflow (transbordamento), e quando \(x_i\) é muito pequeno (negativo), \(e^{x_i}\) causará underflow (subfluxo), resultando em zero.

O princípio central baseia-se na exploração da identidade matemática \[\log \left( \sum_{i=1}^{n} e^{x_i} \right)=m + \log \left( \sum_{i=1}^{n} e^{x_i-m}\right),\] onde \(m\) \(=\) \(\max(x_i)\). Isso equivale a um tipo de normalização dos valores, subtraindo o valor máximo de todos os expoentes, garantindo que o maior valor de entrada para a função exponencial seja \(0\), evitando assim completamente o overflow. O LSE garante que nenhum termo sofra overflow, e o underflow que ocorre atinge apenas termos cuja contribuição é numericamente desprezível.

Esse método é extremamente comum em modelos probabilísticos e algoritmos de aprendizado de máquina, especialmente ao trabalhar com probabilidades em espaço logarítmico. Por exemplo, na função Softmax para classificação em redes neurais, na função de perda de redes neurais, em cálculos de inferência bayesiana e no algoritmo de Viterbi (reconhecimento de fala, bioinformática, processamento de linguagem natural, comunicações digitais e GPS).

A solução estável, considerando que \(m\) \(=\) \(x_{\max}\) \(=\) \(\max(x_1, \dots, x_n)\), é obtida por \[ \begin{aligned} y =& \log \left( \sum_{i=1}^{n} e^{x_i} \right)= \log \left( \sum_{i=1}^{n} e^{x_i- m + m} \right)\\ =& \log \left(e^m \right)+ \log \left( \sum_{i=1}^{n} e^{x_i- m} \right)\\ =& m +\log \left( \sum_{i=1}^{n} e^{x_i- m} \right). \end{aligned} \]

5.4.2 Soma de Kahan (Soma Compensada)

O objetivo é reduzir o erro de arredondamento acumulado ao somar um grande número de números de ponto flutuante. Quando você soma um número muito pequeno a um número muito grande, os bits inferiores do número pequeno são “absorvidos” (perdidos), o que ao acumular gera grandes desvios.

O princípio central é manter uma variável de “compensação” adicional \(c\). Em cada iteração, primeiro subtrai-se a compensação da iteração anterior do valor a ser somado, depois realiza-se a soma normalmente, e em seguida calcula-se a parte de baixa ordem perdida nesta soma (ou seja, o novo valor de compensação) usando a diferença entre a soma e o valor original. Em termos simples, ele “recupera” e armazena temporariamente os pequenos valores que seriam perdidos, para “tentar novamente” na próxima adição.

O método é aplicado em cenários que exigem alta precisão de soma, como bibliotecas de computação científica, processamento de grandes conjuntos de dados (por exemplo, o CompensatedSum do Flink), e em algoritmos de simulação numérica que exigem resultados estáveis.

Imagine que queremos somar um elemento \(x\) ao nosso sum atual:

  1. \(y = x - c\)
    → Subtrai a compensação de erro acumulada do novo elemento \(x\).

  2. \(t = \text{sum} + y\)
    → Realiza a soma no acumulador principal. Aqui ocorre a perda de precisão por arredondamento.

  3. \(c = (t - \text{sum}) - y\)
    → Recupera o erro. Matematicamente, \((t - \text{sum})\) deveria ser igual a \(y\), mas numericamente, devido ao arredondamento, a diferença revela exatamente o que foi descartado. Guardamos isso em \(c\) com o sinal invertido para a próxima iteração.

  4. \(\text{sum} = t\)
    → Atualiza o acumulador principal.

Nota prática: O NumPy não usa a Soma de Kahan em np.mean, np.var, np.std ou np.cov. Em vez disso, utiliza pairwise summation (soma par a par / hierárquica), que reduz o erro de arredondamento de \(O(n)\) para \(O(\log n)\) e é quase tão rápida quanto a soma ingênua. Para a maioria dos casos estatísticos isso é suficiente. A Soma de Kahan oferece precisão ainda maior (erro próximo de \(O(1)\)), mas com um custo computacional um pouco mais alto.

5.4.3 Demonstração Prática em Python

O código abaixo executa ambos os cenários comparando a abordagem ingênua contra a abordagem numericamente estável.

import numpy as np

print("=== CASO 1: LOG-SUM-EXP ===")
# Valores altos que causam estouro de memória (overflow) no cálculo direto
x = np.array([1000.0, 1001.0])

# 1. Abordagem Ingênua
try:
    with np.errstate(over='ignore'):
        lse_ingenuo = np.log(np.sum(np.exp(x)))
except Exception as e:
    lse_ingenuo = "Erro no código"

# 2. Abordagem Estável (LSE)
x_max = np.max(x)
lse_estavel = x_max + np.log(np.sum(np.exp(x - x_max)))

print(f"LSE Ingênuo: {lse_ingenuo} (Resultado inválido por overflow de e^1000)")
print(f"LSE Estável: {lse_estavel}")

print("\n=== CASO 2: SOMA DE KAHAN ===")
# Criando 1.000.000 de elementos iguais a 0.1 usando float32
# O valor esperado exato da soma é 100.000,00
valores = np.ones(1_000_000, dtype=np.float32) * 0.1

# 1. Soma Tradicional
soma_tradicional = np.float32(0.0)
for v in valores:
    soma_tradicional += v

# 2. Soma de Kahan
soma_kahan = np.float32(0.0)
c = np.float32(0.0)  # Variável de compensação do erro
for v in valores:
    y = v - c            # Subtrai a compensação anterior do elemento atual
    t = soma_kahan + y   # Soma no acumulador principal (ocorre o arredondamento)
    c = (t - soma_kahan) - y  # Recupera a perda de precisão real
    soma_kahan = t       # Atualiza o acumulador

print(f"Soma Esperada:    100000.00")
print(f"Soma Tradicional: {soma_tradicional:.2f} (Errou por acúmulo de arredondamento)")
print(f"Soma de Kahan:    {soma_kahan:.2f} (Cravou o valor com precisão)")
=== CASO 1: LOG-SUM-EXP ===
LSE Ingênuo: inf (Resultado inválido por overflow de e^1000)
LSE Estável: 1001.3132616875182

=== CASO 2: SOMA DE KAHAN ===
Soma Esperada:    100000.00
Soma Tradicional: 100958.34 (Errou por acúmulo de arredondamento)
Soma de Kahan:    100000.00 (Cravou o valor com precisão)

5.5 Exercícios

  1. Mostrar a equivalência numérica existente entre as seguintes expressões \[\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}\\ &= W2_{k-1} + k\left(X_k-\bar{X}_{k}\right)^2/(k-1) \end{align*}\] e \[ W2_k = W2_{k-1} + (k-1)\left(X_k-\bar{X}_{k-1}\right)^2/k.\]

    Para isso implementar em Python as duas funções e realizar um comparativo numérico (um exemplo apenas).

  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.