Introdução à Inferência Variacional.

A Inferência Variacional é uma abordagem para se estimar distribuições a posteriori muito difíceis ou complexas.

Frequentemente, em Inferência Bayesiana estamos interessados em computar as medidas $p(z|x)$ para variáveis latentes $z_{1},\dots,z_{n}$ dado as observações $x_{1},\dots,x_{n}$. Quando essa distribuição é intratável podemos aproximar $p(z|x)$ por $q(z)$ onde a medida de proximidade é dada pela divergência de Kullback-Leibler.

Usualmente a Inferência Variacional é realizada através dos seguintes passos:

1. Comece com um modelo para $p(z,x)$.
2. Encontre a aproximação variacional $q(z|\nu)$
onde $\nu$ são os parâmetros variacionais.
3. Encontre a Evidence Lower Bound (ELBO) dada por $L(\nu)=E_{Q}[\log(p(z,x))-\log(q(z|\nu))]$
4. Minimize a quantia $L(\nu)$
5. Volte para o passo 3 e faça isso até haver convergência.

Com base no exposto, vamos incialmente considerar um modelo de Regressão Probit construído da seguinte forma:

?

1 2 3 4 5 6 7 8

#Define a semente set.seed(8333) #Gera os dados x<-rnorm(1000) #Gera a latente z<- rnorm(1000,10+5*x,1) #Gera a variável dependente: y<-ifelse(z<0,0,1)

O objetivo é estmar os parâmetros $a=10,b=5$ do modelo de Regressão Probit. Note que o Processo Gerador de Dados temos que a distribuição a posteriori é dada por: $p(a,b,z|y,x)= \prod_{i=1}^{n}y_{i}^{I(z_{i}>0)}(1-y_{i})^{1-I(z_{i}>0)}\times \frac{1}{\sqrt{2\pi}}\exp\left[-\frac{(z_{i}-(a+bx_{i}))^{2}}{2}\right]$ o qual também pode ser escrito como: $p(a,b,z|y,x)\propto \prod_{i=1}^{n}y_{i}^{I(z_{i}>0)}(1-y_{i})^{I(z_{i}\leq 0)}\times \exp\left[-\frac{(z_{i}-(a+bx_{i}))^{2}}{2}\right]$ para distribuições a priori uniformes para $a,b$. Dado o modelo precisamos definir a aproximação variacional. Nesse caso Blei et. al. (2017) sugere como regra ótima para a aproximação variacional fazer: $q(\theta_{j})\propto \exp\left[E_{\theta_{-j}}(\log(p(\theta|y,x)))\right]$ onde $\theta=(\theta_{1},\dots,\theta_{p})$ são os parâmetros e variáveis latentes do modelo. Nesse sentido, temos: $\log(p(a,b,z|y,x))=\sum_{i=1}^{n}I(z_{i}>0)\log(y_{i})+I(z_{i}\leq 0)\log(1-y_{i})-\sum_{i=1}^{n}\frac{(z_{i}-(a+bx_{i}))^{2}}{2}$ A distribuição varicional ótima é dada por: $q_{j}(z_{j})\propto \exp\left[E_{a,b,z_{-j}}(\log(p(a,b,z|y,x)))\right]$ Note que: $E_{a,b,z_{-j}}(\log(p(a,b,z|y,x)))=I(z_{j}>0)\log(y_{j})+I(z_{j}\leq 0)\log(1-y_{j})-\frac{E_{a,b}(z_{j}-(a+bx_{j}))^{2}}{2}$ Aplicando a exponencial temos que a distribuição ótima para $z_{j}$ é uma distribuição Normal Truncada, isso é $q_{j}(z_{j})\sim N^{+}(E(a)+E(b)x_{j},1)$ se $y_{j}=1$ e $q_{j}(z_{j})\sim N^{-}(E(a)+E(b)x_{j},1)$ se $y_{j}=0$. Similarmente, para $a$ temos $E_{b,z}[\log(p(a,b,z|y,x))]\propto E_{b,z}\left(-\frac{\sum_{i=1}^{n}(z_{i}-(a+bx_{i}))^{2}}{2}\right)$ , removendo os termos que não dependem de $a$ e complentando quadrados temos: $q_{a}(a)\sim N(\frac{\sum_{i=1}^{n}[E(z_{i})-E(b)x_{i}]}{n},\frac{1}{n})$. Finalmente para $b$ temos $E_{a,z}[\log(p(a,b,z|y,x))]\propto E_{a,z}\left(-\frac{\sum_{i=1}^{n}(z_{i}-(a+bx_{i}))^{2}}{2}\right)$, novamente removendo os termos que não dependem de $b$ e complentando quadrados temos: $q_{b}(b)\sim N(\frac{\sum_{i=1}^{n}[E(z_{i})-E(a)]x_{i}}{\sum_{i=1}^{n}x_{i}^{2}},\frac{1}{\sum_{i=1}^{n}x_{i}^{2}})$. Nesse caso, a otimização ocorre atribuindo os valores dos parâmetros a sua respectiva média da distribuição variacional. As equações para atualização dessas distribuições são dadas por:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

#Equações de Update update.newZj = function(newA,newB,j) {   mu = newA + newB*x[j]   if (y[j] == 1) {     return(mu + dnorm(-1*mu)/(1-pnorm(-1*mu)))   } else {     return(mu - dnorm(-1*mu)/(pnorm(-1*mu)))   } } update.newA = function(newZ,newB) {   return(sum(newZ-newB*x)/n) } update.newB = function(newZ,newA) {   return(sum(x*(newZ-newA))/sum(x^2)) }

Finalmente, podemos construir um sistema iterativo baseado nos valores otimizados:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

