3. Funções Matemáticas e Estatísticas

A linha de comando do R funciona como uma potente calculadora. Todas operações aritméticas e funções matemáticas principais estão disponíveis. Exemplo:

> 4 + 9
[1] 13
> 4 - 5
[1] -1
> 4 * 5
[1] 20
> 4 / 5
[1] 0.8
> 4^5
[1] 1024
>

A notação básica de operações algébricas, como a aplicação hierárquica de parênteses, também pode ser utilizada:

> (4 + 5 ) * 7 - (36/18)^3
[1] 55
> (2 * ( 2 * ( 2 * (3-4))))
[1] -8
>

Note que somente os parênteses podem ser utilizados nas expressões matemáticas. As chaves (“{}”) e os colchetes (“[]”) têm outras funções no R:

> (2 * { 2 * [ 2 * (3-4)]})
Error: syntax error in "(2 * { 2 * ["
>

Por que o R é uma potente calculadora? Experimente fazer a seguinte operação matemática na sua calculadora:

> 1 - (1 + 10^(-15))

<box left red | Exercício: Estimador de Pollard>

Pollard (1971) propôs o seguinte estimador para estimar a densidade no método de quadrantes:

<latex>

\widehat{N}_p = \frac{4(4n-1)}{\pi \sum_{i=1}^{n}\sum_{j=1}^4 r_{ij}^2}

</latex>

onde, <latex> r_{ij} </latex> é a distância de árvore do quadrante j no ponto i ao centro do ponto quadrante e n é o número de pontos quadrantes.

A variância desse estimador é:

<latex>

\textrm{Var}\{ \widehat{N}_p \} =  \frac{\widehat{N}_p}{4n-2}.

</latex>

Se o valor da soma do quadrado das distâncias for:

<latex>

\sum_{i=1}^{30}\sum_{j=1}^4 r_{ij}^2} = 2531.794

</latex>

qual a densidade estimada e sua variância?

</box>

<box left red | Exercíco: Área transversal de uma Árvore >

A área transversal de uma árvore é calculada assumindo que a secção transversal do tronco à altura do peito (1,3m) é perfeitamente circular. Se o diâmetro à altura do peito (DAP) de uma árvore for 13.5cm, qual a área transversal?

Se uma árvore possui três fustes com DAPs de: 7cm, 9cm e 12cm, qual a sua área transversal?

</box>

As funções matemáticas comuns também estão disponíveis e podem ser aplicadas diretamente na linha de comando:

> sqrt(9)   # Raiz Quadrada
[1] 3
> abs( - 1 )  # Módulo ou valor absoluto
[1] 1
> abs( 1 )
[1] 1
> log( 10 )   # Logaritmo natural ou neperiano
[1] 2.302585
> log( 10, base = 10) # Logaritmo base 10
[1] 1
> log10(10)   # Também logaritmo de base 10
[1] 1
> log( 10, base = 3.4076) # Logaritmo base 3.4076
[1] 1.878116
> exp( 1 )       # Exponencial
[1] 2.718282
>

As funções trigonométricas:

> sin(0.5*pi)    # Seno
[1] 1
> cos(2*pi)      # Coseno
[1] 1
> tan(pi)        # Tangente
[1] -1.224647e-16
>
> asin(1)       # Arco seno (em radianos)
[1] 1.570796
> asin(1) / pi * 180
[1] 90
>
> acos(0)      # Arco coseno (em radianos)
[1] 1.570796
> acos(0) / pi * 180
[1] 90
> atan(0)    # Arco tangente (em radianos)
[1] 0
> atan(0) / pi * 180
[1] 0
>

Funções para arredondamento:

> ceiling( 4.3478 )
[1] 5
> floor( 4.3478 )
[1] 4
> round( 4.3478 )
[1] 4
> round( 4.3478 , digits=3)
[1] 4.348
> round( 4.3478 , digits=2)
[1] 4.35
>

Funções matemáticas de especial interesse estatístico:

> factorial( 4 )         # Fatorial de 4
[1] 24
> choose(10, 3)          # Coeficientes binomiais: combinação de 10 3-a-3
[1] 120
>  

