6 Miscelânia de funcionalidades do R

6.1 O R como calculadora

Podemos fazer algumas operações matemáticas simples utilizando o R. Vejamos alguns exemplos calculando as seguintes somas:

(a)

\(10^2 + 11^2 + \ldots + 20^2\)
Para obter a resposta devemos

E estes passos correspondem aos seguintes comandos

> (10:20)

 [1] 10 11 12 13 14 15 16 17 18 19 20

> (10:20)^2

 [1] 100 121 144 169 196 225 256 289 324 361 400

> sum((10:20)^2)

[1] 2585

Note que só precisamos do último comando para obter a resposta, mas é sempre útil entender os comandos passo a passo!

(b)

\(\sqrt {\log (1)} + \sqrt {\log (10)} + \sqrt {\log (100)} + \ldots + \sqrt {\log (1000000)} \),
onde \(\log \) é o logaritmo neperiano. Agora vamos resolver com apenas um comando:

> sum(sqrt(log(10^(0:6))))

[1] 16.4365

6.2 Gráficos de funções

Para ilustrar como podemos fazer gráficos de funções vamos considerar cada uma das funções a seguir cujos gráficos são mostrados nas Figuras 12.

(a)

\(f(x) = 1 - \frac {1}{x} \sin (x)\) para \(0 \leq x \leq 50\)

(b)

\(f(x) = \frac {1}{\sqrt {50 \pi }} \; \exp [-\frac {1}{50} (x-100)^2]\) para \(85 \leq x \leq 115\)

A ideia básica é criar um vetor com valores das abscissas (valores de \(x\)) e calcular o valor da função (valores de \(f(x)\)) para cada elemento da função e depois fazer o gráfico unindo os pares de pontos. Vejamos os comandos para o primeiro exemplo.

> x1 <- seq(0,50, l=101) 
> y1 <- 1 - (1/x1) * sin(x1) 
> plot(x1, y1, type="l", xlab = "x", ylab = "f(x)")

Note que este procedimento é o mesmo que aprendemos para fazer esboços de gráficos a mão em uma folha de papel, calculando e unindo ponto a ponto uma sequência de pontos!

PIC

Figura 1: Gráfico da função dada em (a).

Há ainda uma outra maneira de fazer isto no R utilizando plot.function() ou curve() conforme pode ser visto nos comandos a seguir que nada mais fazem do que combinar os três comandos anteriores em apenas um comando.

> plot(function(x) 1 - (1/x) * sin(x), from = 0, to = 50) 
> curve(1 - (1/x) * sin(x), from = 0, to = 50)

Vejamos agora como obter o gráfico para a segunda função.

> x2 <- seq(80, 120, l=101) 
> y2 <- (1/sqrt(50*pi)) * exp(-0.02 * (x2-100)^2) 
> plot(x2, y2, type="l", xlab = "x", ylab = "f(x)")

PIC

Figura 2: Gráfico da função dada em (b).

Note ainda que esta função é a densidade da distribuição normal e o gráfico também poderia ser obtido com:

> y2 <- dnorm(x2, 100, 5) 
> plot(x2, y2, type="l")

ou ainda:

> plot(function(x) dnorm(x, 100, 5), 85, 115)

6.3 Integração numérica

A função integrate() é usada para integração numérica em uma dimensão. Como exemplo vamos considerar resolver a seguinte integral: \begin {equation} I = \int _{-3}^3 x^2 d x. \label {eq:int} \end {equation} Para resolver a integral devemos criar uma função no R com a expressão da função que vamos integrar e esta deve ser passada para integrate() conforme este exemplo:

> fx <- function(x) x^2 
> integrate(fx, -3, 3)

18 with absolute error < 2e-13

A integral acima corresponde à área mostrada no gráfico da Figura 3 obtida com os comandos a seguir. Inicialmente é desenhada a função como nos exemplos anteriores. A área que corresponde à integração é desenhada por polygon() que utiliza uma sequencia de coordenadas \((x,y)\) que define o polígono. A função rev() reverte a ordem dos elementos de uma sequência.

> x <- seq(-4, 4, l=100) 
> x2 <- x^2 
> plot(x, x^2, ty='l') 
> x <- seq(-3, 3, l=100) 
> x2 <- x^2 
> polygon(rbind(cbind(rev(x),0),cbind(x,x2)), col='gray')

PIC

Figura 3: Gráfico onde a área indicada corresponde à integral definida na equação 2.

Vejamos mais um exemplo. Sabemos que para distribuições contínuas de probabilidades a integral está associada a probabilidade em um intervalo. Seja \(f(x)\) uma f.d.p. de uma variável contínua, então \(P(a < X < b) = \int _a^b f(x) dx\). Por exemplo, seja \(X\) uma v.a. com distribuição \(N(100, 81)\) e portanto \(f(x) = \frac {1}{9\sqrt {2\pi }}\exp \{-\frac {1}{162}(x-100)^2\}\). A probabilidade \(P(85 < X < 105)\) pode ser calculada das três formas diferentes que irão retornar os mesmos resultados conforme mostrado a seguir.