#Total de observações n<-length(z) #Chute inicial newA<-3 newB<-9 #Inicializa vetores a.vec<-rep(NA,n) b.vec<-rep(NA,n) newZ<-rep(NA,n) #Inicia o processo iterativo for(it in 1:1000){   #Update Z   for (i in 1:n) {     newZ[i] = update.newZj(newA,newB,i)   }   #Correct Numerical   newZ[is.infinite(newZ)]<-100   #Update A   newA = update.newA(newZ,newB)   newB = update.newB(newZ,newA)   a.vec[it] = newA   b.vec[it] = newB }     plot(a.vec,type="l") summary(a.vec) plot(b.vec,type="l") summary(b.vec)

Usando a função DEoptim em paralelo.

Um dos pacotes mais interessantes para realizar otimização é o pacote DEoptim. Essa abordagem permite a otimização de funções dos mais diversos tipos de uma maneira fácil.

Entretanto, por usar Otimização Evolucionária o processo de otimização pode ser lento, caso utilize-se muitas populações, iterações ou níveis tolerância extremos.

Aqui mostrarei como usar a opção parallelType da função DEoptim para que esse processo seja acelerado por meio da programação em paralelo.

Considere a seguinte função com múltiplos máximos e mínimos locais:

?

1 2 3 4

f<-function(x,y){   res<-(268/((3*x+18)^2+(3*y+7)^2+(18+7)))-(2*(268)/((3*x-21)^2+(3*y+21)^2+(21+21)))-(3*(268)/((3*x+18)^2+(3*y-21)^2+(18+21)))+(4*(268)/((3*x-21)^2+(3*y-7)^2+(21+7)))   return(res) }

Podemos fazer o gráfico dessa função através do seguinte bloco de comandos:

?

1 2 3 4 5 6 7 8

#Chama a biblioteca library(plot3D) #Cria a grade de valores na função z<-f(M[,"x"], M[,"y"]) #Faz o gráfico persp3D(z = z, x = M[,"x"], y = M[,"y"],         expand = 0.3, main = "Função complicada", facets = FALSE, scale = FALSE,         clab = "Valor da função",zlab="Valor", colkey = list(side = 1, length = 0.5))

O próximo passo para realizar a otimização é escrever a função no formato adequado para o uso do método DEoptim. Para isso, basta trocar o argumento da função que ao invés de serem dois parâmetros, será agora um vetor contendo esses dois argumentos, em outras palavras:

?

1 2 3 4 5 6

fEval<-function(parms){   x<-parms[1]   y<-parms[2]   res<-(268/((3*x+18)^2+(3*y+7)^2+(18+7)))-(2*(268)/((3*x-21)^2+(3*y+21)^2+(21+21)))-(3*(268)/((3*x+18)^2+(3*y-21)^2+(18+21)))+(4*(268)/((3*x-21)^2+(3*y-7)^2+(21+7)))   return(res) }

Note que só as três primeiras linhas foram alteradas. A seguir, precisamos configurar os nossos clusters para que a função possa ser otimizada em paralelo, a função busca pelo menor valor da função:

?

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

#Chama o pacote DEoptim e doSNOW library(DEoptim) library(doSNOW) #Guarda o número de clusters disponíveis numSlaves<- detectCores() #Cria os clusters cl <− makeCluster(numSlaves) #Passa todas as bibliotecas que são usadas para o cálculo da função #clusterEvalQ(cl , library(pracma, ... )) #Passa todos os objetos necessários para o cômputo da função #clusterExport(cl , ... ) #Registra os clusters registerDoSNOW (cl) #Realiza a otimização outDEoptim <− DEoptim (fEval, c(-15,-15), c(15,15), DEoptim.control(list(trace=FALSE, NP=50, itermax=50, parallelType = 2))) #Fecha os clusters stopCluster(cl) #Guarda os melhores resultados. fim <- outDEoptim $ optim $ bestmem fim

Alguns comentários são importantes: as linhas 8 e 10 são utilizadas para se passar bibliotecas e objetos do R que são utilizados no cômputo da função que deve ser otimizada. Nesse caso, como não usamos nenhum pacote especial para o cômputo da função e nem bibliotecas esses comandos estão comentados. Caso deseje-se passar esses valores, basta atribuí-los na forma de lista em cada comando. Na função DEoptim o primeiro argumento é a função a ser minimizada, seguido pelos limites inferiores e superiores para os argumentos. O último argumento é uma lista que controla as condições da função de otimização: NP é o número de populações a serem geradas no algoritmo de Otimização Evolucionária, itermax é o número total de iterações e parallelType = 2 define a forma de paralelismo a ser considerada. Observe ainda que se o valor de NP e/ou itermax e o paralelismo não for adotado, esse processo pode ser muito lento computacionalmente, daí a vantagem de se utilizar o processamento em paralelo que pode reduzir esse tempo computacional significativamente. Mote ainda que o comando trace=FALSE é para suprimir os resultados em cada iteração, e caso deseje-se maximizar a função ao invés de minimizar, basta trabalhar com o valor do negativo da função, ou ainda com seu recíproco.

Orientação a objetos no R - S4 Classes e Métodos.

Uma das principais atividades e, também uma boa prática na programação usando o R é a utilização da abordagem de programação orientada a objetos. Essa proposta é útil na construção de pacotes no R bem como na organização e estabilidade das bibliotecas criadas. Nesse tutorial iremos trabalhar com a classe S4 de objetos do R.

Suponha que desejamos criar um pacote que nos ajude a calcular a menção de alunos de uma turma. Nesse caso, o sistema de classes S4 pode ser definido criando-se o objeto aluno:

?

1 2 3 4 5

setClass(   "Aluno",   representation(nome="character", idade="numeric"),   prototype(nome=character(0), idade=numeric(0)) )