<box left red | Exercíco: Área transversal de uma Árvore (Revisitado) >

Se uma árvore possui três fustes com DAPs de: 7cm, 9cm e 12cm, qual o diâmetro (único) que é equivalente à sua área transversal?

</box>

<box left red | Exercício: Cálculo da Biomassa de Árvores do Cerrado>

O modelo alométrico de biomassa ajustado para árvores do Cerradão estabele que a biomassa é dada pela expressão:

<latex>

    \widehat{b_i} = \exp(-1.7953) \, d_i^{2.2974} 

</latex>

onde <latex> b_i </latex> é a biomassa em kg e <latex> d_i </latex> é o DAP em cm.

Já um outro modelo para biomassa das árvores na mesma situação tem a forma:

<latex>

      \widehat{\ln(b_i)} = -2.6464 + 1.9960 \ln(d_i) } + 0.7558 \ln(h_i) 

</latex>

onde <latex> h_i </latex> é a altura das árvores em m.

Para uma árvore com DAP de 15cm e altura de 12m, os modelos resultarão em estimativas muito distintas? </box>

O R também lida com operações matemáticas que envolvem elementos infinitos e elementos indeterminados:

> 1/0
[1] Inf
> -5/0
[1] -Inf
> 500000000000000000/Inf
[1] 0
> 0/0
[1] NaN
> Inf/Inf
[1] NaN
> log(0)
[1] -Inf
> exp(-Inf)
[1] 0
> sqrt(Inf)
[1] Inf
> sqrt( - 1 )
[1] NaN
Warning message:
NaNs produced in: sqrt(-1)
> 2 * NA
[1] NA
> 2 * NaN
[1] NaN
> NA / 10
[1] NA
> NaN / -1
[1] NaN
>

Note que determinadas palavras (além do nome das funções) estão reservadas no R, pois são utilizadas com significado especial:

• pi - constante <latex>\pi = 3.141593 </latex>;
• Inf - infinito;
• NaN - indeterminado (Not a Number), normalmente resultado de uma operação matemática indeterminada;
• NA - indeterminado (Not Avaiable), normalmente caracterizando uma observação perdida (missing value).

Na operações matemáticas, NaN e NA atuam sempre como indeterminado.

<box left red | Exercício Conceitual: O que é uma Observação Perdida>

Como se caracteriza uma observação perdida?

Quando o diâmetro de uma árvore deve ter o valor zero ou o valor NA?

E o peso de um animal? E a biomassa de uma florestal? E a espécie de uma ave?

</box>

Mais do que simples operações aritméticas, o R permite que executemos operações algébricas operando sobre variáveis prédefinidas.

Para definir uma variável, basta escolher um nome (lembre-se das regras de nomes no R) e atribuir a ela um valor:

>
> a = 3.6
> b = sqrt( 35 )
> c = -2.1
> a
[1] 3.6
> b
[1] 5.91608
> c
[1] -2.1
>
> a * b / c
[1] -10.14185
> b^c
[1] 0.02391820
> a + exp(c) - log(b)
[1] 1.944782
>  
> a - b * c / d
Error: object "d" not found
>          

Não esqueça de definir as variáveis previamente!!

<box left red | Exercício Conceitual: Criando Variáveis com Nomes Reservados >

O que acontece se você criar uma variável com o nome pi? Por exemplo,

> pi = 10

O que acontece com a constante <latex>\pi</latex>?

E se for criada uma constante de nome sqrt? O que acontece com a função raíz quadrada (sqrt())?

DICA: O que faz a função search, no comando:

> search()

</box>

Um detalhe importante em linguagens de programação é diferenciar um sinal de igualdade de um sinal de atribuição.

O sinal de igualdade faz uma comparação entre dois elementos. No R o sinal de igualdade é “==“ e será visto em detalhes quando a questão de comparações for tratada.

O R possui vários sinais de atribuição. No tópico acima utilizamos o sinal '=', mas o mesmo efeito pode ser obtido por dois outros meios:

• o sinal ”<-“ (composto pelo sinal de menor ”<“ seguido do hífen ”-“) atribui o valor à direita ao elemento à esquerda:

> k <- 2.46
> k
[1] 2.46

• o sinal ”->“ (composto pelo hífen ”-“ seguido do sinal de maior ”>“) atribui o valor à esquerda ao elemento à direita:

> 2.46 -> k2
> k2
[1] 2.46
>  

Essas três formas de atribuir valores a variáveis são equivalentes e cada usuário trabalha com a forma que julgar mais conveniente.

O R, e a linguagem S, foram criados para operar não apenas número-a-número como uma calculadora convencional.

O R é um ambiente vetorial, isto é, quase todas suas operações atua sobre uma seqüência de números, que genericamente chamaremos de vetores.

Para criar um vetor, podemos usar a função c (c = colar, concatenar). Essa função simplesmente cola todos os argumentos dados a ela, formando um vetor:

> a = c(1, 10, 3.4, pi, pi/4, exp(-1), log( 2.23 ), sin(pi/7) )
> a
[1]  1.0000000 10.0000000  3.4000000  3.1415927  0.7853982  0.3678794  0.8020016  0.4338837
> 

Para criar vetores de números com intervalo fixo unitário (intervalo de 1) se utiliza o operador seqüencial (:):

> b = 1:8
> b
[1] 1 2 3 4 5 6 7 8
> c = 20:32
> c
 [1] 20 21 22 23 24 25 26 27 28 29 30 31 32
> d = 2.5:10
> d
[1] 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5

Uma forma mais flexível de criar seqüências de números (inteiros ou reais) é usando a função 'seq':

> seq(10, 30)
 [1] 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
> seq(10, 30, by=2)
 [1] 10 12 14 16 18 20 22 24 26 28 30
> seq(1.5, 7.9, length=20)
 [1] 1.500000 1.836842 2.173684 2.510526 2.847368 3.184211 3.521053 3.857895
 [9] 4.194737 4.531579 4.868421 5.205263 5.542105 5.878947 6.215789 6.552632
[17] 6.889474 7.226316 7.563158 7.900000

Também é fácil criar uma seqüência de números repeditos utilizando a função 'rep':

