Category

Data Science

Nested loop com o pacote purrr no R

By | Data Science

O pacote purrr é um dos pilares do tidyverse, ao lado do dplyr e do operador pipe em termos de utilidade e generalização. Um de seus principais objetivos é facilitar a intenção de fazer um loop com a sintaxe envolvida. Ao invés de ter que escrever um monte de códigos tediosos do tipo "for-loop" - que pouco diz ao leitor sobre o que o loop faz -, a função map do purrr pode ser lida quase como inglês puro e se encaixa perfeitamente em operações encadeadas (pipe).

Mas isso não é tudo que o purrr tem a oferecer. O pacote vai muito além das funcionalidades map, walk e suas variantes, ajudando ainda mais no trabalho com listas. Neste exercício vamos utilizar algumas dessas funcionalidades do purrr para aplicações do tipo nested loop.

Primeiro vamos entender o nested loop com um exemplo mínimo, sem usar o pacote purrr. Vamos supor que você tenha dois vetores numéricos de tamanhos diferentes e que você queira todas as combinações possíveis dos valores de ambos os vetores, concatenando os valores em um único vetor separado por um hífen. Isso poderia ser feito com o controle de fluxo for:


# Nested loop
for (i in 1:2) {
for (j in 1:4) {
print(paste(i, j, sep = "-"))
}
}

## [1] "1-1"
## [1] "1-2"
## [1] "1-3"
## [1] "1-4"
## [1] "2-1"
## [1] "2-2"
## [1] "2-3"
## [1] "2-4"

Com o purrr poderíamos fazer a mesma operação com a função walk neste contexto, onde o .y do loop de dentro é o vetor numérico 1:2 do loop de fora:


library(purrr)

# Nested loop
walk(1:2, ~walk(1:4, ~print(paste(.y, .x, sep = "-")), .y = .x))

## [1] "1-1"
## [1] "1-2"
## [1] "1-3"
## [1] "1-4"
## [1] "2-1"
## [1] "2-2"
## [1] "2-3"
## [1] "2-4"

Uma outra opção seria, em verdade, evitar o nested loop, já que é de leitura dificultosa e requer um trabalho adicional desnecessário. Portanto, seguindo o conselho do pai da linguagem de programação C++, Bjarne Stroustrup:

To become significantly more reliable, code must become more transparent. In particular, nested conditions and loops must be viewed with great suspicion. Complicated control flows confuse programmers. Messy code often hides bugs.

Poríamos criar antes uma combinação de valores usando a função cross2 e então iterar esses elementos em um loop:


library(magrittr)

cross2(1:2, 1:4) %>%
walk(
~print(paste(.[[1]], .[[2]], sep = "-"))
)

## [1] "1-1"
## [1] "2-1"
## [1] "1-2"
## [1] "2-2"
## [1] "1-3"
## [1] "2-3"
## [1] "1-4"
## [1] "2-4"

Muito melhor, concorda?

Agora vamos a um exemplo prático da vida real. Suponha que você esteja trabalhando com os dados desagregados da inflação brasileira, medida pelo IPCA do IBGE e que precise coletar esses dados usando sua API através do pacote sidrar no R. A coleta dos dados pela API tem uma limitação de 50 mil observações por requisição, portanto você precisa criar uma estratégia para coletar esses dados sem ultrapassar esse limite da API.

Uma maneira de fazer isso, não a melhor, é justamente através de um nested loop. Para operacionalizar vamos coletar os dados do IPCA provenientes das tabelas 1419 e 7060 e suas variáveis do SIDRA/IBGE:


library(sidrar)

# Vetor com códigos das tabelas do IPCA
tables <- c("Tabela 1419" = 1419, "Tabela 7060" = 7060)

# Vetor com códigos das variáveis da tabela
variables <- c(
"IPCA - Variacao mensal (%)" = 63,
"IPCA - Peso mensal" = 66,
"IPCA - Variação acumulada no ano (%)" = 69,
"IPCA - Variação acumulada em 12 meses (%)" = 2265
)

# Coletar dados com nested loop
ipca <- tables %>%
map_dfr(
~map_dfr(
.x = variables,
~get_sidra(
x = .y, # vetor de códigos das tabelas (tables)
variable = .x, # vetor de códigos das variáveis (variables)
period = "all"
),
.y = .x
)
)