Esse objeto contêm, inicialmente, dois slots: a idade do aluno e seu nome. O nome é considerado nesse modelo uma variável caracter e a idade uma variável numérica. O comando prototype define quais devem ser os valores iniciais para esses dois atributos.

Entretanto, na universidade nós temos diversos tipos de aluno, eles podem ser: alunos da graduação, pós-graduação, extensão, aluno especial, ouvinte, etc. Suponha que tenhamos dois tipos de alunos apenas:

?

1 2 3 4 5 6 7 8

#Define as subclasses de alunos setClass("Graduacao",   representation(nota1="numeric",nota2="numeric",nota3="numeric"),   contains="Aluno")   setClass("Posgraduacao",   representation(artigo="character",nota1="numeric",nota2="numeric"),   contains="Aluno")

Na primeira subclasse donominada Graduação o aluno faz três provas representadas pelos slots: nota1, nota2 e nota3. Já na segunda subclasse o aluno é um aluno de Pós-graduação e ele é avaliado de maneira diferente, tem que preparar um artigo o qual pode ser armazenado no slot artigo e faz duas provas cujas notas são armazenadas em nota1 e nota 2. Note que nesse caso, ambos são Alunos e por isso herdam os atributos de nome e idade.

Além do mais, na criação dessas subclasses não atribuímos os valores iniciais desses atributos por opção. A boa prática, no entanto, recomenda que todos esses atributos sejam definidos os valores de inicialização através do método prototype.

Cada tipo de aluno terá sua nota calculada de maneira diferente: para os alunos de graduação a nota final será a média das três notas nas provas, já os alunos de Pós-graduação a nota será 50% da média das duas provas e 50% se ele fez o artigo e zero se não fez. Nesse caso, usaremos a propriedade de Mutabilidade (Mutability) dos sistemas de programação orientada a objetos. Através dessa propriedade um cômputo diferente será realizado para cada tipo de objeto. O primeiro passo é a criação de um método genérico:

?

1 2 3 4 5 6 7 8

#Cria um método generico denominado calculaMencao #o qual é aplicado para o tipo de objeto de entrada. setGeneric(   "calculaMencao",   function(object) {     standardGeneric("calculaMencao")   } )

Em seguida, vamos definir um cálculo de função para cada tipo de aluno:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

#Para alunos de graduação setMethod(   "calculaMencao",   signature("Graduacao"),   function(object) {     #Calcula a média     media<-(object@nota1+object@nota2+object@nota3)/3     #Retorna o resultado     return(media)   } )   #Para alunos de pós-graduação setMethod(   "calculaMencao",   signature("Posgraduacao"),   function(object) {     #Calcula a média     media<-(object@nota1+object@nota2)/2     #Fez o artigo ?     artigo<-ifelse(object@artigo == "",0,10)     #Calcula a nota final     nota<-media*0.5+artigo*0.5     #Retorna o resultado     return(nota)   } )

Note que em ambas as definições dos métodos incluímos o método genérico calculaMencao. O que muda entre eles é a assinatura (signature) e, consequentemente, a maneira como a função é calculada. Vamos agora criar alguns alunos seguindo a classes S4 formada:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

#Aluno de Graduação joao <- new("Graduacao",               nome="João Silva",               idade=22,               nota1=2.75,               nota2=8.78,               nota3=5.72)   #Aluno de Pós-Graduação maria <- new("Posgraduacao",             nome="Maria Silva",             idade=27,             nota1=8.54,             nota2=3.21,             artigo="Como programar em R")   #Aluno de Pós-Graduação alex <- new("Posgraduacao",              nome="Alexandre Silva",              idade=33,              nota1=2.15,              nota2=5.48,              artigo="")

Para não ter que repetir o mesmo código toda vez que desejarmos criar um objeto, podemos construir os constructors:

?

1 2 3 4 5 6 7 8 9

#Cria o constructor: Graduacao <- function(nome, idade, nota1, nota2, nota3){   new("Graduacao", nome=nome, idade=idade, nota1=nota1, nota2=nota2, nota3=nota3) }   #Cria o constructor: Posgraduacao <- function(nome, idade, nota1, nota2, artigo){   new("Graduacao", nome=nome, idade=idade, nota1=nota1, nota2=nota2, artigo=artigo) }

Podemos calcular a menção deles fazendo:

?

1 2 3 4 5 6 7 8 9 10 11

#Calcula a nota do João resultadoJoao<-calculaMencao(joao) resultadoJoao   #Calcula a nota da Maria resultadoMaria<-calculaMencao(maria) resultadoMaria   #Calcula a nota do Alex resultadoAlex<-calculaMencao(alex) resultadoAlex

Podemos ainda ter objetos que são oriundos de mais de uma classe, por exemplo, considere um aluno que faz tanto Graduação quanto Pós-graduação, podemos criar uma classe específica para ele:

?

1 2	setClass("Maluco",          contains=c("Graduacao","Posgraduacao"))

A criação dos objetos e cálculo das notas é feito da seguinte forma:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

#Registro do Sergio na graduação sergioGraduacao <- new("Maluco",                    nome="Sergio Silva",                    idade=22,                    nota1=8.36,                    nota2=9.12,                    nota3=5.64)   #Registro do Sergio na pós-graduação sergioPosgraduacao <- new("Posgraduacao",                       nome="Sergio Silva",                       idade=22,                       nota1=3.45,                       nota2=8.34,                       artigo="Não sei nada")     #Calcula da nota do Sergio resultadoSergio1<-calculaMencao(sergioGraduacao) resultadoSergio1   resultadoSergio2<-calculaMencao(sergioPosgraduacao) resultadoSergio2

É comum ainda criar para os métodos alguns pontos de validação dos valores de entrada, por exemplo, nenhuma nota pode ser negativa ou maior do que 10. Podemos incluir essas restrições fazendo:

?

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