> rep(5, 3)
[1] 5 5 5
> rep(1:5, 3)
 [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
> rep(1:5,each=3)
 [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5
>

<box left red | Exercício: Palmeira com Muitos Fustes I>

Uma palmeira perfilhada possui 10 fustes com os seguintes diâmetros: 5, 6, 7, 5, 10, 11, 6, 8, 9 e 7.

Crie um vetor 'dap' com os diâmetros acimas e uma seqüência que enumera os fustes.

</box>

Todas operações matemáticas aplicadas sobre um vetor, serão aplicadas sobre cada elemento desse vetor:

> 2 * a
[1]  2.0000000 20.0000000  6.8000000  6.2831853  1.5707963  0.7357589  1.6040032
[8]  0.8677675
> sqrt( a )
[1] 1.0000000 3.1622777 1.8439089 1.7724539 0.8862269 0.6065307 0.8955454
[8] 0.6586985
>
> log( a )
[1]  0.0000000  2.3025851  1.2237754  1.1447299 -0.2415645 -1.0000000 -0.2206447
[8] -0.8349787
>             

Se as variáveis que trabalhamos são vetores, operações matemáticas entre variáveis serão realizadas pareando os elementos dos vetores:

> a* b
[1]  1.000000 20.000000 10.200000 12.566371  3.926991  2.207277  5.614011
[8]  3.471070
> a - b
[1]  0.0000000  8.0000000  0.4000000 -0.8584073 -4.2146018 -5.6321206 -6.1979984
[8] -7.5661163
> a^(1/b)
[1] 1.0000000 3.1622777 1.5036946 1.3313354 0.9528356 0.8464817 0.9689709
[8] 0.9008898
>  
> sqrt( a )
[1] 1.0000000 3.1622777 1.8439089 1.7724539 0.8862269 0.6065307 0.8955454
[8] 0.6586985
> log( b )
[1] 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595 1.9459101
[8] 2.0794415
>             

Nesse tipo de operação é importante termos clareza sobre o que pretendemos quando fazemos uma operação matemática entre dois vetores de comprimentos (length) diferentes:

> length( a )
[1] 8
> length( b )
[1] 8
> length( c )
[1] 13
> a * c
 [1]  20.000000 210.000000  74.800000  72.256631  18.849556   9.196986
 [7]  20.852041  11.714861  28.000000 290.000000 102.000000  97.389372
[13]  25.132741
Warning message:
longer object length
        is not a multiple of shorter object length in: a * c
>     

<box left red | Exercício: Palmeira com Muitos Fustes II>

Uma palmeira perfilhada possui 10 fustes com os seguintes diâmetros: 5, 6, 7, 5, 10, 11, 6, 8, 9 e 7.

Calcule a área transversal de cada fuste dessa palmeira.

</box>

<box left red | Exercício: Bits e Bytes> Como construir uma seqüência que representa o aumento do número de bits por byte de computador, quando se dobra o tamanho dos bytes?

Essa seqüência numérica parte do 2 e dobra os valores a cada passo. </box>

Operações matemática sobre vetores, opera elemento-a-elemento do vetor. Já as funções estatísticas operam no vetor como um todo:

> mean( a )
[1] 2.491344
> var( b )
[1] 6
> max( c )
[1] 32
> sd( a )
[1] 3.259248
> sum( c )
[1] 338
> min( b )
[1] 1
> range( c )
[1] 20 32
>      

Algumas funções úteis que não são estatísticas, mas operam no vetor são:

> a
[1]  1.0000000 10.0000000  3.4000000  3.1415927  0.7853982  0.3678794  0.8020016
[8]  0.4338837
> sort(a)
[1]  0.3678794  0.4338837  0.7853982  0.8020016  1.0000000  3.1415927  3.4000000
[8] 10.0000000
> rev(sort(a))
[1] 10.0000000  3.4000000  3.1415927  1.0000000  0.8020016  0.7853982  0.4338837
[8]  0.3678794
> cumsum(sort(a))
[1]  0.3678794  0.8017632  1.5871613  2.3891629  3.3891629  6.5307556  9.9307556
[8] 19.9307556
> cumsum(a)
[1]  1.00000 11.00000 14.40000 17.54159 18.32699 18.69487 19.49687 19.93076
>
> diff(a)
[1]  9.0000000 -6.6000000 -0.2584073 -2.3561945 -0.4175187  0.4341221 -0.3681178
> diff( seq(10, 34, length=15) )
 [1] 1.714286 1.714286 1.714286 1.714286 1.714286 1.714286 1.714286 1.714286
 [9] 1.714286 1.714286 1.714286 1.714286 1.714286 1.714286
> 

<box left red | Exercício: Palmeira com Muitos Fustes III>

Uma palmeira perfilhada possui 10 fustes com os seguintes diâmetros: 5, 6, 7, 5, 10, 11, 6, 8, 9 e 7.

Encontre para essa palmeira as seguintes medidas:

1. a área transversal total (área basal) da palmeira (g);
2. o diâmetro médio dos fustes (<latex>\overline{d}</latex>);
3. o desvio padrão dos diâmetros (s);

Como podemos representar essa palmeira com vários fustes por uma única medida de diâmetro? Calculemos a área transversal a partir de três opções:

A) A partir do diâmetro total, ou soma dos diâmetros:

<latex>

         g_A = \frac{\pi}{4} \left( 10\,\overline{d} \right)^2

</latex>

B) A partir do diâmetro médio e o número de fustes:

<latex>

         g_B = 10\, \frac{\pi}{4} \left( \overline{d} \right)^2

</latex>

C) A partir do diâmetro médio e do desvio padrão:

<latex>

         g_C = \frac{\pi}{4} \left( (10)\overline{d}^2 + (10-1)s^2 \right)

</latex>

Quais dessas opções se aproxima do valor correto de área basal (g)? </box>

Já foi visto que ao se digitar o nome de uma função na linha de comando, o R retorna o código da função. Veja a diferença de:

> ls()

