Abstract:In the contemporary information society, constructing an effective sales prediction model is challenging due to the sizeable amount of purchasing information obtained from diverse consumer preferences. Many empirical cases shown in the existing literature argue that the traditional forecasting methods, such as the index of smoothness, moving average, and time series, have lost their dominance of prediction accuracy when they are compared with modern forecasting approaches such as neural network (NN) and support vector machine (SVM) models. To verify these findings, this paper utilizes the Taiwanese cosmetic sales data to examine three forecasting models: i) the back propagation neural network (BPNN), ii) least-square support vector machine (LSSVM), and iii) auto regressive model (AR). The result concludes that the LS-SVM has the smallest mean absolute percent error (MAPE) and largest Pearson correlation coefficient ( R2 ) between model and predicted values.
Análise Preditiva
Prevendo recessões econômicas usando algoritmos de Machine Learning
Paper bem atual que fala como os autores erraram a crise apenas em relação ao ano mostrando o potencial das Random Forests.
Predicting Economic Recessions Using Machine Learning Algorithms – Rickard Nyman and Paul Ormerod
Implementação de GLM com Grid Search no R usando o H2O.ai como backend
Para quem usa o R não existe nada mais irritante do que ter que lidar com o péssimo gerenciamento de memória da ferramenta, o que limita e muito o uso do R como uma ferramenta séria para a construção de modelos que possam ir para produção e possam permitir a construção de plataformas/sistemas inteligentes.
Vamos aqui em algumas linhas mostrar como usar o H2O.ai como backend de processamento (o que abstraí todos esses problemas de memória e processamento) para a criação de um modelo usando GLM.
O pulo do gato aqui é que o H2O faz todo o gerenciamento de memória, e independente da sua fonte de dados ele faz todo o pipeline do buffer de memória de forma que não há estouro de memória; ou mesmo uma lentidão generalizada no sistema.
Esse exemplo é baseado totalmente na documentação do H2O e tem o objetivo somente de mostrar como essa ferramenta funciona.
Nesse caso eu vou usar em um notebook, mas poderia ser utilizado por exemplo em uma máquina na Amazon usando o comando abaixo no momento da inicialização do cluster:
#Production Cluster (Not applicable) #localH2O <- h2o.init(ip = '10.112.81.210', port =54321, nthreads=-1) # Máquina 1 #localH2O <- h2o.init(ip = '10.112.80.74', port =54321, nthreads=-1) # Máquina 2 - Sim, aqui usamos um cluster com dois nós para processamento! 😉
Primeiramente vamos remover qualquer instalação antiga do H2O.ai da máquina em questão:
# The following two commands remove any previously installed H2O packages for R.
if ("package:h2o" %in% search()) { detach("package:h2o", unload=TRUE) }
if ("h2o" %in% rownames(installed.packages())) { remove.packages("h2o") }
Em seguida vamos fazer o download e instalação de todos os pacotes dos quais o H2O tem alguma dependência direta ou indireta.
# Next, we download packages that H2O depends on.
if (! ("methods" %in% rownames(installed.packages()))) { install.packages("methods") }
if (! ("statmod" %in% rownames(installed.packages()))) { install.packages("statmod") }
if (! ("stats" %in% rownames(installed.packages()))) { install.packages("stats") }
if (! ("graphics" %in% rownames(installed.packages()))) { install.packages("graphics") }
if (! ("RCurl" %in% rownames(installed.packages()))) { install.packages("RCurl") }
if (! ("jsonlite" %in% rownames(installed.packages()))) { install.packages("jsonlite") }
if (! ("tools" %in% rownames(installed.packages()))) { install.packages("tools") }
if (! ("utils" %in% rownames(installed.packages()))) { install.packages("utils") }
Em seguida faremos a instalação da lib do H2O e o instanciamento da lib no R Studio.
# Now we download, install and initialize the H2O package for R.
install.packages("h2o", type="source", repos=(c("http://h2o-release.s3.amazonaws.com/h2o/rel-turing/8/R")))
# Load library library(h2o)
Instalação feita e biblioteca carregada, vamos agora para algumas configurações.
No próprio R Studio você pode escolher o número de processadores no qual o cluster (nesse caso o seu notebook/desktop) vai utilizar. Lembrando que quanto maior for o número de cores utilizados, mais processamento o H2O vai consumir e menos recursos estarão disponíveis para outras tarefas. O padrão é a utilização de 2 cores, mas no meu caso eu vou usar todos os processadores.
# Start instance with all cores. # The -1 is the parameter to use with all cores. Use this carefully. # The default parameter is 2 cores. h2o.init(nthreads = -1)
Agora vamos ver as informações do nosso cluster:
# Cluster Info h2o.clusterInfo() # R is connected to the H2O cluster: # H2O cluster uptime: 3 seconds 267 milliseconds # H2O cluster version: 3.10.0.8 # H2O cluster version age: 2 months and 26 days # H2O cluster name: H2O_started_from_R_flavio.clesio_udy929 # H2O cluster total nodes: 1 # H2O cluster total memory: 1.78 GB # H2O cluster total cores: 4 # H2O cluster allowed cores: 4 # H2O cluster healthy: TRUE # H2O Connection ip: localhost # H2O Connection port: 54321 # H2O Connection proxy: NA # R Version: R version 3.3.2 (2016-10-31)
Como podemos ver dos 4 processadores no total, estou usando todos eles (allowed cores) para o processamento.
Outro fato que podemos ver aqui é o que o H2O também está instanciado para usar a GUI na Web. Para isso, basta entrar no endereço no navegador com o endereço http://localhost:54321/flow/index.html.
Para este exemplo, vamos usar a base de dados Airlines que contém diversas informações reais de voos nos EUA e todas as causas de atraso de 1987 até 2008. A versão completa com 12Gb pode ser encontrada aqui.
Seguindo adiante, vamos agora fazer o carregamento dos dados direto de uma URL e importar em um objeto do R.
# GLM Demo Deep Dive # Path of normalized archive. Can be a URL or a local path airlinesURL = "https://s3.amazonaws.com/h2o-airlines-unpacked/allyears2k.csv"
# We'll create the object .hex (extention of data files in H2O) # and using the importFile property, we'll set the path and the destination frame. # As default behaviour H2O parse the file automatically. airlines.hex = h2o.importFile(path = airlinesURL, destination_frame = "airlines.hex")
Neste caso o objeto airlines.hex é será o dataframe no qual o H2O irá aplicar os algoritmos.
Esse formato .hex é exclusivo do H2O e pode ser usado para inúmeros algoritmos dentro da plataforma, dado que ele já é otimizado para lidar com objetos esparsos e/ou colunas do tipo texto.
Para ver as estatísticas descritivas desse arquivo, basta usar o mesmo summary() do R.
# Let's see the summary summary(airlines.hex)
Para o nosso experimento, vamos dividir a base de treino e teste na proporção 70%/30%.
Uma coisa necessária a se dizer nesse ponto é que devido ao fato do H2O ser uma plataforma projetada para Big Data é utilizado o método de amostragem probabilística. Isso se faz necessário (em termos computacionais), dado que em muitas vezes a operação de seleção/estratificação pode ser custoso.
# Construct test and train sets using sampling # A small note is that H2O uses probabilistic splitting, witch means that resulting splits # can deviate for the exact number. This is necessary when we're talking about a platform that # deals with big data. If you need a exact sampling, the best way is to do this in your RDBMS airlines.split = h2o.splitFrame(data = airlines.hex,ratios = 0.70, seed = -1)
Após criar o objeto do tipo splitFrame, vamos alocar as partições para cada conjunto de dados, sendo que o objeto na primeira posição sempre será a nossa base de treinamento, e na segunda posição a nossa base de teste.
# Get the train dataframe(1st split object) airlines.train = airlines.split[[1]] # Get the test dataframe(2nd split object) airlines.test = airlines.split[[2]]
Vamos sumarizar abaixo cada um desses frames para verificarmos a distribuição de voos cancelados:
# Display a summary using table-like in some sumarized way h2o.table(airlines.train$Cancelled) # Cancelled Count # 1 0 29921 # 2 1 751 h2o.table(airlines.test$Cancelled) # Cancelled Count # 1 0 12971 # 2 1 335
Com as nossas amostras separadas, vamos agora escolher as variáveis que vão entrar no nosso modelo.
Primeiramente, vamos criar dois objetos para passar como parâmetro ao nosso algoritmo.
Como queremos prever se as partidas do voos estão atrasadas, então o objeto Y (variável dependente) será a variável IsDepDelayed (Se o voo de partida está atrasado) e o objeto X (variáveis independentes) serão todos os outros campos do conjunto de dados.
# Set dependent variable (Is departure delayed) Y = "IsDepDelayed"
# Set independent variables
X = c("Origin", "Dest", "DayofMonth", "Year", "UniqueCarrier", "DayOfWeek", "Month", "DepTime", "ArrTime", "Distance")
Agora vamos realizar a criação do modelo usando GLM:
# Define the data for the model and display the results
airlines.glm <- h2o.glm(training_frame=airlines.train
,x=X
,y=Y
,family = "binomial"
,alpha = 0.5
,max_iterations = 300
,beta_epsilon = 0
,lambda = 1e-05
,lambda_search = FALSE
,early_stopping = FALSE
,nfolds = 0
,seed = NULL
,intercept = TRUE
,gradient_epsilon = -1
,remove_collinear_columns = FALSE
,max_runtime_secs = 10000
,missing_values_handling = c("Skip"))
O significado dos parâmetros do modelo são:
x: vetor que contém os nomes das variáveis independentes;
y: índice que contém a variável dependente;
training_frame: Um frame do H2O que contém das variáveis do modelo;
family: Especificação da distribuição do modelo que pode ser gaussiana, binomial, poisson, gamma, e tweedie. Uma ótima explicação de como esses parâmetros podem ser escolhidos está aqui nesse link;
alpha: Um número em [0, 1] especificando a mistura do parâmetro de regularização do elastic-net. Ele que dá o grau de mistura entre os regularizadores Lasso e Ridge. making alpha = 1 penalização via LASSO, alpha = 0 penalização via ridge;
max_iterations: Um inteiro não negativo que especifica o número máximo de interações do modelo;
beta_epsilon: Um inteiro não negativo que especifica a magnitude da diferença máxima entre as estimativas dos coeficientes através de sucessivas interações. Em outras palavras: É esse parâmetro que define a velocidade da convergência do modelo GLM;
lambda: Um parâmetro não negativo para encolhimento do valor da variável através da Elastic-Net, o qual multiplica P(α, β) na função objetivo. Quando lambda = 0, nenhuma penalização é aplicada e o modelo já fica ajustado;
lambda_search: Um valor lógico que indica se haverá algum critério de busca no espaço dos valores de lambda definidos por um parâmetro de mínimo e máximo;
early_stopping: Um valor que indica se haverá uma parada antecipada durante o lambda_search caso o fator de verosimilhança pare de ser alterado na medida que ocorram mais interações;
nfolds: Número de partições em Cross-Validation;
seed: Especifica a semente do random number generator (RNG) para Cross-Validation (garante a reprodutibilidade do experimento);
intercept: Termo constante do modelo que a grosso modo significa o grau de fatores endógenos do modelo;
gradient_epsilon: Critério de convergência. Converge se o gradiente da norma I-Infinito é abaixo de um determinado limite. Se lambda_search = FALSE e lambda = 0, o valor default do gradient_epsilon é igual a .000001, se não for, o valor default é .0001. Se lambda_search = TRUE, os valores condicionais acima são 1E-8 e 1E-6 respectivamente.
remove_collinear_columns: Se não houver nenhum tipo de fator de regularização aplicado, o modelo ignora colunas colineares (o coeficiente será 0);
max_runtime_secs: Número máximo permitido em segundos para o treinamento do modelo. Use 0 para desabilitar; e
missing_values_handling: Contra o que é feito com valores faltantes. Podem ser “MeanImputation” ou “Skip”. MeanImputation substituí os valores faltantes com a média para os atributos numéricos e categórico com a maior frequência. É aplicado durante o treinamento do modelo.
Notem que aqui o céu é o limite em temos de ajustes e/ou parametrizações. O ideal é ter o perfeito entendimento da mecânica de cada um dos parâmetros e utilizar a melhor combinação possível.
Com o nosso modelo ajustado, vamos ver algumas das estatísticas básicas desse modelo.
# View model information: training statistics, performance, important variables summary(airlines.glm) # Model Details: # ============== # # H2OBinomialModel: glm # Model Key: GLM_model_R_1484053333586_1 # GLM Model: summary # family link regularization number_of_predictors_total number_of_active_predictors number_of_iterations training_frame # 1 binomial logit Elastic Net (alpha = 0.5, lambda = 1.0E-5 ) 283 272 5 RTMP_sid_a6c9_1 # # H2OBinomialMetrics: glm # ** Reported on training data. ** # # MSE: 0.2098326 # RMSE: 0.4580749 # LogLoss: 0.607572 # Mean Per-Class Error: 0.3720209 # AUC: 0.7316312 # Gini: 0.4632623 # R^2: 0.1602123 # Null Deviance: 41328.6 # Residual Deviance: 36240.45 # AIC: 36786.45 # # Confusion Matrix for F1-optimal threshold: # NO YES Error Rate # NO 5418 9146 0.627987 =9146/14564 # YES 1771 13489 0.116055 =1771/15260 # Totals 7189 22635 0.366047 =10917/29824 # # Maximum Metrics: Maximum metrics at their respective thresholds # metric threshold value idx # 1 max f1 0.363651 0.711915 294 # 2 max f2 0.085680 0.840380 389 # 3 max f0point5 0.539735 0.683924 196 # 4 max accuracy 0.521518 0.673887 207 # 5 max precision 0.987571 1.000000 0 # 6 max recall 0.040200 1.000000 398 # 7 max specificity 0.987571 1.000000 0 # 8 max absolute_mcc 0.521518 0.348709 207 # 9 max min_per_class_accuracy 0.513103 0.672412 212 # 10 max mean_per_class_accuracy 0.521518 0.674326 207
Aqui neste modelo já temos um resultado de 73,16% de AUC. Nada mal para um modelo que contém poucos ajustes.
Vamos analisar agora a importância de cada uma das variáveis no modelo:
# Get the variable importance of the models h2o.varimp(airlines.glm) # Standardized Coefficient Magnitudes: standardized coefficient magnitudes # names coefficients sign # 1 Origin.TLH 3.233673 NEG # 2 Origin.CRP 2.998012 NEG # 3 Origin.LIH 2.859198 NEG # 4 Dest.LYH 2.766090 POS # 5 Origin.KOA 2.461819 NEG # # --- # names coefficients sign # 278 Dest.JAN 0.000000 POS # 279 Dest.LIT 0.000000 POS # 280 Dest.SJU 0.000000 POS # 281 Origin.LAN 0.000000 POS # 282 Origin.SBN 0.000000 POS # 283 Origin.SDF 0.000000 POS
Alguns valores de atributos são determinantes no atraso dos vôos de partida, principalmente se os aeroportos de origem são TLH (Tallahassee International Airport), CRP (Corpus Christi International Airport), LIH (Lihue Airport), e KOA (Kona International Airport).
Agora vamos usar a função predict para saber como o modelo realizou as classificações.
# Predict using GLM model pred = h2o.predict(object = airlines.glm, newdata = airlines.test)
Após isso, vamos ver o resultado do nosso modelo.
# Look at summary of predictions: probability of TRUE class (p1) summary(pred) # predict NO YES # YES:9798 Min. :0.01186 Min. :0.02857 # NO :3126 1st Qu.:0.33715 1st Qu.:0.37018 # NA : 382 Median :0.48541 Median :0.51363 # Mean :0.48780 Mean :0.51220 # 3rd Qu.:0.62886 3rd Qu.:0.66189 # Max. :0.97143 Max. :0.98814 # NA's :382 NA's :382
Outra forma de usar o modelo GLM no H2O é realizar a escolha de parâmetros via Grid Search.
Grid Search nada mais é do que um método de otimização de parâmetros de um modelo de Machine Learning, sendo que em grande parte das vezes é feita uma combinação com inúmeros parâmetros e a combinação que obtiver o menor erro através de uma determinada função de erro.
O que vamos fazer agora é usar o GLM pra obter o melhor modelo de acordo com uma combinação de parâmetros específica.
Primeiramente, vamos construir uma lista com os parâmetros alpha (recapitulando, esse parâmetro faz uma combinação de Lasso e Ridge via Elastic-Net). Neste caso vamos passar uma lista que vai desde 0.0 até 1 com incremento de 0.05 em cada parâmetro.
# Construct a hyper-parameter space alpha_opts = c(0,0.05,0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1)
Agora vamos criar uma lista com esses parâmetros para passar para a função de Grid posteriormente.
# List of hyperparameters hyper_params_opt = list(alpha = alpha_opts)
Na função de Grid Search, vamos passar alguns parâmetros para o ajuste dos modelos como podemos ver abaixo.
# Grid object with hyperparameters
glm_grid <- h2o.grid("glm"
,grid_id = "glm_grid_1"
,x=X
,y=Y
,training_frame=airlines.train
,hyper_params = hyper_params_opt
,family = "binomial")
Essa etapa pode demorar bastante tempo dependendo do seu volume de dados, e do número de parâmetros escolhidos na lista de search.
Após a finalização do processamento, vamos ordernar a lista de modelos de acordo com o AUC.
# Sort grids by best performance (lower AUC). Little note: As we're dealing with classification # in some probabilistc fashion, we'll use AUC as model selection metric. # If the nature of the problem are cost sensitive (e.g. A delayed departure plane is much expensive for # the airport service than a delayed arrival) precision and recall can be the best choice glm_sorted_grid <- h2o.getGrid(grid_id = "glm_grid_1", sort_by = "auc", decreasing = FALSE)
Para avaliar cada um dos modelos, podemos exibir a ordem dos modelos de acordo com o AUC.
#Print the models print(glm_sorted_grid) # H2O Grid Details # ================ # # Grid ID: glm_grid_1 # Used hyper parameters: # - alpha # Number of models: 21 # Number of failed models: 0 # # Hyper-Parameter Search Summary: ordered by increasing auc # alpha model_ids auc # 1 [D@4800a43e glm_grid_1_model_1 0.7076911403181928 # 2 [D@66030470 glm_grid_1_model_2 0.7122987232329416 # 3 [D@6a4a43d3 glm_grid_1_model_3 0.7145455620514375 # 4 [D@17604a1a glm_grid_1_model_4 0.715989429818657 # 5 [D@21e1e99f glm_grid_1_model_5 0.7169797604977775 # # --- # alpha model_ids auc # 16 [D@78833412 glm_grid_1_model_16 0.720595118360825 # 17 [D@44d770f2 glm_grid_1_model_17 0.7207086912177467 # 18 [D@31669527 glm_grid_1_model_18 0.7208228330257134 # 19 [D@5b376f34 glm_grid_1_model_19 0.7209144533220885 # 20 [D@6acad45e glm_grid_1_model_20 0.7209885192412766 # 21 [D@237ad7de glm_grid_1_model_0 0.7240682725570593
Com esses parâmetros, o melhor modelo é o glm_grid_1_model_0 que teve cerca de 72.40% de AUC. (Nota: Esse modelo está levemente pior do que o modelo padrão, dado que o conjunto de parâmetros do Grid está diferente do que o primeiro modelo).
Para pegar o melhor modelo, basta executar o comando abaixo:
# Grab the model_id based in AUC best_glm_model_id <- glm_grid@model_ids[[1]]
# The best model best_glm <- h2o.getModel(best_glm_model_id)
Vejamos as características desse modelo:
# Summary summary(best_glm) # Model Details: # ============== # # H2OBinomialModel: glm # Model Key: glm_grid_1_model_0 # GLM Model: summary # family link regularization number_of_predictors_total number_of_active_predictors number_of_iterations training_frame # 1 binomial logit Ridge ( lambda = 7.29E-5 ) 283 282 3 RTMP_sid_a6c9_1 # # H2OBinomialMetrics: glm # ** Reported on training data. ** # # MSE: 0.2121424 # RMSE: 0.4605891 # LogLoss: 0.612699 # Mean Per-Class Error: 0.3833898 # AUC: 0.7240683 # Gini: 0.4481365 # R^2: 0.1494395 # Null Deviance: 42448.59 # Residual Deviance: 37585.41 # AIC: 38151.41 # # Confusion Matrix for F1-optimal threshold: # NO YES Error Rate # NO 4993 9601 0.657873 =9601/14594 # YES 1751 14327 0.108907 =1751/16078 # Totals 6744 23928 0.370110 =11352/30672 # # Maximum Metrics: Maximum metrics at their respective thresholds # metric threshold value idx # 1 max f1 0.373247 0.716243 296 # 2 max f2 0.105583 0.846435 391 # 3 max f0point5 0.551991 0.685249 194 # 4 max accuracy 0.513313 0.665949 218 # 5 max precision 0.980714 1.000000 0 # 6 max recall 0.048978 1.000000 399 # 7 max specificity 0.980714 1.000000 0 # 8 max absolute_mcc 0.548278 0.332916 196 # 9 max min_per_class_accuracy 0.524282 0.664324 211 # 10 max mean_per_class_accuracy 0.548278 0.666166 196 # # Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)` # # # # Scoring History: # timestamp duration iteration negative_log_likelihood objective # 1 2017-01-10 11:11:07 0.000 sec 0 21224.29620 0.69198 # 2 2017-01-10 11:11:07 0.066 sec 1 18857.11178 0.61705 # 3 2017-01-10 11:11:07 0.094 sec 2 18795.11788 0.61562 # 4 2017-01-10 11:11:07 0.126 sec 3 18792.70362 0.61559 # # Variable Importances: (Extract with `h2o.varimp`) # ================================================= # # Standardized Coefficient Magnitudes: standardized coefficient magnitudes # names coefficients sign # 1 Origin.MDW 1.915481 POS # 2 Origin.HNL 1.709757 NEG # 3 Origin.LIH 1.584259 NEG # 4 Origin.HPN 1.476562 POS # 5 Origin.AUS 1.439134 NEG # # --- # names coefficients sign # 278 Origin.PHX 0.009111 POS # 279 Dest.PWM 0.008332 POS # 280 Origin.GEG 0.008087 POS # 281 Dest.BOS 0.005105 POS # 282 Dest.MCI 0.003921 NEG # 283 Dest.CHA 0.000000 POS
Para realizar previsões com esse modelo, basta apenas instanciar esse novo objeto e usar a função predict como está abaixo:
# Get model and put inside a object model = best_glm # Prediction using the best model pred2 = h2o.predict(object = model, newdata = airlines.test) # Summary of the best model summary(pred2) # predict NO YES # YES:10368 Min. :0.01708 Min. :0.05032 # NO : 2938 1st Qu.:0.33510 1st Qu.:0.39258 # Median :0.47126 Median :0.52781 # Mean :0.47526 Mean :0.52474 # 3rd Qu.:0.60648 3rd Qu.:0.66397 # Max. :0.94968 Max. :0.98292
Se após isso, você quiser desligar o cluster do H2O basta usar o comando shutdown.
# Shutdown the cluster h2o.shutdown() # Are you sure you want to shutdown the H2O instance running at http://localhost:54321/ (Y/N)? Y # [1] TRUE
Com isso finalizamos esse post/tutorial de como usar o R com o H2O.ai.
Ao longo das próximas semanas, vamos trazer alguns tutoriais e destrinchar um pouco o poder desse H2O.
Forte abraço!
The Predictron: End-To-End Learning and Planning
Um review do resultado em relação à arquitetura:
The predictron is a single differentiable architecture that rolls forward an internal model to estimate values. This internal model may be given both the structure and the semantics of traditional reinforcement learning models. But unlike most approaches to model-based reinforcement learning, the model is fully abstract: it need not correspond to the real environment in any human understandable fashion, so long as its rolled-forward “plans” accurately predict outcomes in the true environment.
The predictron may be viewed as a novel network architecture that incorporates several separable ideas. First, the predictron outputs a value by accumulating rewards over a series of internal planning steps. Second, each forward pass of the predictron outputs values at multiple planning depths. Third, these values may be combined together, also within a single forward pass, to output an overall ensemble value. Finally, the different values output by the predictron may be encouraged to be self-consistent with each other, to provide an additional signal during learning. Our experiments demonstrate that these differences result in more accurate predictions of value, in reinforcement learning environments, than more conventional network architectures.
We have focused on value prediction tasks in uncontrolled environments. However, these ideas may transfer to the control setting, for example by using the predictron as a Q-network (Mnih et al., 2015). Even more intriguing is the possibility of learning an internal MDP with abstract internal actions, rather than the MRP model considered in this paper. We aim to explore these ideas in future work.
Porque intervalos de confiança em previsões de séries temporais não são boas quanto desejamos?
Para quem trabalha com modelos tipicamente de previsão usando ARIMA e Auto-ARIMA um aspecto bem difícil de se estimar é a a incerteza dos termos auto-regressivos (isso é, os seus intervalos de erros).
Na prática, no momento em que temos os termos dos coeficientes auto regressores absorvendo parte da incerteza seja por conta de meta-parametrização (ou a falta dela) ou mesmo devido à natureza dos dados, os modelos de previsão de séries temporais não conseguem captar esse tipo de incerteza, e nesse caso acontecem os problemas dos intervalos de confiança não representarem exatamente um range aceitável/factível.
The problem is that for all but the most trivial time series forecasting method there is no simple way of estimating the uncertainty that comes from having estimated the parameters from the data, and much less so the values of meta-parameters like the amount of differencing needed, how many autoregressive terms, how many moving average terms, etc (those example meta-parameters come from the Box-Jenkins ARIMA approach, but other forecasting methods have their own meta-parameters to estimate too).
Previsão de Séries Temporais usando XGBoost – Pacote forecastxgb
Para quem já teve a oportunidade de trabalhar com previsão de variáveis categóricas em Machine Learning sabe que o XGBoost é um dos melhores pacotes do mercado, sendo largamente utilizado em inúmeras competições no Kaggle.
A grande diferença feita pelo Peter Ellis foi realizar algumas adaptações para incorporar algumas variáveis independentes através do parâmetro xreg ao modelo preditivo de séries temporais.
Para quem trabalha com análise de séries temporais, esse trabalho é muito importante até porque o forecastxgb tríade Média-Móvel/ARIMA (ARMA)/(S)ARIMA em que tanto estatísticos/Data Miners/Data Scientists ficam presos por comodidade ou falta de meios.
Um exemplo da utilização do pacote está abaixo:
# Install devtools to install packages that aren't in CRAN
install.packages("devtools")
# Installing package from github
devtools::install_github("ellisp/forecastxgb-r-package/pkg")
# Load the libary
library(forecastxgb)
# Time Series Example
gas
# Model
model <- xgbts(gas)
# Summary of the model
summary(model)
# Forecasting 12 periods
fc <- forecast(model, h = 12)
# Plot
plot(fc)
Previsão de Churn no Warcraft
Direto do Dataiku
Previsão do número de homicídios no Rio de Janeiro (Sim, eu voltei a postar!!!)
Um pequeno disclaimer: Depois de muito relutar, estou voltando com as atividades do blog. Isso significa de maneira prática, que apesar do clipping ser o grande driver dos acessos, sinto que há uma lacuna em relação a posts técnicos aqui na blogosfera brasileira.
Dessa forma, vou fazer mais posts técnicos com assuntos ligados à Machine Learning, Data Mining e afins.
Disclaimer feito, bora pro post original.
*******
Fala pessoal, tudo bem?
Um dos principais problemas em utilizar processadores numéricos como o Excel ou mesmo o LibreOffice é que se perde um pouco a flexibilidade de parametrizar os modelos de previsão.
Além disso, há um problema grave nesses processadores de reprodutibilidade, dado que pode haver inúmeras versões, planilhas voando pra todo o lado via e-mail, e por aí vai.
Nesse post eu vou utilizar dados reais da quantidade de homicídios que houveram no Rio de Janeiro de 1991 até 2011.
Os dados foram extraídos e consolidados pelo CESeC – Centro de Estudos de Segurança e Cidadania e o link está no final do script.
Antes de mais nada, vamos carregar o package Forecast:
# Load library
if(!require(forecast)){
install.packages("forecast")
library(forecast)
}
Depois disso, vamos realizar a carga dos dados. Nesse caso eu fiz uma transformação na planilha original e passei os dados para csv.
# Dataset
rj_murder <- read.csv("https://raw.githubusercontent.com/fclesio/learning-space/master/Datasets/rj-homicidios-1991-2011.csv")
Vamos dar uma olhada nos dados para nos certificar que tudo está OK.
# Check data head(rj_murder)
Com os dados OK, para incluir um conjunto de dados no package de Forecast e ARIMA, temos que antes de mais nada, pegar todos os datapoints e construir uma série através de um objeto de timeseries como no bloco de código abaixo:
# Transform in time series dataset all metrics inside the dataset rj_murder_total <- ts(rj_murder$total, start=c(1991, 1), end=c(2011, 12), frequency=12) rj_murder_capital <- ts(rj_murder$capital, start=c(1991, 1), end=c(2011, 12), frequency=12) rj_murder_baixada <- ts(rj_murder$baixada, start=c(1991, 1), end=c(2011, 12), frequency=12) rj_murder_interior <- ts(rj_murder$interior, start=c(1991, 1), end=c(2011, 12), frequency=12)
Para entendermos o que aconteceu, vamos destrinchar erssa função TS:
ts() função para a construção de série temporal
rj_murder$total conjunto de datapoints do dataframe da coluna total
start=c(1991, 1) início da série é janeiro de 1991
end=c(2011, 12) o fim da série é dezebro de 2011
frequency=12 como estamos falando de dados mensais, a frequencia da série é de 12 datapoints por ano
Agora que entendemos o que fizemos, vamos dar uma olhada nos nossos objetos ts() e ver se os mesmos estão corretos.
#Plot series par(mfrow=c(2,2)) plot(rj_murder_total ,main = "Total de Homicidios no RJ de 1991-2011" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5 ) plot(rj_murder_capital ,main = "Total de Homicidios no RJ de 1991-2011 (Capital)" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5 ) plot(rj_murder_baixada ,main = "Total de Homicidios no RJ de 1991-2011 (Baixada)" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5 ) plot(rj_murder_interior ,main = "Total de Homicidios no RJ de 1991-2011 (Interior)" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5 )
Agora que vimos que as nossas séries estão corretas, vamos criar alguns objetos usando a função auto.arima().
# Fit ARIMA models fit_rj_murder_total <- auto.arima(rj_murder_total) fit_rj_murder_capital <- auto.arima(rj_murder_capital) fit_rj_murder_baixada <- auto.arima(rj_murder_baixada) fit_rj_murder_interior <- auto.arima(rj_murder_interior)
Essa função retorna o melhor modelo ARIMA de acordo com os valores de AIC, AICc e BIC.
Mas deixando um pouco a teoria de lado, vamos ver como ficaram os ajustes dos modelos ARIMA de acordo com cada uma das séries.
#Plot ARIMA Models par(mfrow=c(2,2)) plot(forecast(fit_rj_murder_total,h=12) ,main = "Total de Homicidios no RJ de 1991-2011 \n Previsão usando ARIMA" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5) plot(forecast(fit_rj_murder_capital,h=12) ,main = "Total de Homicidios no RJ de 91-2011 (Capital) \n Previsão usando ARIMA" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5) plot(forecast(fit_rj_murder_baixada,h=12) ,main = "Total de Homicidios no RJ de 91-2011 (Baixada) \n Previsão usando ARIMA" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5) plot(forecast(fit_rj_murder_interior,h=12) ,main = "Total de Homicidios no RJ de 91-2011 (Interior) \n Previsão usando ARIMA" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5)
Um bom resultado, mas agora vamos ver o que conseguimos fazer usando o package forecast default, com o horizonte de previsão de 12 meses:
#Forecasting using confidence intervals in 80%, 95% and 99% with 12 months ahead ahead_rj_murder_total <- forecast(rj_murder_total, level=c(80, 95, 99), h=12) ahead_rj_murder_capital <- forecast(rj_murder_capital, level=c(80, 95, 99), h=12) ahead_rj_murder_interior <- forecast(rj_murder_interior, level=c(80, 95, 99), h=12) ahead_rj_murder_baixada <- forecast(rj_murder_baixada, level=c(80, 95, 99), h=12)
Modelos ajustados, vamos plotar os resultados:
#Plot forecasting par(mfrow=c(2,2)) plot(ahead_rj_murder_total ,main = "Total de Homicidios no RJ de 91-2011 (Total) \n Previsão usando Forecast package" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5 ,shadecols="oldstyle") plot(ahead_rj_murder_capital ,main = "Total de Homicidios no RJ de 91-2011 (Capital) \n Previsão usando Forecast package" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5 ,shadecols="oldstyle") plot(ahead_rj_murder_baixada ,main = "Total de Homicidios no RJ de 91-2011 (Baixada) \n Previsão usando Forecast package" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5 ,shadecols="oldstyle") plot(ahead_rj_murder_interior ,main = "Total de Homicidios no RJ de 91-2011 (Interior) \n Previsão usando Forecast package" ,xlab = "Ano" ,ylab = "Qtde Homicidios" ,type = "l" ,lwd = 2.5 ,shadecols="oldstyle")
Podemos ver que tanto usando o modelo ARIMA quanto com o uso do package forecast vimos que haverá uma queda no número de homicídios. Com isso podem ser elaboradas (usando mais dados) políticas públicas de segurança.
Refs:
[1] – http://www.inside-r.org/packages/cran/forecast/docs/auto.arima
[2] – http://www.ucamcesec.com.br/wordpress/wp-content/uploads/2011/04/RJ_2011_evolucao.xls
[3] – http://www.statmethods.net/advstats/timeseries.html
Escolha de Variáveis
Um dos maiores desafios na construção de modelos de classificação e de modelos preditivos é saber escolher as variáveis corretas para inclusão no modelo.
Como foi falado aqui em inúmeras vezes, antes de cometer o erro clássico de escalar hardware e software antes de analisar os dados, mesmo que superficialmente; veja se as variáveis do modelo estão adequadas.
O Dan Steinberg neste vídeo fala da importância de se saber importância das variáveis no modelo, no qual não somente essas variáveis vão dar o aspecto de compreensão relativa a qual o espectro de dados são pertinentes para as tarefas do algoritmo; como também, entender essa importância pode dar subsídio para outras análises que por ventura venham a eliminar a fragilidade do modelo.
Um dos aspectos levantados foi que antigamente para levantamento da importância dessas variáveis eram usados técnicas de regressão as quais de acordo os respectivos coeficientes regressores eram utilizados como maneira de ranquear as variáveis.
Contudo, com as técnicas mais modernas de análise de dados, e em especial com novos algoritmos, e a necessidade de modelos que além de terem um alto poder de classificação e predição devem ser compreensíveis, conhecer a importância de cada uma das variáveis ajuda entender o grau de especificidade do modelo.
Em outras palavras, essa atividade auxilia no entendimento do papel – ou força – de cada uma das variáveis no modelo.
Uma heurística interessante que foi explicada no vídeo é conhecida como Leave-One-Variable-Out (LOVO).
A técnica de LOVO consiste em retirar sistematicamente uma variável por vez do modelo, e após isso o modelo preditivo é gerado sem essa variável e de acordo com a variância, isto é, a degradação dos resultados, esse processo auxilia em medir o quanto o modelo perde se aquela variável sair.
Essa heurística é extremamente válida em casos em que se trabalha com heurísticas como Redes Neurais Artificiais, no qual muitas vezes mesmo com alterações em parâmetros de arquitetura (Hidden Layers, Neurônios de Entrada, Neurônios de Saída, Momentum, Taxa de Aprendizado, etc) não há uma visão tão nítida da influência da variável na convergência do modelo.
Medindo a Acurácia das Previsões
Neste paper do Rob Hyndman é apresentado algumas formas de medir a acurácia de modelos de predição. Obrigatório para quem trabalha/estuda modelos de classificação e regressão.
Porque o Overfitting é mais perigoso do que uma acurácia baixa?
Após a leitura desse artigo, fica mais evidente que modelos de dados devem ser testados, se possíveis, com amostras separadas dos conjuntos de dados de treinamento e teste (Holdout).
Data Mining / Machine Learning / Data Analysis 