setClass("Graduacao",          representation(            nota1="numeric",            nota2="numeric",            nota3="numeric"),          validity = function(object) {            if (object@nota1<0 | object@nota1>10 )            {              stop("A nota na prova 1 tem que ser válida.")            }            if (object@nota2<0 | object@nota2>10 )            {              stop("A nota na prova 2 tem que ser válida.")            }            if (object@nota3<0 | object@nota3>10)            {              stop("A nota na prova 3 tem que ser válida.")            }            return(TRUE)          },          contains="Aluno"          )

Nas classes S4 de objetos temos algumas funções especiais:

• slotNames: fornece o nome dos slots.
• getSlots: fornece o nome dos slots e seus tipos.
• getClass: fornece o nome dos slots, seus tipos e das classes antecessoras.

CAViaR - Conditional Autoregressive Value at Risk by Regression Quantiles

O modelo CAViaR - Conditional Autoregressive Value at Risk by Regression Quantiles é um modelo utilizado para se mensurar o Value at Risk para períodos diários e, assim, monitorar o risco de maneira autogregressiva.

Engle e Manganelli (2004) sugerem quatro possíveis funções para a estimação do CAViaR:

• Adaptive: $VaR_{t}=VaR_{t-1}+\beta_{1}\left\{\left[1+\exp\left(G\left[y_{t-1}-VaR_{t-1}\right]\right)^{-1}\right]-\tau\right\}$.
• Symmetric absolute value: $VaR_{t}=\beta_{1}+\beta_{2}VaR_{t-1}+\beta_{3}|y_{t-1}|$.
• Asymmetric slope: $VaR_{t}=\beta_{1}+\beta_{2}VaR_{t-1}+\beta_{3}(y_{t-1})^{+}+\beta_{4}(y_{t-1})^{-}$.
• Indirect GARCH(1, 1): $VaR_{t}=\left(\beta_{1}+\beta_{2}VaR_{t-1}^{2}+\beta_{3}y_{t-1}^{2}\right)^{1/2}$

onde $y_{t}$ é o retorno do ativo no tempo $t$ com $t=1,\dots, T$. $G$ é uma constante positiva (Engle e Manganelli (2004) sugerem $G=10$), $\tau$ é o nível do VaR desejado (usualmente $\tau=0.01$ ou $\tau=0.05$) e $(x)^{+}=\max(x,0)$ e $(x)^{-}=-\min(x,0)$

Como as funções são recursvivas é necessário uma primeira estimativa para o $VaR_{t=1}$. Considerando os mesmos dados de Engle e Manganelli (2004) os quais trabalharam com os retornos dos preços das açoes da GM, IBM e SP500 para Abril 7, 1986 até Abril 7, 1999, podemos obter a primeira estimativa de $VaR_{t=1}$ como:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14

#Limpa os workspace rm(list=ls()) #Lê os dados stocks<-read.table("dataCAViaR.txt") colnames(stocks)<-c("GM","IBM","SP500") #Nível do VaR desejado tau<-0.05 #Divide a base em treinamento e validação train<-stocks[1:2892,] valid<-stocks[2893:3392,] #Criar o vetor de VaR VaR<-rep(NA,nrow(train)) #Primeira estimativa VaR[1]<- -qnorm(tau)*sd(train$GM)

Em seguida, Engle e Manganelli (2004) sugerem estimar os parâmetros dos modelos por meio da estrutura de Regressão Quantílica: $\min_{\boldsymbol\beta}\frac{1}{T}\sum_{t=1}^{T}\left[\tau-1\left(y_{t} \leq VaR_{t}\right)\right]\left[y_{t}-VaR_{t}\right]$ onde $1(\cdot)$ é uma função indicadora que assume valor igual a 1 quando seu argumento é verdadeiro $y_{t}\lt VaR_{t}$ e 0 caso contrário. Assim, a função objetivo para o modelo Adaptive é dada por:

?

1 2 3 4 5 6 7 8 9 10 11 12

#Sugestão Engle e Manganelli (2004) G<-10 #Function adpatative<-function(x){   beta1<-x[1]   for(i in 2:nrow(train)){     VaR[i]<- VaR[i-1] + beta1*((1+(exp(G*(train[i-1,"GM"] - VaR[i-1]))^(-1)))-tau)   }   #Objective Function   res<-sum((tau-(train$GM < VaR))*(train[,"GM"]-VaR))/nrow(train)   return(res) }

O próximo passo é otimizar essa função, isso pode ser feito utilizando-se a função optim:

?

1 2 3 4

#Realiza a otimização usando o métood de Brent res<-optim(par=c(0), adpatative, method="Brent", lower=c(-10), upper=c(10)) #Resultado do parâmetro de interesse beta<-res$par

A função optim tem como argumentos nesse exemplo: valor inicial para o(s) parâmetro(s) de interesse, função a ser minimizada, método utilizado e limites inferiores e superiores para o(s) parâmetro(s) de interesse. O último passo é a previsão do VaR:

?

1 2 3 4 5 6 7 8 9 10 11

VaR<-rep(NA,nrow(dados)) VaR[1]<- -qnorm(tau)*sd(dados[,"GM"]) #Previsão do CAViaR adpatativeForecast<-function(x){   beta1<-x[1]   for(i in 2:nrow(dados)){     VaR[i]<- VaR[i-1] + beta1*((1+(exp(G*(dados[i-1,"GM"]-VaR[i-1]))^(-1)))-tau)   }   #Objective Function   return(VaR) }

onde o objeto dados pode ser as bases de treinamento e validação:

?

1 2 3 4 5 6

#Base de treinamento dados<-train #Faz a previsão VaR.train<-adpatativeForecast(beta) #Medida de acuracia hits mean(dados$GM< -VaR.train)-tau

O mesmo processo pode ser realizado para as demais funções.