> fx <- function(x){(1/(9*sqrt(2*pi))) * exp(-(1/162)*(x-100)^2)} 
> integrate(fx, 85, 105)
0.6629523 with absolute error < 7.4e-15
> integrate(function(x) dnorm(x, 100, 9), 85, 105)
0.6629523 with absolute error < 7.4e-15
> pnorm(105, 100, 9) - pnorm(85, 100, 9)
[1] 0.6629523

6.4 Matemática simbólica no R

Embora o R seja um programa predominantemente para operações numéricas, é possível obter alguns resultados simbólicos, em particular para expressões de derivadas que podem ser informadas para algoritmos de otimização numérica. A forma básica de utilização consiste em:
(i) defina a expressão desejada dentro de quote(),
(ii) use D() para obter a expressão da derivada desejada informando a expressão e o termo em relação ao qual deseja-se derivar a expressão,
(iii) use eval() caso queira obter o valor numérico de uma determinada expressão. A documentação help(D) fornece mais detalhes.
Vejamos um exemplo.

> f <- quote(sin(x^2 + log(y+z))) 
> f

sin(x^2 + log(y + z))

> df.dx <- D(f, "x") 
> df.dx

cos(x^2 + log(y + z)) * (2 * x)

> df.dy <- D(f, "y") 
> df.dy

cos(x^2 + log(y + z)) * (1/(y + z))

> eval(f, list(x=1, y=2, z=3))

[1] 0.5073913

> eval(df.dx, list(x=1, y=2, z=3))

[1] -1.723432

Existem programas computacionais especializados em matemática simbólica dentre os quais destacam-se os projetos axiom, maxima, yacas e SymPy e ainda o site WolframAlfa.

Todos estes programas estão disponíveis para várias plataformas/sistemas operacionais. Também estão incluídos nas distribuições de sistemas LINUX podendo ser instalados a partir de versões pré-compiladas. Alguns ainda oferecem interfaces gráficas.

Para o YACAS há o pacote do R Ryacas ( cran.r-project.org/package=Ryacas) que permite acessar suas funcionalidades diretamente de dentro de uma sessão. O pacote rSymPy, que fazia o mesmo para o SymPy, foi removido do CRAN e não é mais mantido.

6.5 Exercícios

1.

Calcule o valor das expressões abaixo

(a)

Seja \(x = (12, 11, 14, 15, 10, 11, 14, 11)\).
Calcule \(E = -n\lambda + (\sum _1^n x_i) \log (\lambda ) - \sum _1^n \log (x_i !)\), onde \(n\) é o número de elementos do vetor \(x\) e \(\lambda = 10\).
Dica: o fatorial de um número pode ser obtido utilizando a função prod. Por exemplo o valor de \(5!\) é obtido com o comando prod(1:5).
Há ainda uma outra forma usando a função Gama e lembrando que para \(a\) inteiro, \(\Gamma (a+1) = a!\). Portanto podemos obter o valor de \(5!\) com o comando gamma(6).

(b)

\(E = (\pi )^2 + (2\pi )^2 + (3\pi )^2 + \ldots + (10\pi )^2\)

(c)

\(E = \log (x+1) + \log (\frac {x+2}{2}) + \log (\frac {x+3}{3}) + \ldots + \log (\frac {x+20}{20})\), para \(x=10\)

2.

Obtenha o gráfico das seguintes funções:

(a)

\(f(x) = x^{12} (1-x)^{8}\) para \(0 < x < 1\)

(b)

Para \( \phi =4 \),

\begin {eqnarray*} \rho (h) = \left \{ \begin {array}{ll} 1 - 1.5\frac {h}{\phi } + 0.5(\frac {h}{\phi })^3 \mbox { , se $h < \phi $} \cr 0 \mbox { , caso contrário} \end {array} \right . \end {eqnarray*}
3.

Considerando as funções acima calcule as integrais a seguir e indique a área correspondente nos gráficos das funções.

(a)

\(I_1 = \int _{0.2}^{0.6} f(x) d x\)

(b)

\(I_2 = \int _{1.5}^{3.5} \rho (h) d h\)

4.

Mostre os comandos para obter as seguintes sequências de números

(a)

1 11 21 31 41 51 61 71 81 91

(b)

1 1 2 2 2 2 2 3 3 3

(c)

1.5 2.0 2.5 3.0 3.5 1.5 2.0 2.5 3.0 3.5 1.5 2.0 2.5 3.0 3.5

5.

Escreva a sequência de comandos para obter um gráfico \(x\) versus \(y\), onde \(x\) é um vetor com 100 valores igualmente espaçados no intervalo \([-1,1]\) e \(y = \sin (x)*\exp (-x)\).

6.

Escreva uma sequência de comandos no R para calcular a soma dos 80 primeiros termos das séries:

(a)

\(1 + 1/3^2 + 1/5^2 + 1/7^2 + 1/9^2 + \ldots \)

(b)

\(1 - 1/2^2 + 1/3^2 - 1/4^2 + 1/5^2 - 1/6^2 + \ldots \)