A outra possibilidade é gerar primeiro uma lista com combinações entre os códigos de tabelas e códigos das variáveis e, então, iterar estes elementos em um único loop:


library(dplyr)

# Coletar dados iterando combinações no loop
ipca2 <- cross2(tables, variables) %>%
map_dfr(
~get_sidra(
x = .[[1]],
variable = .[[2]],
period = "all"
)
)

# Verificar se resultados são equivalentes
all.equal(
arrange(ipca, `Mês`),
arrange(ipca2, `Mês`)
)

## [1] TRUE

Vale enfatizar que eu apenas arranhei a superfície do assunto neste texto, mas espero ter pelo menos instigado você a fazer sua própria investigação sobre o que o purrr pode oferecer, além de tornar as iterações mais legíveis, principalmente para os casos de nested loop.

 

Hackeando o R: Utilizando o RJAGS

By | Hackeando o R

No Hackeando o R de hoje, vamos continuar a desenvolver o ferramental de estatística bayesiana, e para isso vamos mostrar como utilizar a interface do R para o programa JAGS (Just Another Gibbs Sampler), que é especializado em métodos de MCMC. No post da semana passada, apresentamos o algoritmo de Metropolis, que define um passo-a-passo para criar novas observações da distribuição desconhecida dos parâmetros que queremos estimar. Existem outros métodos de gerar essa amostragem, porém não iremos nos preocupar a fundo com eles, e deixar que o JAGS escolha o mais apropriado. Com isso, nosso trabalho passa a ser ter em mãos a distribuição a priori dos parâmetros, e a distribuição dos dados condicionada aos parâmetros. Para reproduzir os códigos desse post, você precisará baixar o JAGS, disponível aqui.

Conforme vimos, para fazer a estimação do parâmetro de interesse, precisamos de uma distribuição a priori, a distribuição condicional de uma amostra, e uma amostra que irá atualizar nossas crenças sobre o parâmetro. A estimação do MCMC será feita pelo JAGS, através do pacote rjags. Esse sistema trabalha com listas, logo devemos gerar uma lista da nossa amostra a partir dos dados dela. Os dados estão disponíveis aqui.

dados <- readr::read_csv('z15N50.csv')

dataList = list(y=dados$y)

Note que geramos uma lista com um único elemento, que é nossa amostra. O próximo passo é definir o modelo, na sintaxe do JAGS. Isso é feito com uma string, que é salva em um arquivo na mesma pasta em que estamos trabalhando. A sintaxe é de modo geral a mesma para os modelos, definindo para cada elemento da amostra (dentro do for) a sua distribuição em função do parâmetro, e a distribuição a priori do parâmetro. Abaixo, definimos uma priori Beta(1, 1), ou seja, uma distribuição uniforme.


modelString = "
model {
for ( i in 1:50 ) {
y[i] ~ dbern( theta )
}
theta ~ dbeta(1,1)
}
"
writeLines( modelString , con="TEMPmodel.txt" )

O último passo que resta para rodarmos o algoritmo é definirmos o ponto inicial. Isso é arbitrário, porém se escolhermos pontos iniciais ruins, a convergência pode demorar, de modo que é preferível que o ponto inicial esteja em uma região de alta probabilidade da posterior. Um método possível é selecionar diversos pontos e rodar múltiplos MCMCs, porém para reduzir o código, vamos utilizar a estimativa de MLE da amostra, que é a média. Novamente, precisamos colocar o elemento em uma lista, e o nome dele deve ser igual ao nome passado ao modelo acima.


thetaInit = sum(dados$y)/length(dados$y)
initsList = list(theta=thetaInit)

Com isso, basta carregar o pacote rjags e rodar a função jags.model. Antes de gerarmos nossa cadeia de Markov, tomamos duas medidas iniciais: a adaptação e o burn-in. A adaptação é uma série de amostras que o JAGS extrai para calibrar suas estimativas das condicionais, e também escolher a magnitude ótima para as tentativas de salto que são feitas a cada passo do algoritmo. Esse processo não gera uma cadeia de Markov, logo é descartado do resultado final.