RcppEigen: Operações Matriciais e Rcpp - Parte 3

A biblioteca RcppEigen é uma das mais eficientes para o cômputo numérico de operações algébricas uma vez que utiliza algortimos optimizados para a memória cache e SIMD.

Por exemplo, considere o desafio de estimar uma Matriz de Variâncias e Covariâncias com base em um conjunto de dados. Essa estimação levaria em conta pelo menos $N\times(N-1)/2$ operações (calculando todas as variâncias e covariâncias na matriz triangular uma vez que a matriz é simétrica). Esse procedimento levaria em consideração o uso de loops aninhados (nested loops) os quais podem ser muito ineficientes e demorados quando implementados diretamente em R.

Nesse sentido, podemos usar as operações do pacote RcppEigen também documentadas na página do Eigen C++:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

// [[Rcpp::export]] #include <rcppeigen.h> // [[Rcpp::depends(RcppEigen)]] #include <cmath> using namespace Rcpp;   // [[Rcpp::export]] Eigen::MatrixXd covCalcula(Eigen::MatrixXd mat) {   //Centra as variáveis com respeito a média   Eigen::MatrixXd centered = mat.rowwise() - mat.colwise().mean();   //Calcula a covariância   Eigen::MatrixXd cov = (centered.adjoint() * centered) / double(mat.rows() - 1);   return cov; } </cmath></rcppeigen.h>

Observe que as operações .rowwise() e .colwise() se referem, respectivamente, aos cômputos para cada linha e para cada coluna. Já o método .adjoint() calcula a Matriz Adjunta.

A invocação desse método no R é realizada da seguinte forma:

?

1 2 3 4 5 6 7 8 9

#Limpa o Working Directory #Chama um conjunto de dados data("ChickWeight") #Converte todas as variáveis para numeric dat2 <- data.frame(lapply(ChickWeight, function(x) as.numeric(as.character(x)))) #Transforma para matriz mat<-as.matrix(dat2) #Calcula a matriz de covariância covCalcula(mat)

RcppEigen: Operações Matriciais e Rcpp - Parte 2

Uma das vantagens de se trabalhar com o RcppEigen e Rcpp, além da velocidade, é o fato de que não é necessário implementar todas as funções de interesse em C++. Podemos importar funções já existentes do R para que essas sejam executadas no ambiente Rcpp.

Por exemplo, um problema comum na formação de portfólios é o fato das matrizes de correlações estimadas entre os retornos dos ativos não ser Positivas Definidas. Nesses casos, problemas teóricos e numéricos podem surgir impedindo a obtenção de uma solução ótima para o problema de alocação de portfólios. Higham, N. (2002) propôs uma abordagem numérica para se encontrar a matriz Positiva Definida mais próxima da matriz de ineteresse.

Essa função está disponível no pacote Matrix e é denominada nearPD. Podemos importá-la no ambiente Rcpp da seguinte forma:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

// [[Rcpp::export]] Eigen::MatrixXd nearPDefinite(Eigen::MatrixXd X, int maxit=1e+6, double eigtol = 1e-06, double conv_tol = 1e-07, double posd_tol = 1e-08, bool keepDiagonal=false){   //Cria o objeto com o Environment.   Rcpp::Environment G = Rcpp::Environment::global_env();   //Objeto com as informações do pacote Matrix   Rcpp::Environment Matrix("package:Matrix");   //Importa a função de ineteresse   Function f = Matrix["nearPD"];   //Obtêm o resultado da função nearPD   Rcpp::List res = f(X, false, keepDiagonal, true, false, true, false, false,  eigtol, conv_tol, posd_tol, 100, "I", false);   //Pega o primeiro elemento da lista e transforma para matrix   Rcpp::NumericMatrix mat = internal::convert_using_rfunction(wrap(res[0]), "as.matrix");   //Converte de Rcpp::NumericMatrix para Eigen::MatrixXd   Eigen::MatrixXd Xpd=convertMatrix<eigen::matrixxd,rcpp::numericmatrix>(mat);   //Retorna o resultado   return(Xpd); } </eigen::matrixxd,rcpp::numericmatrix>

A função anterior, após compilada, pode ser utilizada no ambiente R de maneira direta pelo nome nearPDefinite. Entretanto, é necessário que a biblioteca Matrix tenha sido habilitada e esteja disponível no environment do R para que seja utilizada. Essa estratégia pode ser para outras funções de outros pacotes também.

RcppEigen: Operações Matriciais e Rcpp - Parte 1

Uma das grandes críticas que os usuários do R fazem é que por vezes algumas operações são muito demoradas no R. Em especial, a execução de loops.

Há, no entanto alguns truques que podem ser utilizados para se aumentar a velocidade dos códigos produzidos em R, quais sejam: usar operações algébricas vetoriais ao invés de loops, utilizar os membros da família apply, utilizar processamento em paralelo , etc.

Entretanto, algumas vezes é necessário programar funções mais complexas e que necessitam do uso de loops clássicos. Quando esse for o caso, uma boa solução é utilizar o pacote RcppEigen que é uma versão em R da biblioteca em C++ Eigen.

Nesse grupo de posts acerca do RcppEigen vamos mostrar como esse pacote pode ser usado e também como novos pacotes para o R podem ser criados. O primeiro passo é construir um projeto de pacote no RStudio:

Definido a criação de um novo pacote o próximo passo é determinar:

Podemos criar um projeto com base em um diretório já existente ou em um novo diretório. Escolhendo um novo diretório temos as seguintes opções:

Escolhemos nesse caso a opção R Package. Em seguida, definimos o local e nome do pacote que desejamos criar:

Automaticamente o R cria alguns arquivos e pastas. Para usar o RcppEigen precisamos de alguns arquivos específicos. Esses arquivos são obtidos executando uma função padrão do pacote:

?

1 2 3 4 5 6 7

#Limpa o Workspace #Habilita o pacote library(RcppEigen) #Define o local temporário para criação dos arquivos setwd("C:\\LocalTemorario") #Cria os arquivos necessários RcppEigen.package.skeleton("ArquivosRcppEigen")

Na pasta definida pelo comando setwd um conjunto de arquivos é criado.

Deve-se copiar os arquivos DESCRIPTION e NAMESPACE e também a pasta src (todas elas dentro da pasta "ArquivosRcppEigen") para a pasta do pacote criado no RStudio (nesse exemplo denominada "NomeDoPacote"). Após essas etapas o pacote estará pronto para receber as funções escritas em RcppEigen, as quais veremos nas próximas partes desse tutorial.

taskscheduleR: Agendando tarefas no R.

A automação de tarefas é importante quando deseja-se executar um mesmo script de maneira regular. Nesse sentido o pacote taskscheduleR pode apoiar no agendamento de scripts feitos em R.

Como exemplo vamos supor que estamos interessados em decidir se devemos comprar ou vender uma ação no dia seguinte baseado no preço de fechamento dos dias anteriores usando um modelo de Máquinas de Suporte Vetorial.

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

#Limpa o Workspace rm(list=ls()) #Chama a biblioteca library(devtools) library(quantmod) library(kernlab) #Instala o pacote taskscheduleR devtools::install_github("jwijffels/taskscheduleR") #Chama a biblioteca taskscheduleR library(taskscheduleR) #Especifica as datas de interesse startDate = as.Date("2013-01-01") endDate = as.Date("2015-12-31") #Obtêm os dados do ativo PETR4 e PETR3 getSymbols("PETR4.SA", src="yahoo",from=startDate,to=endDate) #Faz o Candle Chart candleChart(PETR4.SA,theme='white', type='candles') reChart(major.ticks='months',subset='first 16 weeks')

O Candle Chart é apresentado a seguir:

Em seguida vamos treinar uma máquina para tentar apontar a direção do resultado futuro, isto é, se a ação irá apresentar um resultado positivo ou negativo:

?

1 2 3 4 5 6 7 8

#Retorno do Preço de Fechamento ret<-na.omit(ClCl(PETR4.SA)) #Direção do mercado y<-as.numeric(ifelse(ret<0,0,1)) #Cria a base de dados df<-data.frame("Y"=c(y[-1],NA),"X"=ret) #Renomeia as colunas colnames(df)<-c("Y","X")

Nesse bloco a ideia é tentar prever a direção $Y=sinal(r_{t})$ usando o retorno do período anterior $X=r_{t-1}$. Supondo que o processo de treinamento tenha sido realizado e que tenhamos descobertos $C=0.1$ e $\sigma=16$ para o modelo C-SVC e kernel Gaussiano, a máquina é construída fazendo-se:

?

1 2 3 4 5 6 7 8 9

#Constrói a máquina svm <- ksvm(Y~X,data=df, type="C-svc",C=0.1,             kpar=list(sigma=16),cross=3) #Faz a previsão pred<- kernlab::predict(svm,df) #Compara com os valores reais table(pred,df$Y) #Salva a máquina treinada save.image("Maquina.RData")

Supondo que essa acurácia seja a desejada, agora vamos programar para que o R todo dia baixe a cotação do dia anterior e faça a previsão de qual deveria ser nossa ação no dia atual (vender se a previsão é que o retorno seja positivo ou comprar se a previsão seja do retorno negativo). As opções do pacote taskscheduleR são: ‘ONCE’, ‘MONTHLY’, ‘WEEKLY’, ‘DAILY’, ‘HOURLY’, ‘MINUTE’, ‘ONLOGON’, ‘ONIDLE’. O primeiro passo é montar um script geral que será executado todo dia:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

#Chama a biblioteca library(quantmod) library(kernlab) #Lê a máquina criada load("Maquina.RData") #Obtêm a cotação atual getSymbols("PETR4.SA", src="yahoo") #Último retorno observado ret<-tail(ClCl(PETR4.SA),1) #Faz a previsao para o dia previsao<-kernlab::predict(svm,ret) #Monta o arquivo de saida fileConn<-file("Decisao.txt") writeLines(ifelse(previsao==1,"Vender","Comprar"), fileConn) close(fileConn)

Esse script deve ser salvo em um arquivo de extensão .R (por exemplo, ScriptPrevisao.R) e, no nosso exemplo, ele deverá ser executado todo dia as 8:00 AM. Após a execução, ele criará um arquivo denominado Decisao.txt que contêm a decisão do dia: Comprar ou Vender. O último passo é programar o computador para executar esse script todo dia no horário determinado, para isso, basta executarmos uma única vez no R:

?

1 2 3 4 5 6

#Faz o agendamento library(taskscheduleR) myscript <- system.file("saida", "ScriptPrevisao.R", package = "taskscheduleR") #O computador irá executar todo dia as 8:00AM taskscheduler_create(taskname = "DecisaoPETR4", rscript = myscript,                      schedule = "DAILY", starttime = "08:00")

Caso queira deletar o agendamento para que ele não se repita mais, basta fazer:

?

1 2 3

#Deleta o agendamento library(taskscheduleR) taskscheduler_delete(taskname = "DecisaoPETR4")

Introdução ao RCpp - Parte 4.

Frequentemente é necessário inicializar um NumericVector com algum valor de constante fixo, nesse caso podemos usar a função fill(.):

?

1 2 3 4 5 6 7 8 9 10

#include <rcpp.h> using namespace Rcpp; // Função para repetir o valor no vetor // [[Rcpp::export]] NumericVector inicializaVetor(double valor, int tamanho) {   NumericVector x(tamanho);   x.fill(valor);   return x; } </rcpp.h>