para:

> ls

A maioria das funções precisa de certas informações para orientar o seu procedimento, tais informações são chamados de argumentos.

Os argumentos de qualquer função são detalhadamente explicados nas páginas de ajuda sobre a função. Mas para uma rápida consulta dos argumentos de uma função podemos usar a função 'args':

> args(ls)
function (name, pos = -1, envir = as.environment(pos), all.names = FALSE,
    pattern)
NULL
> args(q)
function (save = "default", status = 0, runLast = TRUE)
NULL
> args(save.image)
function (file = ".RData", version = NULL, ascii = FALSE, compress = !ascii,
    safe = TRUE)
NULL
>    

Algumas funções, entretanto, são primitivas ou internas e seus argumentos não são apresentados. Geralmente, nesses casos os argumentos são bastante óbvios:

> args(sin)
NULL
> sin
.Primitive("sin")
>

Outras funções simplesmente não possuem argumentos:

> args(getwd)
function ()
NULL
> getwd
function ()
.Internal(getwd())
<environment: namespace:base>
>                                

Ao observar o resultado da função 'args', você notará que alguns argumentos são seguindos de uma expressão que se inicia com o sinal de igualdade ('='). A expressão após o sinal de igualdade é chamada de valor default do argumento. Se o usuário não informar o valor para um dado argumento, a função usa o valor default. Como exemplo veja a função 'save.image':

> args(save.image)
function (file = ".RData", version = NULL, ascii = FALSE, compress = !ascii,
    safe = TRUE)
NULL
> 

Se o usuário simplesmente evocar a função 'save.image()', sem informar o nome do arquivo onde á área de trabalho deve ser gravada, o R gravará as informações num arquivo com nome '.RData'.

<box left red | Exercício: Argumentos de Funções Estatísticas > Quais são os argumentos (e seus valores default) das seguintes funções:

• mean,
• sd,
• range, e
• cumsum.

</box>

<box left red | Exercício: Argumentos de Funções de Uso Comum > Quais são os argumentos (e seus valores default) das funções:

• sort,
• log, e
• seq.

O que é o argumento ”. . .”? </box>

Sendo um ambiente para análise de dados, o R dispõe de um grande conjunto de funções para trabalhar com Distribuições Estatísticas. Essas funções ajudam não só na análise de dados, como também permitem a simulação de dados.

A distribuição Normal é a distribuição central da teoria estatística. Para gerar uma amostra de observações de uma distribuição normal utilizamos a função 'rnorm':

> args( rnorm )
function (n, mean = 0, sd = 1)
NULL
> vn1 = rnorm( 1000, mean = 40,  sd = 9 )
> mean( vn1 )
[1] 39.47248
> sd( vn1 )
[1] 8.523735
> range( vn1 )
[1] 14.93126 62.11959
>
> vn2 = rnorm( 100000, mean = 40,  sd = 9 )
> mean( vn2 )
[1] 40.02547
> sd( vn2 )
[1] 9.025218
> range( vn2 )
[1]  3.40680 78.25496
>     

Se quisermos saber a probabilidade acumulada até um certo valor de uma variável com distribuição normal utilizamos a função 'pnorm':