Se nossa estimativa inicial é mal colocada, é possível que ela dê foco demasiado para regiões de baixa probabilidade da posterior nas primeiras observações, até que ela encontre regiões de maior probabilidade. Para evitar que isso contamine a distribuição estimada final, é padrão eliminarmos o trecho inicial do caminho percorrido, chamado de burn-in. Abaixo, deixamos 500 amostras para a adaptação, e 500 amostras de burn-in. Após isso, geramos 3334 da cadeia final com a função coda.samples(), resultando em cerca de 10000 observações, dado que estamos criando 3 cadeias.


library(rjags)

jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList , inits=initsList ,
n.chains=3 , n.adapt=500 )

update( jagsModel , n.iter=500 )

codaSamples = coda.samples( jagsModel , variable.names=c("theta") ,
n.iter=3334 )

Para visualizar o resultado final, vamos utilizar as visualizações do coda, pacote que gerou as amostras do MCMC. O gráfico da esquerda mostra o caminho de cada simulação, enquanto o da direita é uma combinação das distribuições estimadas de cada cadeia.


plot(codaSamples)

Criando uma versão "janela-móvel" de qualquer função

By | Data Science

O que são cálculos em janela móvel e por quê são importantes em análise de dados e de séries temporais? Neste texto exploramos a criação e aplicação de funções em janelas móveis usando o R.

Em análise de séries temporais, nada é estático. Uma correlação pode existir em um subconjunto temporal ou uma média pode variar de um dia para o outro. Os cálculos em janelas móveis simplesmente aplicam funções a um subconjunto de largura fixa (ou não) desses dados (também conhecido como uma janela móvel), indexando uma observação a cada cálculo.

Existem alguns motivos comuns pelos quais você pode querer usar um cálculo em janela móvel em análise de séries temporais:

  • Obter a tendência central ao longo do tempo (média, mediana)
  • Obter a volatilidade ao longo do tempo (desvio-padrão, variância)
  • Detectar mudanças de tendência (médias móveis rápidas vs. lentas)
  • Obter uma relação entre duas séries temporais ao longo do tempo (covariância, correlação)

O tipo de cálculo em janela móvel mais simples existente no R é o base::cumsum() que faz a soma acumulada de valores de um vetor, no qual cada janela móvel é definida como sendo todos os valores anteriores ao atual. Um exemplo visual de soma acumula de uma sequência de números de 1 a 10:

Entendida a importância e intuição do funcionamento dessas operações, passemos a alguns exemplos práticos. Demonstraremos como operacionalizar esses cálculos em janelas móveis com as principais opções de pacotes de R disponíveis atualmente para essa finalidade:

  • timetk
  • slider
  • runner

Essa lista com certeza não esgota todas as possibilidades, existem outras opções e fica ao critério do usuário escolher o que melhor se adequa ao seu contexto.

Pacotes

Para reproduzir os códigos aqui expostos é necessário alguns pacotes de R, o código a seguir verifica se os mesmos estão instalados, instala caso necessário e carrega os mesmos para a memória:

# Instalar/carregar pacotes
if(!require("pacman")) install.packages("pacman")
pacman::p_load(
"magrittr",
"dplyr",
"ggplot2",
"ggthemes",
"tidyquant",
"timetk",
"tidyr",
"slider",
"lubridate",
"runner"
)

O pacote {timetk}

O pacote timetk é um toolkit para visualização, tratamento e transformação de séries temporais e oferece as funções slidify() e slidify_vec() para operações em janela móvel.

Como exemplo inicial, utilizamos os dataset FANG que traz os preços diários de 4 ações de tecnologia listadas na NASDAQ:

# Dados de preços diários de ações (FANG)
df_fang <- tidyquant::FANG
df_fang
# # A tibble: 4,032 x 8
## symbol date open high low close volume adjusted
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FB 2013-01-02 27.4 28.2 27.4 28 69846400 28
## 2 FB 2013-01-03 27.9 28.5 27.6 27.8 63140600 27.8
## 3 FB 2013-01-04 28.0 28.9 27.8 28.8 72715400 28.8
## 4 FB 2013-01-07 28.7 29.8 28.6 29.4 83781800 29.4
## 5 FB 2013-01-08 29.5 29.6 28.9 29.1 45871300 29.1
## 6 FB 2013-01-09 29.7 30.6 29.5 30.6 104787700 30.6
## 7 FB 2013-01-10 30.6 31.5 30.3 31.3 95316400 31.3
## 8 FB 2013-01-11 31.3 32.0 31.1 31.7 89598000 31.7
## 9 FB 2013-01-14 32.1 32.2 30.6 31.0 98892800 31.0
## 10 FB 2013-01-15 30.6 31.7 29.9 30.1 173242600 30.1
## # ... with 4,022 more rows

Calculando uma média móvel de 30 períodos para o preço ajustado (adjusted):

df_fang %>%
dplyr::select(symbol, date, adjusted) %>%
dplyr::group_by(symbol) %>% 
# Criar média móvel
dplyr::mutate(
roll_avg_30 = timetk::slidify_vec(
.x = adjusted,
.f = mean, 
.period = 30, 
.align = "right", 
.partial = TRUE
)
) %>%
tidyr::pivot_longer(cols = c(adjusted, roll_avg_30)) %>%
# Gerar gráfico
timetk::plot_time_series(
date, 
value, 
.title = "Média móvel de 30 períodos",
.color_var = name,
.facet_ncol = 2, 
.smooth = FALSE, 
.interactive = FALSE
)

Para mais operações simultaneamente pode ser interessante primeiro definir as funções utilizando o slidify() para depois calcular a operação em janelas móveis:

# Média móvel
roll_avg_30 <- timetk::slidify(.f = mean, .period = 30, .align = "right", .partial = TRUE)
# Soma acumulada
roll_sum_30 <- timetk::slidify(.f = sum, .period = 30, .align = "right", .partial = TRUE)
# Desvio padrão móvel
roll_sd_30 <- timetk::slidify(.f = sd, .period = 30, .align = "right", .partial = TRUE)
# Regressão móvel (rolling regression)
roll_lm_90 <- timetk::slidify(~ lm(..1 ~ ..2 + ..3), .period = 90, .unlist = FALSE, .align = "right")

E então aplicar as funções criadas:

df_fang %>%
dplyr::select(symbol, date, adjusted, volume) %>%
dplyr::group_by(symbol) %>% 
# Operações em janelas móveis
dplyr::mutate(
numeric_date = as.numeric(date),
m_avg_30 = roll_avg_30(adjusted),
m_sum_30 = roll_sum_30(adjusted),
m_sd_30 = roll_sd_30(adjusted),
m_lm_90 = roll_lm_90(adjusted, volume, numeric_date)
) %>% 
dplyr::filter(!is.na(m_lm_90))

## # A tibble: 3,676 x 9
## # Groups: symbol [4]
## symbol date adjusted volume numeric_date m_avg_30 m_sum_30 m_sd_30
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FB 2013-05-10 26.7 30847100 15835 26.8 805. 0.842
## 2 FB 2013-05-13 26.8 29068800 15838 26.9 807. 0.805
## 3 FB 2013-05-14 27.1 24930300 15839 26.9 808. 0.756
## 4 FB 2013-05-15 26.6 30299800 15840 27.0 809. 0.747
## 5 FB 2013-05-16 26.1 35499100 15841 26.9 808. 0.762
## 6 FB 2013-05-17 26.2 29462700 15842 26.9 807. 0.766
## 7 FB 2013-05-20 25.8 42402900 15845 26.9 806. 0.794
## 8 FB 2013-05-21 25.7 26261300 15846 26.8 805. 0.822
## 9 FB 2013-05-22 25.2 45314500 15847 26.7 802. 0.863
## 10 FB 2013-05-23 25.1 37663100 15848 26.6 799. 0.880
## # ... with 3,666 more rows, and 1 more variable: m_lm_90 <list>

O pacote {slider}

O pacote slider oferece uma sintaxe familiar ao purrr para calcular operações em janelas móveis.

O exemplo mais simples de sua utilização é o cálculo de uma média móvel dado um vetor numérico. A janela móvel é indicada pelo argumento .before:

# Média móvel (2 "períodos") com alinhamento a direita
slider::slide_dbl(
.x = 1:5, 
.f = ~mean(.x), 
.before = 1, # valor atual + 1 valor anterior
.complete = TRUE 
)

O pacote oferece algumas funções derivadas de slide_dbl() para os cálculos mais comuns: slide_mean(), slide_sum(), slide_prod(), slide_min() e slide_max().

Para séries temporais o pacote oferece ainda um recurso interessante de realizar a operação com janela móvel relativa a um índice com as funções slide_index(). Se você já quis calcular algo como uma “média móvel de 3 meses”, em que o número de dias em cada mês é irregular, você pode gostar desta função.

No exemplo abaixo, primeiro é calculado uma média móvel de 90 períodos sem considerar o índice e depois é calculado outra média móvel, mas considerando um índice de 3 meses. A função slide_index_mean() realiza o cálculo da média móvel de 3 meses considerando o valor atual em relação a mesma data aproximada de 3 meses anteriores (não necessariamente 90 observações anteriores). O resultado é uma média móvel diária relativa a uma janela móvel de 3 meses:


# Calcular médias móveis
df_amzn_roll <- df_fang %>%
dplyr::filter(symbol == "AMZN") %>%
dplyr::select(date, symbol, adjusted) %>%
dplyr::mutate(
roll_avg_90 = slider::slide_mean(adjusted, before = 89, complete = TRUE),
roll_idx_3m = slider::slide_index_mean(
x = adjusted, # vetor numérico
i = date, # vetor índice (datas)
before = ~ .x %m-% base::months(3), # para evitar erro gerado por NAs
complete = TRUE
)
)

# Gerar gráfico
df_amzn_roll %>%
dplyr::filter(date <= lubridate::as_date("2014-06-30")) %>%
tidyr::pivot_longer(cols = c(adjusted, roll_avg_90, roll_idx_3m)) %>%
timetk::plot_time_series(
date,
value,
.title = "AMZN",
.color_var = name,
.smooth = FALSE,
.interactive = FALSE
)

Perceba que a média móvel que considera o índice de 3 meses tem início anterior, pois considera que os meses possuem número de dias irregulares e se ajusta a tal.

O pacote {runner}

O pacote runner oferece um framework arrojado para calcular operações com janelas móveis de tamanho fixo ou variante, opções de inclusão de lags, tratamento de NAs, computação paralela, e é desenhado para dados em formato de séries temporais ou longitudinais.

Exemplos simples de médias móveis com e sem lag:


# Média móvel de 2 "períodos"
runner::runner(1:5, f = mean, k = 2, na_pad = TRUE)

## [1] NaN 1.5 2.5 3.5 4.5

# Média móvel de 2 "períodos" com 1 lag
runner::runner(1:5, f = mean, k = 2, lag = 1, na_pad = TRUE)

## [1] NaN NaN 1.5 2.5 3.5

Às vezes, as observações de um conjunto de dados não são espaçados igualmente (ausência de fins de semana, feriados, etc.) e, portanto, o tamanho da janela móvel deve variar para manter o período de tempo esperado. Para essa situação, podemos especificar o argumento idx para que a função em execução seja aplicada nas janelas móveis dependendo da data. Neste exemplo abaixo a janela móvel possui o índice semanal variante conforme a data:


# Criar vetores de dados
dates <- Sys.Date() + cumsum(sample(1:5, 20, replace = TRUE)) # série temporal irregular
values <- 1:20

dates

## [1] "2021-08-23" "2021-08-24" "2021-08-26" "2021-08-30" "2021-09-04"
## [6] "2021-09-08" "2021-09-12" "2021-09-13" "2021-09-16" "2021-09-18"
## [11] "2021-09-23" "2021-09-26" "2021-09-27" "2021-09-30" "2021-10-03"
## [16] "2021-10-07" "2021-10-11" "2021-10-15" "2021-10-17" "2021-10-20"

values

## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

# Calcular média móvel com janela variante
runner::runner(
x = values,
k = "7 days",
idx = dates,
f = mean
)

## [1] 1.0 1.5 2.0 3.0 4.5 5.5 6.5 7.0 8.0 8.5 10.5 11.5 12.0 13.0 14.0
## [16] 15.5 16.5 17.5 18.0 19.0

Espero que esta exploração com exemplos de cálculos com funções em janelas móveis tenha sido proveitosa, apesar de breve e resumida. Recomendo fortemente consultar as documentações dos pacotes, pois os mesmos fornecem todas as informações, detalhes, motivação e mais exemplos.

 

Hackeando o R: MCMC - algoritmo de Metropolis

By | Hackeando o R

No Hackeando o R de hoje, vamos apresentar o método de MCMC, ferramenta de grande importância para a estimação dentro da estatística bayesiana. O exemplo que apresentamos na semana passada era de fácil computação (como pode ser visto pelos códigos utilizados), muito por causa do fato de que a equação encontrada possuía solução fechada. Quando possuímos poucos parâmetros para estimar, e o espaço de possibilidades dos parâmetros é pequeno, também podemos resolver 'manualmente' calculando as probabilidades, pois é questão de percorrer cada uma das possibilidades e encontrar o valor do Teorema de Bayes. Para problemas mais complexos, como os que encontramos na vida real, essas soluções fáceis muitas vezes não estão disponíveis. Pense no caso de uma probabilidade condicional sem fórmula fechada, para a estimação de 6 parâmetros, que podem assumir 1000 valores cada um. Estaríamos olhando para o cálculo de 1000^6 combinações, o que é pesado até para computadores modernos.

A solução para esse problema é o chamado MCMC (Monte-Carlo Markov Chain). Apesar dele ser mais utilizado em problemas com múltiplos parâmetros, vamos mostrar como ele funciona com apenas um, para facilitar a compreensão. O ponto importante do método é que ele só depende da estimação dos termos no numerador do Teorema de Bayes (mostrado abaixo), eliminando qualquer preocupação com a integral do denominador, que é o grande problema para encontrarmos soluções fechadas.

A ideia básica do MCMC é gerar uma estimativa da distribuição a posteriori a partir de amostras dela, sem que nós precisemos construir ela per se. Hoje, vamos mostrar o algoritmo de Metropolis, que gera uma aproximação da distribuição a posteriori através de regras de decisão simples. O algoritmo utiliza apenas a razão entre duas probabilidades a posteriori, logo há estimação da integral acaba sendo desnecessária. Dada uma posição inicial para os parâmetros, os passos do algoritmo são simples:

1: geramos um movimento aleatório dos parâmetros no seu espaço, que será testado;
2: Verificamos se a posteriori na nova posição tem valor maior que na original. Se sim, o movimento ocorre. Caso contrário, o movimento ocorrerá com probabilidade igual à razão entre as posterioris;
3: geramos uma observação uniforme de 0 a 1, e comparamos ela com a razão calculada, de modo a validar o movimento acima.

Com isso, vamos testar o algoritmo para uma distribuição de eventos Bernoulli. Sabemos que essa distribuição possui fórmula fechada, logo iremos utilizá-la para comparar com a distribuição 'aproximada' encontrada de modo empírico.

dados = c(rep(0,6),rep(1,14))

probParcial = function(theta, data){

z = sum(data)
N = length(data)
p_Xi_dado_Theta = theta^z * (1-theta)^(N-z)
p_Xi_dado_Theta[theta > 1 | theta < 0] = 0

pTheta = dbeta(theta , 1, 1)
pTheta[theta > 1 | theta < 0] = 0

parcial = p_Xi_dado_Theta * pTheta
return(parcial)
}

n = 50000
trajetoria = rep(0, n)
trajetoria[1] = 0.5

set.seed(12334)

for ( t in 1:(n-1) ) {
posicao = trajetoria[t]

choque = rnorm(1, mean = 0, sd = 0.2)

probabilidade = min(1,
probParcial(posicao + choque, dados)
/ probParcial(posicao, dados))

if ( runif(1) < probabilidade ) {

trajetoria[t+1] = posicao + choque

} else {

trajetoria[t+1] = posicao

}
}

df = data.frame(passo = 1:50000, valor = trajetoria)

library(ggplot2)

ggplot(df, aes(x=valor)) + geom_histogram(aes(y = stat(density))) +
stat_function(fun = function(x) dbeta(x, 15, 7), color = "red",
size = 1)

Como paralelizar código no R

By | Data Science

Computação paralela é um tema complexo e vasto, mas ainda pouco conhecido. Graças ao esforço da comunidade de R já foram criados diversos pacotes que tornam a paralelização de códigos muito mais fácil e segura. Neste texto abordamos alguns exemplos práticos e introdutórios para demonstrar os ganhos dessa abordagem.

O objetivo é apresentar instrumentos e demonstrar o ganho imediato de paralelizar o seu código de R. Para detalhes mais técnicos e uma discussão mais aprofundada do tema recomendamos o capítulo 4 do livro Data Science for Economists and Other Animals de McDermott & Rubin.

Os pacotes utilizados podem ser instalados/carregados conforme o código abaixo:

if(!require("pacman")) install.packages("pacman")
pacman::p_load(
"tictoc",
"purrr",
"future",
"parallelly",
"furrr"
)

Exemplo: raiz quadrada super lenta

Nosso exemplo prático vai ser o cálculo de uma raiz quadrada realizado de forma super lenta. Para isso criaremos a função abaixo que calcula a raiz quadrada de um número mas em uma versão "lenta", ou seja, obriga o computador a esperar 2 segundos antes de fazer qualquer coisa. Isso servirá apenas para fins didáticos de demonstração da paralelização que faremos à seguir:

raiz_lenta <- function(x) {
Sys.sleep(2)
x_sqrt <- sqrt(x)
return(x_sqrt)
}

Execução serial

Agora vamos usar essa função criada iterando uma sequência de números de forma serial, ou seja, os cálculos não serão executados em paralelo. Utilizaremos aqui abordagem do pacote purrr do tidyverse para iterações no R. A função map_dbl() retorna um vetor numérico com os resultados dos cálculos. Para contabilizar o tempo de execução do código utilizaremos o pacote tictoc:


tictoc::tic()
exemplo_serial <- purrr::map_dbl(1:10, raiz_lenta)
tictoc::toc(log = TRUE)

## 20.25 sec elapsed

Conforme esperado, a iteração durou 20 segundos para ser concluída, pois forçamos uma parada de 2 segundos na função raiz_lenta() a cada iteração serial. Isso significa que podemos facilmente acelerar esse código se rodarmos a iteração em paralelo.

Execução paralela

Para rodar o código em paralelo primeiro é importante saber que o procedimento depende do número de núcleos de processador disponíveis no computador. Uma maneira fácil de obter essa informação no R é através do código abaixo:


parallelly::availableCores()

## system
## 4

Em segundo lugar, é preciso configurar o R para executar o código em paralelo, e é neste ponto que o pacote future, criado pelo Henrik Bengtsson, entra em cena.

Neste exemplo sugerimos a utilização do pacote future pois o mesmo oferece uma sintaxe simples e intuitiva para quem já está acostumado com o tidyverse. Graças a este pacote, o único procedimento que devemos fazer é especificar uma estratégia, ou seja, executar o código de maneira serial, paralela, etc., através da função plan(). Como queremos executar em paralelo, especificamos a estratégia multisession() que cria várias sessões do R no computador para executar o código.

Geralmente recomenda-se alocar para utilização em paralelo a quantidade total de núcleos de processador deixando um livre para outras tarefas do computador - como, por exemplo, navegação na internet, editor de texto, etc. -, mas aqui utilizaremos 2 núcleos apenas. Para outras opções consulte a documentação da função multisession().


future::plan(
strategy = future::multisession(workers = 2)
)

Realizada a configuração do R para execução em paralelo, vamos voltar ao nosso exemplo de iteração do cálculo de raiz quadrada. Executaremos novamente essa iteração, mas dessa vez em paralelo. Para isso faremos uma modificação marginal no código que é utilizar uma variante da função purrr::map_dbl() específica para iteração em paralelo, que é a função furrr::future_map_dbl().


tictoc::tic()
exemplo_paralelo <- furrr::future_map_dbl(1:10, raiz_lenta)
tictoc::toc(log = TRUE)

## 10.58 sec elapsed

Um resultado empolgante: o tempo de execução caiu pela metade! Isso tudo com pouquíssimas alterações em nosso código, bastou configurar o R (através do future::plan) e alterar a função de iteração (furrr::future_map_dbl). Um ganho de performance fenomenal, não?!

Vamos verificar se o output é idêntico nos dois casos (serial e paralelo):


all.equal(exemplo_serial, exemplo_paralelo)

## [1] TRUE

Agora você já pode iniciar um olhar mais crítico sobre os seus códigos para identificar pontos de melhoria de performance.

 

Receba diretamente em seu e-mail gratuitamente nossas promoções especiais
e conteúdos exclusivos sobre Análise de Dados!

Assinar Gratuitamente
{"cart_token":"","hash":"","cart_data":""}