Após compilar o código, a invocação no R é quase imediata:

?

1 2 3 4 5 6 7

#Invoca a biblioteca RCpp library(Rcpp) #Define o endereço do arquivo Cpp setwd("C:\\Blog\\Source") #Compila o código em C++ sourceCpp('Teste.cpp') inicializaVetor(2,5)

O pacote Rcpp tem muitas funcionalidades. Por meio do Rcpp podemos por exemplo, usar funções já implementadas no R, como:

?

1 2 3 4 5 6

// [[Rcpp::export]] NumericVector calculaCurtose(NumericVector x) {   Rcpp::Environment funcoesPacote("package:e1071");   Rcpp::Function funcaoCurtose = funcoesPacote["kurtosis"];   return funcaoCurtose (x); }

Ou por exemplo, gerar números com a distribuição normal:

?

1 2 3 4 5 6

// [[Rcpp::export]] NumericVector geraNormal(){   Rcpp::Environment stats("package:stats");   Rcpp::Function rnorm = stats["rnorm"];   return rnorm(10, Rcpp::Named("sd", 100.0)); }

Comandos esses que após compilados são invocados na forma:

?

1 2 3 4 5 6 7 8 9 10

#Invoca a biblioteca RCpp library(Rcpp) #Define o endereço do arquivo Cpp setwd("C:\\Blog\\Source") #Compila o código em C++ sourceCpp('Teste.cpp') #Executa as funções: x<-rnorm(1919) calculaCurtose(x) geraNormal()

Entre os objetos mais flexíveis do Rcpp está o objeto GenericVector o qual funciona como uma lista permitindo a inserção de qualquer outro objeto em seus slots. Além do mais, é possível usar os seguintes métodos nesse objeto:

• push_back(x) - Insere o elemento x no final da lista.
• push_front(x) - Insere o elemento x no início da lista.
• insert(i, x) - Insere o elemento x na posição i.
• erase(i) - Apaga o elemento existente na posição i.

Todas essas funções são dinâmicas e não há a necessidade de criar uma nova lista para que os novos elementos sejam inseridos.

Introdução ao RCpp - Parte 3.

Aqui neste post continuaremos com a análise dos possíveis métodos desenvovidos usando Rcpp. O primeiro exemplo trata do uso de matrizes, especificamente, vamos criar uma função para fazer a transposição de uma matriz:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

#include <rcpp.h> using namespace Rcpp; // Função para transposição de matriz // [[Rcpp::export]] NumericMatrix transpose(NumericMatrix x) {   //Obtêm o número de linhas e colunas:   int nrow = x.nrow(), ncol = x.ncol();   //Cria a matriz com resultado   NumericMatrix resultado(ncol,nrow);   //Faz a transposição:   for(int i=0;i < nrow; i++){     for(int j=0;j < ncol; j++){       resultado(j,i)=x(i,j);     }   }   return resultado; } </rcpp.h>

Após compilar o código, a invocação no R é quase imediata:

?

1 2 3 4 5 6 7 8 9 10

#Invoca a biblioteca RCpp library(Rcpp) #Define o endereço do arquivo Cpp setwd("C:\\Blog\\Source") #Compila o código em C++ sourceCpp('Teste.cpp') #Cria uma matriz de exemplo mat<-matrix(c(2,8,9,0,0,0,4,3,2), nrow=3,ncol=3) #Transpoem a matriz transpose(mat)

Métodos mais sofisticados usando funções implementadas em R e invocadas em C++ também podem ser utilizadas:

?

1 2 3 4 5 6

//Usando funções do R: // [[Rcpp::export]] NumericVector ExemploKernel(NumericVector x, NumericVector y, Function f) {   NumericVector res = f(x,y);   return res; }

Método esse que recebe como argumento dois vetores numéricos e uma função bivariada:

?

1 2 3 4 5 6 7 8 9

#Cria uma função no R gaussian<-function(x,y){   return(sum((x-y)^2)) } #Gera dois vetores x<-rnorm(1000) y<-rnorm(1000) #Invoca a função em Rcpp: ExemploKernel(x,y,gaussian)

Introdução ao RCpp - Parte 2.

Dando continuidade ao post de Introdução ao RCpp - Parte 1. Aqui eu mostrarei como os objetos do Rcpp se relacionam com os objetos em R. Os principais objetos do Rcpp são:

1. NumericVector - Objeto representando um vetor numérico no R.
2. CharacterVector - Objeto representando um vetor caracter no R.
3. NumericMatrix - Objeto representando uma matriz numérica no R.
4. CharacterMatrix - Objeto representando uma matriz caracter no R.
5. LogicalVector - Objeto representando um vetor lógico no R.
6. IntegerVector - Objeto representando um vetor inteiro no R.
7. LogicalMatrix - Objeto representando uma matriz lógica no R.
8. IntegerMatrix - Objeto representando uma matriz inteira no R.
9. List - Objeto representando uma lista no R.

Aqui trabalharemos com exemplos de vetores e nos próximos posts o uso de matrizes e listas. Suponha que temos um vetor de retornos e um vetor com os nomes dos ativos representando esses vetores, algo como:

?

1 2 3 4 5 6 7 8 9 10

#Invoca a biblioteca RCpp library(Rcpp) #Define o endereço do arquivo Cpp setwd("C:\\Blog\\Source") #Define a semente para simulação dos dados set.seed(19393) #Gera um vetor numérico de retornos: ret<-rnorm(1000, 0,0.1) #Gera um vetor caracter de nomes: nomes<-replicate(1000, paste(sample(LETTERS, 5, replace=TRUE), collapse=""))

Com base nesses dois inputs podemos pensar em algumas funções escritas em Rcpp que utilizem esses insumos. Por exemplo, considere que desejamos obter os 5 menores valores do vetor ret e os respectivos nomes dessas empresas.

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

//Exemplo de função para obter os k menores valores. // [[Rcpp::export]] NumericVector kMenores(int k, NumericVector ret, CharacterVector nomes) {   //Cria um vetor numérico de tamanho k   NumericVector vetMenoresRetornos(k);   //Cria um vetor caracter de tamanho k   CharacterVector vetMenoresNomes(k);   //Faz o loop para achar os k menores   for (int i = 0; i < ret.size(); i++){     //Verifica qual menor é:     int iMenor = qualMenor(ret[i],vetMenoresRetornos);     if(iMenor>-1){       vetMenoresRetornos[iMenor]=ret[i];     }         }   return vetMenoresRetornos; }

Onde a função qualMenor é criada no topo do mesmo arquivo em C++ da seguinte forma:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14

#include <rcpp.h> using namespace Rcpp;   //Função para dizer qual posição esse menor é: int qualMenor(double valor,NumericVector vetMenoresRetornos){   //Inicializa a variável   int menor=-1;   //Procura qual o menor:   for(int i=0;i< vetMenoresRetornos.size();i++){     if(valor<=vetMenoresRetornos[i])  return i;   }   return(menor); } </rcpp.h>

Observe que as linhas 1 e 2 do código anterior são mandatórias e devem constar em todos os arquivo com extensão .cpp, ademais como o código anterior não contêm a linha // [[Rcpp::export]] isso significa que a função é interna e não pode ser invocada pelo R. Finalmente, a invocação no R é feito da seguinte forma:

?

1 2 3 4 5 6 7 8 9 10

#Invoca a biblioteca RCpp library(Rcpp) #Define o endereço do arquivo Cpp setwd("C:\\Blog\\Source") #Compila o código em C++ sourceCpp('Teste.cpp') #Inova a função compilada teste<-kMenores(5,ret,nomes) #Confere com a função do R head(sort(teste),5)

Repare que em nenhum momento utilizamos o objeto nomes. Para imprimir o nome das empresas com menor retorno no console podemos alterar a função kMenores da seguinte forma:

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

//Exemplo de função para obter os k menores valores. // [[Rcpp::export]] NumericVector kMenores(int k, NumericVector ret, CharacterVector nomes) {   //Cria um vetor numérico de tamanho k   NumericVector vetMenoresRetornos(k);   //Cria um vetor caracter de tamanho k   CharacterVector vetMenoresNomes(k);   //Faz o loop para achar os k menores   for (int i = 0; i < ret.size(); i++){     //Verifica qual menor é:     int iMenor = qualMenor(ret[i],vetMenoresRetornos);     if(iMenor>-1){       //Guarda os retornos mínimos       vetMenoresRetornos[iMenor]=ret[i];       //Guarda os nomes das firmas       vetMenoresNomes[iMenor]=nomes[i];     }   }   //Antes de retornar os resultados imprimir no console o nome   for(int j=0;j < k;j++){     std::cout << "Firma: " << vetMenoresNomes[j] << " Retorno: " << vetMenoresRetornos[j] << std::endl;   }   return vetMenoresRetornos; }

A invocação procede de maneira equivalente a apresentada anteriormente, a diferença é que agora no console do R o nome das firmas e seus retornos simulados são apresentados. Alguns links interessantes para mais detalhes sobre o RCpp:

1. http://adv-r.had.co.nz/Rcpp.html
2. http://gallery.rcpp.org/articles/r-function-from-c++/
3. http://pages.stat.wisc.edu/~gvludwig/Rcpp/

Introdução ao RCpp - Parte 1.

O R possui uma limitação devido a demora de alguns procedimentos realizados nele. O software foi criado para ser ótimo do ponto de vista de processamento computacional quando o usuário utiliza as informações de maneira vetorizada. Quando o processo é realizado de elemento a elemento alguma demora no processamento surge.

Nesse sentido, aqui daremos uma primeira introdução quanto ao uso do pacote RCpp (http://www.rcpp.org/), o qual permite otimizar alguns procedimentos em R usando linguagem C++.

O primeiro passo para desenvolver suas próprias funções em RCpp é ter instalado:

1. R
2. RTools - É preciso que ele esteja instalado na raiz (C:\) e com todas as opções habilitadas.
3. RStudio

Em seguida, abrimos o RStudio e escolhemos a opção File → New File → C++ File:

O RStudio fornecerá um exemplo de como construir funções usando RCpp:

Onde o primeiro bloco de código, da linha 1 até a 24 é código escrito em C++, o que está contido na linha 26 é o código em R usando para a invocação da função em C++. Para compilar a função é necessário clicar no botão Source destacado na figura anterior. Ao clicar nesse botão a função é compilada e fica disponível para ser executada.

É comum no entanto ter dois arquivos: um arquivo com as funções em C++ e outro arquivo com as funções em R, por isso considere o seguinte código para o arquivo com extensão .cpp:

?

1 2 3 4 5 6 7 8

#include <rcpp.h> using namespace Rcpp; // This is a simple example of exporting a C++ function to R. // [[Rcpp::export]] NumericVector timesTwo(NumericVector x) {   return x * 2; } </rcpp.h>

E no arquivo em R (com extensão .R):

?

1 2 3 4 5 6 7 8

#Invoca a biblioteca RCpp library(Rcpp) #Define o endereço do arquivo Cpp setwd("C:\\Blog\\Source") #Compila o código em C++ sourceCpp('Teste.cpp') #Inova a função compilada timesTwo(42)

Note que nesse exemplo, Teste.cpp é o nome do arquivo com o código em RCpp armazenado na pasta C:\Blog\Source. No próximo post mostrarei como fazer funções mais complexas envolvendo os pacotes do R e seus objetos com o Rcpp.