> args(pnorm )
function (q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
NULL
>
> pnorm( 1.96, mean = 0 , sd = 1 )
[1] 0.9750021
> pnorm( 1.96 )
[1] 0.9750021
>
> pnorm( 27, mean = 20, sd = 7 )
[1] 0.8413447
> pnorm( 13, mean = 20, sd = 7 )
[1] 0.1586553
> 

Se quisermos obter o valor de um quantil da distribuição normal utilizamos a função 'qnorm':

> args( qnorm )
function (p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
NULL
> qnorm( 0.90 )
[1] 1.281552
> qnorm( 0.30 )
[1] -0.5244005
>
> qnorm( 0.90, 20, 7)
[1] 28.97086
> qnorm( 0.30, 20, 7)
[1] 16.32920
>  

A função 'dnorm' fornece a densidade probabilística para cada valor de uma variável Normal:

> args( dnorm )
function (x, mean = 0, sd = 1, log = FALSE)
NULL
> x = seq(-4, 4, length=10000)             # Seqüência de -4 a 4 com 10.000 valores
>
> plot(x, dnorm(x))                        # Curva da Dist. Normal com média 0 e desvio padrão 1
> points(x, dnorm(x, sd=2))                # Curva da Dist. Normal com média 0 e desvio padrão 2 (adicionada ao gráfico)
>

<box left red | Exercício: Amplitude Normal>

Tomando uma variável que segue a Distribuição Normal, o que acontence com a amplitude de variação dos dados à medida que o tamanho da amostra cresce (por exemplo n= 100, 1000, 10000)?

</box>

<box left red | Exercício: Intervalo Normal I > Qual o intervalo da Distribuição Normal Padronizada que têm a média no centro e contem 50% das observações? </box>

<box left red | Exercício: Intervalo Normal II > Qual a probabilidade de uma observação da variável Normal Padronizada estar no intervalo [-1.96 , 1.96]? </box>

O que foi apresentado para Distribuição Normal pode ser generalizado para todas as distribuições que o R trabalha.

Há quatro funções para se trabalhar com distribuições estatísticas:

• ddistrib - retorna a densidade probabilística para um dado valor da variável;
• pdistrib - retorna a probabilidade acumulada para um dado valor da variável;
• qdistrib - retorna o quantil para um dado valor de probabilidade acumulada;
• rdistrib - retorna valores (números aleatórios) gerados a partir da distribuição;

No caso da Distribuição Normal: distrib = norm. Para outras distribuições temos:

<box 35% orange centered | DISTRIBUIÇÕES ESTATÍSTICAS NO R>

</box>

<box left red | Exercício: Teste t > Você realizou um teste t de Student bilateral e obteve o valor t = 2.2 com 19 graus de liberdade.

O teste é significativo ao nível de probabilidade de 5%? E se o valor observado fosse t = 1.9? </box>

<box left red | Exercício: Teste F > Você realizou um teste F e obteve o valor F = 2.2 com 19 graus de liberdade no numerador e 24 graus de liberdade no denominador.

O teste é significativo ao nível de probabilidade de 5%? E se o valor observado fosse F = 2.5? </box>

<box left red | Exercício: Padrão Espacial I > Gere duas amostras (p.ex.: x e y) de tamanho 1000 (n=1000) de números da distribuição Uniforme.

Faça um gráfico plotando uma amostra contra a outra (plot(x,y)). Qual o padrão espacial observado?

Você consegue explicá-lo? </box>

<box left red | Exercício: Padrão Espacial II > Gere duas amostras (p.ex.: xp e yp) de tamanho 10 (n=10) de números da distribuição Uniforme.

Gere duas amostras (p.ex.: xf e yf) de tamanho 1000 (n=1000) de números da distribuição Normal com média zero e desvio padrão 0.2)

Faça um gráfico plotando a soma das amostras X (xp+xf) contra a soma das amostras Y (yp+yf) (plot(xp+xf,yp+yf)).

Qual o padrão espacial observado? Você consegue explicá-lo? </box>

<box left red red | Exercício: Gráfico Quantil-Quantil > Construa uma seqüência ordenada de 1000 números entre 0 e 1:

> p = seq(0, 1, length=1000)

O vetor 'p' representa um vetor de probabilidades acumuladas.

Gere 1000 números aleatórios da distribuição Normal Padronizada e coloque os números em ordem:

> x = sort( rnorm(1000) )

Faça um gráfico dos quantis da distribuição Normal, tomando o vetor 'p' de probabilidades, contra os valores de 'x':

> plot( qnorm(p), x )

Como é o gráfico resultante?

Repita o mesmo processo para as distribuições Exponencial ( 'rexp' ) e Lognormal ( 'rlnorm' ). Como são os gráficos resultantes? Por que? </box>

<box left red | Exercício: Padrão Temporal > Suponha que terremotos ocorram de forma completamente aleatória no decorrer do tempo.

Qual a melhor distribuição estatística para representar o intervalo de tempo entre dois terremotos sucessivos? </box>