11

Suponha a seguinte fórmula para calcular a variância:

variancia <- function(x) {
  n <- length(x)
  (1/(n^2-n))*(n*(sum(x^2))-(sum(x)^2))
}

Veja que ela é equivalente à função var na maior parte dos casos:

teste <- 1:5
var(teste)
[1] 2.5
variancia(teste)
[1] 2.5
all.equal(var(teste),variancia(teste))
[1] TRUE

Ou neste outro exemplo:

set.seed(1)
x1 <- rnorm(100, 10, 100)
var(x1)
[1] 8067.621
variancia(x1)
[1] 8067.621
all.equal(variancia(x1), var(x1))
[1] TRUE

Entretanto, no caso abaixo, ela resulta em um valor impossível (valor negativo):

set.seed(1)
x2 <- runif(1000) + 10^12
variancia(x2)
[1] -140878367
var(x2)
[1] 0.08316728

Por que a diferença entre as duas funções? Como assegurar que a função variancia obtenha o valor correto no último exemplo?

2
  • 10^12 não extrapola o maior inteiro representável em R?
    – user4552
    Commented 21/02/2014 às 22:10
  • No caso não são inteiros, são pontos flutuantes do tipo double. Eles têm 64 bits, 1 para o sinal, 11 para o expoente e 52 para o significante, o que dá uma precisão de aproximadamente 16 dígitos. Commented 21/02/2014 às 23:15

2 Respostas 2

11

A sua função foi vítima do cancelamento catastrófico. Isso pode acontecer quando se subtrai dois números próximos e de mesmo sinal, no caso da sua função:

sum(x2^2)
[1] 1e+27

sum(x2)^2 / length(x2)
[1] 1e+27

No caso da fórmula usada na função variancia isso ocorre normalmente quando a variância do vetor é muito menor que sua média.

Vou propor duas soluções pouco eficientes, mas simples:

  • Utilize outra fórmula:
variancia2 <- function(x) {
  n <- length(x)
  media <- mean(x)
  sum((x - media)^2) / (n - 1)
}

variancia2(x2)
[1] 0.08316728
  • Utilize a sua fórmula, mas retire antes a média do vetor, isso não altera o valor da variância.
variancia(x2 - mean(x2))
[1] 0.08316727
4

Complementando a resposta do Marcos Banik.

Um número do ponto flutuante do tipo double (64 bits) pode ser, grosso modo, resumido em 3 partes:

Ponto flutuante tipo double: sinal (1bit), ordem de magnitude (11 bits) e precisão (52 bits)

Isso dá para representar ordens de magnitude de cerca de 10^308 mas com uma precisão de cerca de 16 dígitos (detalhes sobre como o pacote base do R lida com números podem ser vistos na ajuda ?.Machine), além disso números irracionais ou cujo denominador não sejam potência de 2 são aproximados.

Então veja que um número muito grande pode ser representado pelo double, mas não com tanta precisão. Isto pode gerar grandes problemas com operações como soma e subtração. Os números calculados na fórmula variancia são da ordem de (10^12)^2=10^24 no terceiro exemplo, e só temos 52 bits para representar dígitos significantes (os demais são imprecisos). Quando subtraímos um do outro, eliminamos os "dígitos bons" e só sobram os "dígitos ruins", causando o resultado absurdo.

Uma forma de solucionar o problema é buscar algoritmos mais estáveis para ponto flutuante, como os que o Marcos propôs. Mas supondo que isso não seja possível, você pode utilizar números de precisão arbitrária.

No R o pacote Rmpfr (Multiple Precision Floating-Point Reliable) fornece números com precisão arbitrária (ao custo de gastar mais memória e tempo de execução, então dependendo do seu computador e do problema, nem sempre é possível).

Então, caso não fosse possível de alguma forma melhorar o algoritmo de cálculo da fórmula variancia, poderíamos utilizar o Rmpfr. Precisaríamos de mais de 30 dígitos de precisão, o que daria mais de log2(10^30)=99.65 bits. Arredondando para 128 bits:

library(Rmpfr)
x2.mpfr <- mpfr(x2, 128)
variancia(x2.mpfr)
1 'mpfr' number of precision  128   bits 
[1] 0.0831672741323709434253943823576340475867

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .