text:
Sunday, October 21, 2018

Trying out R's keras API - Can we predict the ETH price with recurrent neural networks?

rnn.utf8.md
Hi folks!
In this post, I will combine two of my passions: data science and cryptocurrencies. To be more specific, we will try to predict the minimum and maximum price of Ethereum (ETH) traded on binance for Bitcoins (BTC) within the next 10 minutes. As data, we will use the kline/candlestick data of the previous 180 minutes. Analyzing candlestick data and using it for predictions is a science in its own right, but I do not plan to make a deep dive into that. Instead, we will try to make our predictions with recurrent neural networks (RNN). Like all neural networks, RNNs learn the features on their own and reduce the need for manual feature engineering.
We will try out different RNN architectures and do some hyperparameter training. We will do so using R - I know, some time ago I stated that I prefer python for deep learning. The reason why we still use R instead of python is that R finally got a keras API and I am more than thrilled to take it for a spin.
Before we start I like to mention that I have not seen too many working crypto price prediction models out there - maybe it cannot be done, perhaps it can, and people keep them for themselves to get rich. Anyway, let’s try our luck!

Preparations

General

library(tidyverse)
library(magrittr)
library(keras)
library(mlrMBO)
Next, we set a seed. I do this to make everything reproducible, but honestly, I doubt that it will work. The seed controls all “randomness” within R, but hardly anything will be computed in R but only piped through the API to keras which runs in python. Hence, I am not sure whether the seed has any effect at all (if you follow along on your own machine, then please let me know whether your results differ).
set.seed(1234)

Getting the data

I harvested some data a time ago that you can download from here. We will use it for model training. Later in this post, we will use binance’s API to get more current data and test our models.
After you have downloaded the data we load it into R. It comprises candlestick data of every minute of the trading pairs BTC-ETH and BTC-USDT. I have included BTC-USDT because BTC is still the leading cryptocurrency and its price in $(or USDT) for sure has an impact on the cost of other currencies like ETH. dat_eth <- read_csv("binance_eth.csv") %>% select(-ignore) dat_usdt <- read_csv("binance_btcusdt.csv") %>% select(-ignore) Let’s take a quick glimpse at the data to see with what we are going to work. dat_eth %>% glimpse() ## Observations: 324,977 ## Variables: 11 ##$ timeopen          <dbl> 1.509491e+12, 1.509491e+12, 1.509491e+12, 1....
## $open <dbl> 0.047613, 0.047819, 0.047700, 0.047703, 0.04... ##$ high              <dbl> 0.047729, 0.047819, 0.047700, 0.047703, 0.04...
## $low <dbl> 0.047613, 0.047819, 0.047700, 0.047636, 0.04... ##$ close             <dbl> 0.047686, 0.047819, 0.047700, 0.047684, 0.04...
## $volume <dbl> 22.662, 2.722, 1.859, 0.246, 10.838, 12.593,... ##$ timeclose         <dbl> 1.509491e+12, 1.509491e+12, 1.509491e+12, 1....
## $quotedVolume <dbl> 1.08098804, 0.13016330, 0.08867430, 0.011725... ##$ trades            <int> 22, 7, 1, 3, 6, 11, 37, 18, 38, 6, 14, 3, 19...
## $takerBuyVol <dbl> 15.052, 0.000, 1.859, 0.136, 8.826, 1.592, 2... ##$ takerBuyQuotedVol <dbl> 0.71822354, 0.00000000, 0.08867430, 0.006485...
We merge the data sets with an inner join on their timestamps (timeopen - we remove timeclose because it is just 60 seconds later and therefore redundant).
dat_raw <- inner_join(dat_eth, dat_usdt, by=c("timeopen")) %>%
select(-contains("timeclose"))

Data preparations

I know I have just told you that neural networks reduce the need for manual feature engineering, but a little is still necessary or at least sensible.
One problem is that in some cases the number of trades within the minute is missing. We replace those missings with 0s, and the network will learn how to deal with them. A further problem is that some variables (volumes, number of trades, and rates of BTC-USDT) have quite a positive skew. We normalize them a little by taking their log - this emphasizes differences between lower values while reducing the differences between high values. This step is not necessary for deep learning, but in my limited experience it pays off. One problem with this transformation, however, is that there is no log of values smaller or equal to 0. Hence, we need to replace all zeros with small positive values.
dat_log <- dat_raw %>%
mutate_at(vars(contains("trades")), function(x) ifelse(is.na(x), 0, x)) %>%
function(x) ifelse(x < 0.000001, (unique(sort(x))[2])/5, x) ) %>%
mutate_at(vars(contains("ume"), contains("Vol"), contains("trades"), ends_with("y")), log)
Finally, we remove the timestamp, because we want the model to learn how the momentary trading patterns effect the price in the near future and not just to remember when it has seen which price (if our model relied mainly on its “price memory”, then it would not be much use for predictions with data it has not seen). An alternative to removing the timestamp completely would be to convert it, e.g., in day of week and time of day. Maybe there are some regularities which our model can pick up (e.g., a big trader or trading bot that always buys/sells at the same time), but this might easily change in the future. Hence, relying only on the trading patterns contained in the candle stick data seems more robust to me.
data <- data.matrix(dat_log[,-1])

Train and validation data

Notice: Normally, we would split the data into three sets: training, validation, and test set. However, for now, we will only divide into training and validation set. For testing purposes, we will use data created in the time needed to train our models, which we will retrieve when our models are ready.
Let’s take a look at the overall dimension of the available data.
dim(data)
## [1] 301868     18
So we have a little over 300k minutes and 18 variables.
Let’s use the first 200 days (200*24*60 = 288000 minutes) for training and the rest for validation. Further, we perform one last data preparation step and standardize all variables to have mean = 0 and standard deviation = 1. Please notice that we use the mean and standard deviation of the training data.
Note: Each variable is standardized on its own. In most cases, this is the way to go, but in our case, it is a tiny error. Later we will predict the high AND low values, but only one loss value can be evaluated. Hence, we will pool the error on high and low predictions. The problem is that high and low have similar, but not equal standard deviations. Since both values are scaled to have a standard deviation of one, the error in one weights a little bit more than in the other. A more balanced solution, e.g., scaling both with their pooled SD, would be preferable, but this small error should not impact our results too much.
train_data <- data[1:288000,]
mean <- apply(train_data, 2, mean)
std <- apply(train_data, 2, sd)
data <- scale(data, center = mean, scale = std)
Even if 288k is not that much data (and would easily fit into the memory as a whole), we choose an approach applicable to more substantial amounts of data and slice it into nice little batches which will be processed one after another.
In keras, a generator achieves that. We define a general version that provides us with a batch of inputs and the associated outputs every time it is called.
generator <- function(data, lookback, delay, min_index, max_index,
shuffle = FALSE, batch_size = 640, steps=1){

if (is.null(max_index))
max_index <- nrow(data) - delay - 1
i <- min_index + lookback
function() {
if (shuffle) {
rows <- sample(c((min_index+lookback):max_index), size = batch_size)
} else {
if (i + batch_size >= max_index)
i <<- min_index + lookback
rows <- c(i:min(i+batch_size, max_index))
i <<- i + length(rows)
}

samples <- array(0, dim = c(length(rows),
lookback / steps,
dim(data)[[-1]]))
targets <- array(0, dim = c(length(rows),2))

for (j in 1:length(rows)) {
indices <- seq(rows[[j]] - lookback, rows[[j]],
length.out = dim(samples)[[2]])
samples[j,,] <- data[indices,]

##targets
##high
targets[j,1] <- max(data[(rows[[j]] + 1):(rows[[j]] + delay),2])
##low
targets[j,2] <- min(data[(rows[[j]] + 1):(rows[[j]] + delay),3])
}

list(samples, targets)
}
}
Next, we define a generator that can be used for new data and only includes inputs but no outputs.
generator_predict <- function(data, lookback, delay, min_index, max_index,
shuffle = FALSE, batch_size = 640, steps=1){
if (is.null(max_index))
max_index <- nrow(data) - delay - 1
i <- min_index + lookback
function() {
if (shuffle) {
rows <- sample(c((min_index+lookback):max_index), size = batch_size)
} else {
if (i + batch_size >= max_index)
i <<- min_index + lookback
rows <- c(i:min(i+batch_size, max_index))
i <<- i + length(rows)
}

samples <- array(0, dim = c(length(rows),
lookback / steps,
dim(data)[[-1]]))

for (j in 1:length(rows)) {
indices <- seq(rows[[j]] - lookback, rows[[j]],
length.out = dim(samples)[[2]])
samples[j,,] <- data[indices,]

}

list(samples)
}
}
At last, we define the specific generators that we will use. We use lookback = 180 setting the considered time to three hours, steps = 1 indicating that two adjacent three hours time series are only shifted by one minute, delay = 10 trying to predict 10 minutes ahead, and batch_size = 640 fixing the batch size to 640 three hours time series. train_gen will be used for training, val_gen for validating, and val_gen_pred is needed for predictions with the validation data.
lookback = 180
steps = 1
delay = 10
batch_size <- 640

train_gen <- generator(
data,
lookback = lookback,
delay = delay,
min_index = 1,
max_index = lookback*1600,
shuffle = TRUE,
steps = steps,
batch_size = batch_size
)
val_gen = generator(
data,
lookback = lookback,
delay = delay,
min_index = lookback*1600+1,
max_index = NULL,
steps = steps,
batch_size = batch_size
)

val_gen_pred = generator_predict(
data,
lookback = lookback,
delay = delay,
min_index = lookback*1600+1,
max_index = NULL,
steps = steps,
batch_size = batch_size
)

val_steps <- (nrow(data) - lookback*1600+1 - lookback) / batch_size

Benchmark

Before we can start building the models, we need one last thing: a benchmark against which we compare the predictions of the models.
We will use the mean absolute error resulting from just using the values of the last record as our predictions. This seemingly easy and naive approach is often quite hard to beat. On a practical level, we iterate through the validation set with val_gen, calculate all resulting absolute errors, and average them.
evaluate_naive_method <- function() {
batch_maes <- c()
for (step in 1:val_steps) {
c(samples, targets) %<-% val_gen()
preds <- samples[,dim(samples)[[2]],2:3]
mae <- mean(abs(preds - targets))
batch_maes <- c(batch_maes, mae)
}
mean(batch_maes)
}
(benchmark <- evaluate_naive_method())
## [1] 0.004545522
So our models need to have a MAE lower than 0.00455 to be any good at all.

Training

Now, all preparations are done, and the exciting part starts. We train models and try to optimize their performance.

Single layer RNN

We start with the most basic RNN architecture - a single layer net. But even with a simple RNN, there are some options from which we need to choose. E.g.:
• What is the optimal size?
• Which cell type should we use (LSTM or GRU)?
• Which optimization algorithm performs best?
Experience helps to find proper starting configurations, but you never really know what configuration works best in a given situation until you try. There are several approaches to find the optimal hyperparameters. The most exhaustive method is a grid search, which iterates through all parameter combinations and chooses the best in the end. The downside is that it can take quite a long time to iterate through all combinations (especially if several hyperparameters need optimization). Since we do not have too much data and the models will not be too complex, grid search would probably be viable in a bearable amount of time. Still, we choose a more efficient approach (with the downside that we probably do not end up with the optimal but with a good solution). We will use bayesian optimization provided by the mlrMBO package.
To apply mlrMBO’s optimization (or any other optimization provided by other packages) we need to wrap the model building into a function that takes the hyperparameters we like to tune as arguments and returns a single numeric value that we want to minimize or maximize.
As arguments, we choose the number of units, their type, and the optimization algorithm and return the mean absolute error in the validation sample.
Note: While the hyperparameter training tries to optimize the MAE, the model training per se is set to optimize the mean squared error. This is somewhat inconsistent, but to me both types of error are important, and I like the compromise of using both. If you disagree, then please let me know in the discussion section.
Note: The different cell types and optimization algorithms have further parameters on their own and including them into the hyperparameter training could be sensible, but we stick to the defaults for now because each additional parameter multiplies hyperparameter training time.
Note: For more information on model building and training check out the later section of this post where we will train models with “optimal” hyperparameters.
simpleRNN_fit <- function(unit_count, cell_type, opt_type){

uc <- unit_count*8

## define callbacks
callbacks_list <- list(
callback_early_stopping(
monitor = "val_loss",
patience = 5
),
callback_model_checkpoint(
filepath = "01_simpleRNN.h5",
monitor = "val_loss",
save_best_only = TRUE
)
)
## build model
model_simpleRNN <- keras_model_sequential()

if(cell_type == "gru"){
model_simpleRNN %<>%
layer_gru(units = uc,
input_shape = list(NULL, dim(data)[[-1]]))
} else if(cell_type == "lstm"){
model_simpleRNN %<>%
layer_lstm(units = uc,
input_shape = list(NULL, dim(data)[[-1]]))
}

model_simpleRNN %<>% layer_dense(units = 2)

## compile model
keras::compile(
model_simpleRNN,
loss = "mse",
metrics = c("mse", "mae")
)
} else if (opt_type == "nadam"){
keras::compile(
model_simpleRNN,
loss = "mse",
metrics = c("mse", "mae")
)
} else if (opt_type == "rmsprop"){
keras::compile(
model_simpleRNN,
optimizer = optimizer_rmsprop(),
loss = "mse",
metrics = c("mse", "mae")
)
keras::compile(
model_simpleRNN,
loss = "mse",
metrics = c("mse", "mae")
)
}

## fit model
model_simpleRNN %>% fit_generator(
train_gen,
callbacks = callbacks_list,
steps_per_epoch = 450,
epochs = 2000,
validation_data = val_gen,
validation_steps = val_steps
)

## evaluate model
eval_results <- evaluate_generator(model_simpleRNN, val_gen, val_steps)

evaluate_generator(model_simpleRNN,val_gen, val_steps)$mean_absolute_error %>% as.numeric() } Next, we define the hyperparameter space of possible values … simpleRNN_ps = makeParamSet( makeIntegerParam("unit_count", lower = 1, upper = 12), makeDiscreteParam("cell_type", values = c("gru", "lstm")), makeDiscreteParam("opt_type", values = c("adam", "nadam", "rmsprop", "adadelta")) ) … and randomly choose some points within it as initial data for the optimization. For the blog post, we take only 14 points, but a higher number might be sensible. Further, please be aware that you can manually add combinations as well. So, if you believe to have a good guess, then include it. simpleRNN_des = generateDesign(n = 14, par.set = simpleRNN_ps) Next, we evaluate the chosen configurations. simpleRNN_des %<>% mutate( y = pmap(list(unit_count, cell_type, opt_type), simpleRNN_fit) ) ## Error in mutate_impl(.data, dots): Evaluation error: AttributeError: 'Sequential' object has no attribute 'metrics'. OK, that is weird. But if the tidy approach is not working, then let’s do it the old-fashioned way (if you know what I have done wrong - please let me know). for (i in 1:dim(simpleRNN_des)[1]){ simpleRNN_des$y[i] <- simpleRNN_fit(simpleRNN_des$unit_count[i], simpleRNN_des$cell_type[i],
simpleRNN_des$opt_type[i]) } Let’s look at the resulting starting values. simpleRNN_des ## unit_count cell_type opt_type y ## 1 6 lstm adam 0.006573573 ## 2 1 lstm nadam 0.004850281 ## 3 7 gru rmsprop 0.006180876 ## 4 11 lstm rmsprop 0.007254474 ## 5 9 gru adam 0.005405560 ## 6 10 gru nadam 0.004455248 ## 7 2 gru adam 0.005086179 ## 8 2 gru nadam 0.005570578 ## 9 3 lstm adadelta 0.006913974 ## 10 5 gru adadelta 0.006262964 ## 11 8 lstm rmsprop 0.005818588 ## 12 10 lstm adadelta 0.004891048 ## 13 4 gru adam 0.004516485 ## 14 12 lstm adadelta 0.005799765 Next, we use the starting values to initialize an OptState object, which links the problem to the result and keeps track of the optimization progress. simpleRNN_ctrl = makeMBOControl() simpleRNN_ctrl = setMBOControlInfill(simpleRNN_ctrl, crit = crit.ei) simpleRNN_opt.state = initSMBO( par.set = simpleRNN_ps, design = simpleRNN_des, control = simpleRNN_ctrl, minimize = TRUE, noisy = TRUE) We update the OptState 16 times (a completely arbitrary number and probably too small if it was not just for a blog post). for (i in 1:16){ prop = proposePoints(simpleRNN_opt.state) y = simpleRNN_fit(prop$prop.points$unit_count, prop$prop.points$cell_type, prop$prop.points$opt_type) updateSMBO(simpleRNN_opt.state, x = prop$prop.points,
y = y)
}
Let’s have a look at which combination would be evaluated next.
proposePoints(simpleRNN_opt.state)
## $prop.points ## unit_count cell_type opt_type ## 34 11 gru nadam ## ##$propose.time
## [1] 0.219
##
## $prop.type ## [1] "infill_ei" ## ##$crit.vals
##               [,1]
## [1,] -2.135247e-07
##
## $crit.components ## se mean ## 1 0.0003345841 0.005409518 ## ##$errors.model
## [1] NA
##
## attr(,"class")
## [1] "Proposal" "list"
We could easily go on and evaluate this combination or others and add them to the OptState, but for now, we conclude the hyperparameter training.
finalizeSMBO(simpleRNN_opt.state)
## Recommended parameters:
## Objective: y = 0.004
##
## Optimization path
## 14 + 16 entries in total, displaying last 10 (or less):
##    unit_count cell_type opt_type           y dob eol error.message
## 21          1      lstm    nadam 0.004612412  21  NA          <NA>
## 22          1      lstm    nadam 0.004866440  22  NA          <NA>
## 24          1      lstm    nadam 0.004694078  24  NA          <NA>
## 25          6      lstm     adam 0.004983207  25  NA          <NA>
## 26          6      lstm     adam 0.004760775  26  NA          <NA>
## 28          6      lstm     adam 0.005319725  28  NA          <NA>
## 29         10       gru    nadam 0.004878603  29  NA          <NA>
## 30          6      lstm     adam 0.005167965  30  NA          <NA>
##    exec.time            ei error.model train.time prop.type propose.time
## 21        NA -1.045612e-06        <NA>          0    manual           NA
## 22        NA -3.552289e-06        <NA>          0    manual           NA
## 23        NA -2.116553e-06        <NA>          0    manual           NA
## 24        NA -7.931060e-07        <NA>          0    manual           NA
## 25        NA -1.027754e-06        <NA>          0    manual           NA
## 26        NA -2.580085e-06        <NA>          0    manual           NA
## 27        NA -6.329430e-07        <NA>          0    manual           NA
## 28        NA -2.318679e-06        <NA>          0    manual           NA
## 29        NA -2.471307e-07        <NA>          0    manual           NA
## 30        NA -2.435425e-07        <NA>          0    manual           NA
##              se        mean
## 21 0.0003982259 0.005415668
## 22 0.0003481275 0.005127366
## 23 0.0006486102 0.005972417
## 24 0.0002483238 0.005037904
## 25 0.0005510765 0.005845696
## 26 0.0005269609 0.005614215
## 27 0.0004980626 0.005772202
## 28 0.0004598631 0.005461888
## 29 0.0001859679 0.004944392
## 30 0.0003434029 0.005423978
The resulting recommended configuration is …
simpleRNN_opt.state$opt.result$mbo.result$x ##$unit_count
## [1] 10
##
## $cell_type ## [1] "gru" ## ##$opt_type
## [1] "nadam"
… resulting in an estimated MAE of
simpleRNN_opt.state$opt.result$mbo.result$y ## [1] 0.004455248 Hm, that is slightly better than the naive benchmark of 0.00455. Hopefully, more complex architectures fare better. Stacked RNN One contender is the stacked RNN architecture. Instead of using just one layer of cells, it stacks several layers over one another. Again, we need to wrap the model building into a function. As arguments, we use the number of layers, the number of units, and the shape (are further layers getting smaller, greater, or stay the same as the previous layers). We will not re-evaluate cell type and optimization algorithm and use the optimal selection in the single layer model. Doing so is a pragmatic decision, but not necessarily the best one. Since all parameters interact with one another, other settings might be preferable for this architecture. So if you have time - check it out. Note: All but the top layer are set to return_sequences = TRUE. Otherwise, they would return only the state of the last cell which would not give consecutive layers much information to train on. stackedRNN_fit <- function(layer_count, unit_count, shape){ ## check input parameter Check <- ArgumentCheck::newArgCheck() if (layer_count < 2) ArgumentCheck::addError( msg = "'layer_count' must be >= 2", argcheck = Check ) if (unit_count < 1) ArgumentCheck::addError( msg = "'unit_count' must be >= 2", argcheck = Check ) if (shape < -2 | shape > 2) ArgumentCheck::addError( msg = "'shape' must be within [-2;2]", argcheck = Check ) ArgumentCheck::finishArgCheck(Check) ## set-up unit_base <- 8 uc <- unit_count*unit_base ## define callbacks callbacks_list <- list( callback_early_stopping( monitor = "val_loss", patience = 5 ), callback_model_checkpoint( filepath = "02_stackedRNN.h5", monitor = "val_loss", save_best_only = TRUE ) ) ## build model model_stackedRNN <- keras_model_sequential() %>% layer_gru(units = uc, input_shape = list(NULL, dim(data)[[-1]]), return_sequences = TRUE) growth <- 1/(shape*2) growth[shape == 0] <- 0 for (i in 2:layer_count){ if (i < layer_count){ model_stackedRNN %<>% layer_gru(units = uc, return_sequences = TRUE) } else { model_stackedRNN %<>% layer_gru(units = uc) } uc <- uc + uc * growth uc <- c(uc, unit_base) %>% max() } model_stackedRNN %<>% layer_dense(units = 2) ## compile model keras::compile( model_stackedRNN, optimizer = optimizer_nadam(), loss = "mse", metrics = c("mse", "mae") ) ## fit model model_stackedRNN %>% fit_generator( train_gen, callbacks = callbacks_list, steps_per_epoch = 450, epochs = 2000, validation_data = val_gen, validation_steps = val_steps ) ## load model model_stackedRNN <- load_model_hdf5("02_stackedRNN.h5") eval_results <- evaluate_generator(model_stackedRNN, val_gen, val_steps) evaluate_generator( model_stackedRNN,val_gen, val_steps)$mean_absolute_error %>%
as.numeric()
}
Again, we define the parameter space, …
stackedRNN_ps = makeParamSet(
makeIntegerParam("layer_count", lower = 2, upper = 5),
makeIntegerParam("unit_count", lower = 1, upper = 10),
makeIntegerParam("shape", lower = -2, upper = 2)
)
… choose a number of starting combinations, …
stackedRNN_des = generateDesign(n = 15, par.set = stackedRNN_ps)
… and evaluate them.
for (i in 1:dim(stackedRNN_des)[1]){
stackedRNN_des$y[i] <- stackedRNN_fit(stackedRNN_des$layer_count[i],
stackedRNN_des$unit_count[i], stackedRNN_des$shape[i])
}
Let’s look at the starting values.
stackedRNN_des
##    layer_count unit_count shape           y
## 1            3          7     2 0.005768042
## 2            2          3     2 0.006272834
## 3            2          2     1 0.004910406
## 4            4          8     0 0.006819058
## 5            2          1     0 0.006527344
## 6            4          7    -1 0.006985018
## 7            4          9     1 0.006050816
## 8            2          6     1 0.005260754
## 9            3          5    -2 0.005264357
## 10           5         10    -1 0.005853208
## 11           5          1    -2 0.005304723
## 12           5         10     0 0.006696024
## 13           3          3    -1 0.005411024
## 14           3          4    -2 0.005591315
## 15           5          5     2 0.005284705
Next, we create the OptState
stackedRNN_ctrl = makeMBOControl()
stackedRNN_ctrl = setMBOControlInfill(stackedRNN_ctrl, crit = crit.ei)
stackedRNN_opt.state = initSMBO(
par.set = stackedRNN_ps,
design = stackedRNN_des,
control = stackedRNN_ctrl,
minimize = TRUE,
noisy = TRUE)
… and let it train for a while.
for (i in 1:15){
prop = proposePoints(stackedRNN_opt.state)
y = stackedRNN_fit(prop$prop.points$layer_count,
prop$prop.points$unit_count,
prop$prop.points$shape)
x = prop$prop.points, y = y) } The next configuration to be evaluated would be … proposePoints(stackedRNN_opt.state) ##$prop.points
##   layer_count unit_count shape
## 9           4          4     1
##
## $propose.time ## [1] 0.242 ## ##$prop.type
## [1] "infill_ei"
##
## $crit.vals ## [,1] ## [1,] -0.002232555 ## ##$crit.components
##            se        mean
## 1 0.002886951 0.002986154
##
## $errors.model ## [1] NA ## ## attr(,"class") ## [1] "Proposal" "list" … but we finish hyperparameter training for now. finalizeSMBO(stackedRNN_opt.state) ## Recommended parameters: ## layer_count=2; unit_count=1; shape=1 ## Objective: y = 0.005 ## ## Optimization path ## 15 + 15 entries in total, displaying last 10 (or less): ## layer_count unit_count shape y dob eol error.message ## 21 2 2 1 0.005336026 21 NA <NA> ## 22 2 1 1 0.004736808 22 NA <NA> ## 23 2 1 1 0.005203352 23 NA <NA> ## 24 2 1 1 0.005320675 24 NA <NA> ## 25 2 1 1 0.005275674 25 NA <NA> ## 26 3 3 0 0.005190356 26 NA <NA> ## 27 3 10 0 0.005690224 27 NA <NA> ## 28 5 10 2 0.022936167 28 NA <NA> ## 29 5 5 0 0.005759698 29 NA <NA> ## 30 2 6 2 0.005913317 30 NA <NA> ## exec.time ei error.model train.time prop.type propose.time ## 21 NA -1.005071e-05 <NA> 0 manual NA ## 22 NA -9.477561e-06 <NA> 0 manual NA ## 23 NA -1.081139e-05 <NA> 0 manual NA ## 24 NA -8.118479e-06 <NA> 0 manual NA ## 25 NA -6.718298e-06 <NA> 0 manual NA ## 26 NA -5.815763e-06 <NA> 0 manual NA ## 27 NA -1.293808e-05 <NA> 0 manual NA ## 28 NA -6.261678e-06 <NA> 0 manual NA ## 29 NA -1.479915e-03 <NA> 0 manual NA ## 30 NA -1.072348e-03 <NA> 0 manual NA ## se mean ## 21 0.0001833585 0.005012616 ## 22 0.0002095543 0.005063587 ## 23 0.0002212725 0.005017067 ## 24 0.0002191364 0.005042574 ## 25 0.0002273631 0.005077032 ## 26 0.0006086373 0.005927065 ## 27 0.0004542407 0.005423839 ## 28 0.0004028942 0.005448340 ## 29 0.0037524021 0.004771086 ## 30 0.0018982689 0.004173024 Now, the optimal configuration guess is … stackedRNN_opt.state$opt.result$mbo.result$x
## $layer_count ## [1] 2 ## ##$unit_count
## [1] 1
##
## $shape ## [1] 1 … and is estimated to yield a MAE of stackedRNN_opt.state$opt.result$mbo.result$y
## [1] 0.004736808
That is a bummer. Even after some hyperparameter training stacked RNNs perform worse on the problem than our naive benchmark.

1d-conv RNN

Let’s try one more architecture. This time we start with a 1d-convolutional layer before we move to recurrent layers. Again, we need to wrap model building into a function to be able to optimize the hyperparameters. Specific to the 1d conv layer are the number of filters and the kernel size.
convRNN_fit <- function(filters, kernel_size, rnn_count, rnn_units){

## check input parameter
Check <- ArgumentCheck::newArgCheck()

if (filters < 1)
msg = "'filters' must be >= 1",
argcheck = Check
)

if (kernel_size < 2)
msg = "'kernel_size' must be >= 2",
argcheck = Check
)

if (rnn_count < 1)
msg = "'rnn_count' must be >= 1",
argcheck = Check
)

if (rnn_units < 1)
msg = "'rnn_units' must be >= 1",
argcheck = Check
)

ArgumentCheck::finishArgCheck(Check)

## set-up

unit_base <- 8

uc <- rnn_units*unit_base

fc <- filters*unit_base

## define callbacks
callbacks_list <- list(
callback_early_stopping(
monitor = "val_loss",
patience = 5
),
callback_model_checkpoint(
filepath = "03_1dconvRNN.h5",
monitor = "val_loss",
save_best_only = TRUE
)
)
## build model
model_convRNN <- keras_model_sequential() %>%
layer_conv_1d(filters = fc,
kernel_size = kernel_size,
activation = "relu",
input_shape = list(NULL, dim(data)[[-1]]))

for (i in 1:rnn_count){
if (i < rnn_count){
model_convRNN %<>%
layer_gru(units = uc,
return_sequences = TRUE)
uc <- uc + uc * 0.5
} else {
model_convRNN %<>%
layer_lstm(units = uc)
}

}

model_convRNN %<>% layer_dense(units = 2)

## compile model
keras::compile(
model_convRNN,
loss = "mse",
metrics = c("mse", "mae")
)

## fit model
model_convRNN %>% fit_generator(
train_gen,
callbacks = callbacks_list,
steps_per_epoch = 450,
epochs = 2000,
validation_data = val_gen,
validation_steps = val_steps
)

eval_results <- evaluate_generator(model_convRNN,val_gen, val_steps)

evaluate_generator(model_convRNN,val_gen, val_steps)$mean_absolute_error %>% as.numeric() } Next, we define the hyperparameter space, … convRNN_ps = makeParamSet( makeIntegerParam("filters", lower = 1, upper = 8), #gets multiplied by 8 makeIntegerParam("kernel_size", lower =3, upper = 10), makeIntegerParam("rnn_count", lower = 1, upper = 2), makeIntegerParam("rnn_units", lower = 1, upper = 10) #gets multiplied by 8 ) … choose random starting combinations, … convRNN_des = generateDesign(n = 14, par.set = convRNN_ps) … and evaluate them. for (i in 1:dim(stackedRNN_des)[1]){ convRNN_des$y[i] <- convRNN_fit(convRNN_des$filters[i], convRNN_des$kernel_size[i],
convRNN_des$rnn_count[i], convRNN_des$rnn_units[i])
}
Let’s look at the starting values.
convRNN_des
##    filters kernel_size rnn_count rnn_units           y
## 1        5           6         2         3 0.006306844
## 2        3           6         2         8 0.004858218
## 3        4           9         2        10 0.006397791
## 4        8           9         1         1 0.013997144
## 5        3           4         1         9 0.005559878
## 6        2          10         2         8 0.005239270
## 7        5           3         2         5 0.005935301
## 8        6           7         1         6 0.006564416
## 9        2           4         1         6 0.005154384
## 10       7           4         1         1 0.006090563
## 11       7           9         1         4 0.006378257
## 12       8           7         2         7 0.009060436
## 13       1           5         1         3 0.004755455
## 14       4           8         2         4 0.008622950
Next, we create the OptState
convdRNN_ctrl = makeMBOControl()
convdRNN_ctrl = setMBOControlInfill(convdRNN_ctrl, crit = crit.ei)
convRNN_opt.state = initSMBO(
par.set = convRNN_ps,
design = convRNN_des,
control = convdRNN_ctrl,
minimize = TRUE,
noisy = TRUE)
… and let it train for a while.
for (i in 1:16){
prop = proposePoints(convRNN_opt.state)
y = convRNN_fit(prop$prop.points$filters,
prop$prop.points$kernel_size,
prop$prop.points$rnn_count,
prop$prop.points$rnn_units)
x = prop$prop.points, y = y) } The next configuration to be evaluated would be … proposePoints(convRNN_opt.state) ##$prop.points
##    filters kernel_size rnn_count rnn_units
## 37       1           4         1         5
##
## $propose.time ## [1] 0.283 ## ##$prop.type
## [1] "infill_ei"
##
## $crit.vals ## [,1] ## [1,] -0.0004018811 ## ##$crit.components
##            se        mean
## 1 0.001935427 0.005617492
##
## $errors.model ## [1] NA ## ## attr(,"class") ## [1] "Proposal" "list" … but we finish hyperparameter training for now. finalizeSMBO(convRNN_opt.state) ## Recommended parameters: ## filters=1; kernel_size=4; rnn_count=1; rnn_units=1 ## Objective: y = 0.005 ## ## Optimization path ## 14 + 16 entries in total, displaying last 10 (or less): ## filters kernel_size rnn_count rnn_units y dob eol ## 21 3 10 2 10 0.014055169 21 NA ## 22 4 5 2 2 0.006309814 22 NA ## 23 2 6 2 4 0.005546650 23 NA ## 24 3 5 2 5 0.006806981 24 NA ## 25 1 5 2 5 0.005864668 25 NA ## 26 1 4 1 1 0.004710464 26 NA ## 27 2 6 1 8 0.005526419 27 NA ## 28 1 3 1 3 0.005252020 28 NA ## 29 1 7 1 1 0.005396571 29 NA ## 30 2 8 2 6 0.005885977 30 NA ## error.message exec.time ei error.model train.time prop.type ## 21 <NA> NA -0.0002779504 <NA> 0 manual ## 22 <NA> NA -0.0002977433 <NA> 0 manual ## 23 <NA> NA -0.0002868495 <NA> 0 manual ## 24 <NA> NA -0.0006456876 <NA> 0 manual ## 25 <NA> NA -0.0005032000 <NA> 0 manual ## 26 <NA> NA -0.0005196650 <NA> 0 manual ## 27 <NA> NA -0.0004569890 <NA> 0 manual ## 28 <NA> NA -0.0004400800 <NA> 0 manual ## 29 <NA> NA -0.0004613463 <NA> 0 manual ## 30 <NA> NA -0.0003792393 <NA> 0 manual ## propose.time se mean ## 21 NA 0.001043168 0.005069259 ## 22 NA 0.002566205 0.006859780 ## 23 NA 0.002509819 0.006834781 ## 24 NA 0.001832746 0.004933277 ## 25 NA 0.001711494 0.005150923 ## 26 NA 0.001616271 0.005023404 ## 27 NA 0.001674741 0.005186315 ## 28 NA 0.001154470 0.004752034 ## 29 NA 0.001883278 0.005386079 ## 30 NA 0.001512320 0.005228841 Now, the optimal configuration guess is … convRNN_opt.state$opt.result$mbo.result$x
## $filters ## [1] 1 ## ##$kernel_size
## [1] 4
##
## $rnn_count ## [1] 1 ## ##$rnn_units
## [1] 1
… and is estimated to yield a MAE of
convRNN_opt.state$opt.result$mbo.result$y ## [1] 0.004710464 Hm … that is slightly better than the stacked RNN but still worse than our naive benchmark. Learn “optimal” models After hyperparameter optimization I am not overly confident that we will end up with great or even good models. Anyhow, let’s use the optimal hyperparameter configurations to train our models. We start by clearing keras’ cache. k_clear_session() single layer RNN Again, we start with the single layer RNN. But before we build the models we set two callbacks. The first indicates that we want to abort training when the validation loss does not decrease for 10 consecutive training epochs. The second saves the model every time a new minimum validation loss is found. callbacks_simple <- list( callback_early_stopping( monitor = "val_loss", patience = 10 ), callback_model_checkpoint( filepath = "01_simple.h5", monitor = "val_loss", save_best_only = TRUE ) ) Next, we define the architecture of the RNN. model_simple <- keras_model_sequential() %>% layer_gru(units = 80, input_shape = list(NULL, dim(data)[[-1]])) %>% layer_dense(units = 2)  summary() gives us an overview of the model. summary(model_simple) ## ___________________________________________________________________________ ## Layer (type) Output Shape Param # ## =========================================================================== ## gru_1 (GRU) (None, 80) 23760 ## ___________________________________________________________________________ ## dense_1 (Dense) (None, 2) 162 ## =========================================================================== ## Total params: 23,922 ## Trainable params: 23,922 ## Non-trainable params: 0 ## ___________________________________________________________________________ Then we need to compile the model, defining which optimization algorithm and loss to use and which metrics should be tracked. The compile model can then be trained on our data. model_simple %>% compile( optimizer = optimizer_nadam(), loss = "mse", metrics = c("mse", "mae") ) history_simple <- model_simple %>% fit_generator( train_gen, callbacks = callbacks_simple, steps_per_epoch = 450, epochs = 100, validation_data = val_gen, validation_steps = val_steps ) After training has finished we can visualize how the metrics have changed during training. plot(history_simple) The plot is a little bit hard to decipher. plot.keras_training_history() uses ggplot2 as the graphical backbone, so we can use ggplot commands to amend the plot. Let’s try it out by scaling the y-axis to log10. plot(history_simple) + scale_y_log10() This is a little bit more informative. From my point of view the training process seems quite normal. Let’s check how well the model performs in the end. model_simple <- load_model_hdf5("01_simple.h5") evaluate_generator(model_simple,val_gen, val_steps) ##$loss
## [1] 5.61021e-05
##
## $mean_squared_error ## [1] 5.61021e-05 ## ##$mean_absolute_error
## [1] 0.004650629
OK, that is worse than our naive benchmark of 0.004546 - nothing to be excited yet.

stacked RNN

Next, we train the “optimal” stacked RNN.
We use the same callbacks as before, but save the best model to another file.
callbacks_stacked <- list(
callback_early_stopping(
monitor = "val_loss",
patience = 10
),
callback_model_checkpoint(
filepath = "02_stacked.h5",
monitor = "val_loss",
save_best_only = TRUE
)
)
Then we build the model …
model_stacked <- keras_model_sequential() %>%
layer_gru(units = 8,
input_shape = list(NULL, dim(data)[[-1]]),
return_sequences = TRUE) %>%
layer_gru(units = 12) %>%
layer_dense(units = 2)
… and check its summary().
summary(model_stacked)
## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #
## ===========================================================================
## gru_2 (GRU)                      (None, None, 8)               648
## ___________________________________________________________________________
## gru_3 (GRU)                      (None, 12)                    756
## ___________________________________________________________________________
## dense_2 (Dense)                  (None, 2)                     26
## ===========================================================================
## Total params: 1,430
## Trainable params: 1,430
## Non-trainable params: 0
## ___________________________________________________________________________
Next, we can compile and train it.
model_stacked %>% compile(
loss = "mse",
metrics = c("mse", "mae")
)

history_stacked <- model_stacked %>% fit_generator(
train_gen,
callbacks = callbacks_stacked,
steps_per_epoch = 450,
epochs = 100,
validation_data = val_gen,
validation_steps = val_steps
)
Let’s have a look at the training process …
plot(history_stacked) +
scale_y_log10()

… and how the model performs in the end.
model_stacked <- load_model_hdf5("02_stacked.h5")

evaluate_generator(model_stacked,val_gen, val_steps)
## $loss ## [1] 5.549375e-05 ## ##$mean_squared_error
## [1] 5.549375e-05
##
## $mean_absolute_error ## [1] 0.004800438 Hm … no reason to be happy, either. conv RNN Let’s train the last model, the 1d convolutional + 2 layer RNN net. Again, we define the callbacks, … callbacks_conv <- list( callback_early_stopping( monitor = "val_loss", patience = 10 ), callback_model_checkpoint( filepath = "03_conv.h5", monitor = "val_loss", save_best_only = TRUE ) ) … build the model, … model_conv <- keras_model_sequential() %>% layer_conv_1d(filters = 8, kernel_size = 4, activation = "relu", input_shape = list(NULL, dim(data)[[-1]])) %>% layer_gru(units = 8) %>% layer_dense(units = 2)  … check its summary(), … summary(model_conv) ## ___________________________________________________________________________ ## Layer (type) Output Shape Param # ## =========================================================================== ## conv1d_1 (Conv1D) (None, None, 8) 584 ## ___________________________________________________________________________ ## gru_4 (GRU) (None, 8) 408 ## ___________________________________________________________________________ ## dense_3 (Dense) (None, 2) 18 ## =========================================================================== ## Total params: 1,010 ## Trainable params: 1,010 ## Non-trainable params: 0 ## ___________________________________________________________________________ … compile, and train it. model_conv %>% compile( optimizer = optimizer_nadam(), loss = "mse", metrics = c("mse", "mae") ) history_conv <- model_conv %>% fit_generator( train_gen, callbacks = callbacks_conv, steps_per_epoch = 450, epochs = 100, validation_data = val_gen, validation_steps = val_steps ) Let’s have a look at the training process … plot(history_conv) + scale_y_log10() … and how the model performs in the end. model_conv <- load_model_hdf5("03_conv.h5") evaluate_generator(model_conv,val_gen, val_steps) ##$loss
## [1] 5.487677e-05
##
## $mean_squared_error ## [1] 5.487677e-05 ## ##$mean_absolute_error
## [1] 0.004629169
Hm … it does not get better.

pool models

At the moment it seems like we have wasted our time so far. But there is still hope. Maybe none of the models is any good on its own, but we can combine them to an ensemble which might perform better.
To do so we need to combine the predictions of all three models. We use predict_generator() to get those …
pred_simple <- predict_generator(model_simple, val_gen_pred, val_steps)
pred_stacked <- predict_generator(model_stacked, val_gen_pred, val_steps)
pred_conv <- predict_generator(model_conv, val_gen_pred, val_steps)
… and merge them into array.
pred_all <- abind::abind(pred_simple, pred_stacked, pred_conv, along = 3)
Let’s check the size of the array as a quick control if everything worked.
dim(pred_all)
## [1] 13461     2     3
This seems right - we have 13461 measurement points, and 3 predictions of each of the 2 targets.

mean

The most common way to pool multiple predictions is to average them.
pred_mean <- apply(pred_all, c(1,2), mean)
After we have done so, we check the size of the array again.
dim(pred_mean)
## [1] 13461     2
The resulting array has only two dimensions, because the three predictions are combined into their average.
Next, we need the targets with which we will compare the predictions.
targets_list <- list()
for (step in 1:val_steps) {
c(samples, targets) %<-% val_gen()
targets_list[[step]] <- targets
}

targets <- do.call(rbind, targets_list)
dim(targets)
## [1] 13461     2
The target array and the prediction array have the same size - so probably everything has worked so far.
Now, we can evaluate the pooled prediction by calculating the mean absolute error.
mae <- mean(abs(pred_mean - targets))

mae
## [1] 0.00445963
Hm … at last the naive benchmark is beaten. However, only by a little. Maybe a weighted average that puts more emphasis on the better performing individual models would work better. If you try it out, please let us know. We move on to other pooling options.

max/min (liberal)

One of those is choosing the minimum or maximum. Since we have two targets we even can combine both approaches. Let’s start with the most liberal prediction, choosing the maximum high and minimum low.
pred_maxmin <- reshape2::melt(pred_all) %>%
group_by(Var1) %>%
summarize(
high = max(1),
low = min(2)
) %>%
ungroup() %>%
select(-Var1) %>%
as.matrix()
Again, we calculate the MAE.
mae <- mean(abs(pred_maxmin - targets))

mae
## [1] 0.005214722
Being too liberal does not seem to be the way to go in this case.

min/max (conservative)

Let’s try the opposite and be as conservative as possible with the minimum high and maximum low.
pred_minmax <- reshape2::melt(pred_all) %>%
group_by(Var1) %>%
summarize(
high = min(1),
low = max(2)
) %>%
ungroup() %>%
select(-Var1) %>%
as.matrix()
One more time, we calculate the MAE.
mae <- mean(abs(pred_minmax - targets))

mae
## [1] 0.004352053
That looks more like it. Let’s check how many percent this prediction is better than our naive benchmark.
((1 - mae/benchmark) * 100) %>%
round(2)
## [1] 4.26
4.26% - I am not overly thrilled, but it demonstrates that some prediction of short term crypto price development is possible. At the end of the post I will discuss some ideas to improve this result even further. But first we need to test it on truly new data.

More data

Preparations

First we need to check the time-frame that is already covered in our data.
(start_old <- dat_raw$timeopen %>% min()) ## [1] 1.509491e+12 (end_old <- dat_raw$timeopen %>% max())
## [1] 1.528989e+12
Cool - for those of us that cannot readily understand UNIX timestamps we convert it into human readable dates.
as.Date(as.POSIXct(start_old/1000, origin="1970-01-01"))
## [1] "2017-10-31"
as.Date(as.POSIXct(end_old/1000, origin="1970-01-01"))
## [1] "2018-06-14"
OK, our data includes the crypto bubble of December 2017 and its burst in early 2018. Accordingly, prices behaved quite differently in our training data than in time span that we will download in a moment (in the last months the trend was sideways to slightly downwards without any big jumps). Further, binance, although the #1 crypto-exchange now, is quite young and has grown considerable in the time-frame covered by the training data. I am intrigued how our models will fare with more current data.
First, we need a helper function to connect to the binance API.
get_chartdata_unix <- function(symbol,
interval = "1m",
start = NULL,
limit = 500){

raw <- httr::GET(
url = "https://api.binance.com",
path = "/api/v1/klines",
query = list (
symbol = symbol,
interval = interval,
limit = limit,
startTime = start
)) %>%
httr::content(as="text", encoding="UTF-8") %>%
jsonlite::fromJSON()

as_tibble(raw) %>%
mutate_all(as.numeric) %>%
setNames(c("timeopen", "open", "high", "low", "close", "volume", "timeclose", "quotedVolume", "trades",
}
Next, we check the time now as an upper margin.
(now <- as.numeric(Sys.time())*1000)
## [1] 1.539112e+12
For those of us (including me) that still cannot parse UNIX timestamps.
as.Date(as.POSIXct(now/1000, origin="1970-01-01"))
## [1] "2018-10-09"
And add one minute (60.000 milliseconds) to the last available point in our data as a lower margin.
start_new <- end_old + 60000
stamp <- start_new
binance_eth_new <- NULL

while(stamp < now){
binance_eth_new <- binance_eth_new %>%
## and slowly grow the data set
bind_rows(get_chartdata_unix("ETHBTC", start = stamp))

stamp <- stamp + 60000*500
## add a little sleep to avoid getting banned
Sys.sleep(0.5)
}

## write to disk
write_csv(binance_eth_new, "binance_eth_new.csv")
… and the new BTC-USDT data.
stamp <- start_new
binance_btcusdt_new <- NULL

while(stamp < now){
binance_btcusdt_new <- binance_btcusdt_new %>%
bind_rows(get_chartdata_unix("BTCUSDT", start = stamp))

stamp <- stamp + 60000*500
Sys.sleep(0.5)
}

write_csv(binance_btcusdt_new, "binance_btcusdt_new.csv")
Combine the new data sets …
dat_raw_new <- inner_join(binance_eth_new,
binance_btcusdt_new,
by=c("timeopen")) %>%
select(-contains("timeclose"))
… and perform all other preparations steps.
dat_new <- dat_raw_new %>%
mutate_at(vars(contains("trades")), function(x) ifelse(is.na(x), 0, x)) %>%
function(x) ifelse(x < 0.000001, (unique(sort(x))[2])/5, x) ) %>%
mutate_at(vars(contains("ume"), contains("Vol"), contains("trades"), ends_with("y")), log) %>%
select(-timeopen, -starts_with("ignore")) %>%
data.matrix() %>%
scale(center = mean, scale = std)
As a last preparation step we build the required generators.
lookback = 180
steps = 1
delay = 10
batch_size <- 640

test_gen <- generator(
dat_new,
lookback = lookback,
delay = delay,
min_index = 1,
max_index = NULL,
shuffle = FALSE,
steps = steps,
batch_size = batch_size
)

test_gen_pred = generator_predict(
dat_new,
lookback = lookback,
delay = delay,
min_index = 1,
max_index = NULL,
shuffle = FALSE,
steps = steps,
batch_size = batch_size
)

test_steps <- (nrow(dat_new) - 1 - lookback) / batch_size

Test on new data

Before we test our models we check how our naive approach performs to get a reasonable benchmark.
test_naive_method <- function() {
batch_maes <- c()
for (step in 1:test_steps) {
c(samples, targets) %<-% test_gen()
preds <- samples[,dim(samples)[[2]],2:3]
mae <- mean(abs(preds - targets))
batch_maes <- c(batch_maes, mae)
}
mean(batch_maes)
}
(benchmark_test <- test_naive_method())
## [1] 0.003478366

single layer RNN

evaluate_generator(model_simple, test_gen, test_steps)
## $loss ## [1] 0.003198003 ## ##$mean_squared_error
## [1] 0.003198003
##
## $mean_absolute_error ## [1] 0.03449668 Wow - that is catastrophic. The prediction quality is an order of magnitude word worse than the naive benchmark. stacked RNN Will the stacked model perform better? evaluate_generator(model_stacked, test_gen, test_steps) ##$loss
## [1] 0.001995372
##
## $mean_squared_error ## [1] 0.001995372 ## ##$mean_absolute_error
## [1] 0.02614605
Clearly, NO!

conv RNN

My hopes are low, but for the sake of completeness - let’s try the convolutional RNN.
evaluate_generator(model_conv, test_gen, test_steps)
## $loss ## [1] 0.0009481604 ## ##$mean_squared_error
## [1] 0.0009481604
##
## \$mean_absolute_error
## [1] 0.01946304
Nah - that is not good.

pooled

OK, all the individual models were useless. Maybe, their ensemble performs better.
Let’s prepare by getting the targets …
targets_list <- list()
for (step in 1:test_steps) {
c(samples, targets) %<-% test_gen()
targets_list[[step]] <- targets
}

targets_test <- do.call(rbind, targets_list)
… and the model predictions …
pred_simple_test <- predict_generator(model_simple, test_gen_pred, test_steps)
pred_stacked_test <- predict_generator(model_stacked, test_gen_pred, test_steps)
pred_conv_test <- predict_generator(model_conv, test_gen_pred, test_steps)
… and merge them into array.
pred_all_test <- abind::abind(pred_simple_test,
pred_stacked_test,
pred_conv_test,
along = 3)

mean

As a first pooling option we try the mean of the predictions.
pred_mean_test <- apply(pred_all_test, c(1,2), mean)
(mae_mean_test <- mean(abs(pred_mean_test - targets_test)))
## [1] 0.02560716
Hm … that did not help too much.

min/max (conservative)

As a last straw we will try the conservative pooling approach, which worked best on the validation set.
pred_minmax_test <- reshape2::melt(pred_all_test) %>%
group_by(Var1) %>%
summarize(
high = min(1),
low = max(2)
) %>%
ungroup() %>%
select(-Var1) %>%
as.matrix()
(mae_minmax_test <- mean(abs(pred_minmax_test - targets_test)))
## [1] 0.02212562
Oh, my - let’s just consider our attempt of predicting ETH prices with RNNs as failed.

Conclusion

In a nutshell: Forecasting is hard and The Times They Are a-Changin’.
We were able to map some regularities within the training sample, which also seemed to reproduce reasonably in the adjacent time-frame of the validation sample. The crypto-space, however, still evolves rapidly and the last few months have not much in common with the period covered by the training sample. Further, we used data from binance, a very young company with lots of development during the last year, which might have impaired the comparability of training and test data even further.
On an abstract level, the lesson to learn is that systems (might) change and formerly good predictive models can become lousy. Hence, it is crucial to evaluate performance constantly and to update your models.

Some ways to improve fit (not exhaustive)

Our approach failed, but that does not mean that yours must too. Some things that might help you are:
• more data
• time has passed, and more data is available now - use it! (in my opinion, the most important shortcoming of our training data was that it did not cover a sideways market as we can observe in the test data. Hence our models were fully unprepared for this scenerio.)
• we only used a minimal set of features. You can add further useful ones (e.g., data from other exchanges, sentiment in tweets or Reddit, …)
• regularization
• Dropout: amend layer_gru() with dropout and recurrent_dropout. Note: adding dropout to the first recurrent layer does not pay off in my experience
• Max Pooling: add a layer_max_pooling_1d() after layer_conv_1d()
• more models
• the ensemble worked best, there is no reason that you cannot add further models, e.g., different types and architectures, models that use intervals greater than one minute to catch long-running trends, …
• train models to predict the error of your current model instead of the target

Closing Remarks

Although we failed, in the end, I hope you have seen how awesome the keras API for R is and how comfortable you can build even complex neural networks with it. Further, we peaked a little into different RNN architectures and hyperparameter training in general. Maybe, this can be a starting point for your own projects.
If you have any questions or comments, please post them in the comments section. If something is not working as outlined here, please check the package versions you are using. The system I have used was:
sessionInfo()
## R version 3.4.4 (2018-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
##  [3] LC_TIME=de_AT.UTF-8        LC_COLLATE=en_US.UTF-8
##  [5] LC_MONETARY=de_AT.UTF-8    LC_MESSAGES=en_US.UTF-8
##  [7] LC_PAPER=de_AT.UTF-8       LC_NAME=C
## [11] LC_MEASUREMENT=de_AT.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  base
##
## other attached packages:
##  [1] rgenoud_5.8-2.0   bindrcpp_0.2.2    mlrMBO_1.1.2
##  [4] smoof_1.5.1       checkmate_1.8.5   BBmisc_1.11
##  [7] mlr_2.13          ParamHelpers_1.11 keras_2.2.0
## [10] magrittr_1.5      forcats_0.3.0     stringr_1.3.1
## [16] tidyr_0.8.1       tibble_1.4.2      ggplot2_3.0.0
## [19] tidyverse_1.2.1   conflicted_0.1.0  rmarkdown_1.10
##
## loaded via a namespace (and not attached):
##  [1] httr_1.3.1         jsonlite_1.5       viridisLite_0.3.0
##  [4] splines_3.4.4      modelr_0.1.2       assertthat_0.2.0
##  [7] cellranger_1.1.0   yaml_2.2.0         pillar_1.3.0
## [10] backports_1.1.2    lattice_0.20-35    glue_1.3.0
## [13] reticulate_1.10    digest_0.6.16      RColorBrewer_1.1-2
## [16] rvest_0.3.2        colorspace_1.3-2   htmltools_0.3.6
## [19] Matrix_1.2-12      plyr_1.8.4         pkgconfig_2.0.2
## [22] lhs_0.16           broom_0.5.0        misc3d_0.8-4
## [25] haven_1.1.2        scales_1.0.0       parallelMap_1.3
## [28] DiceKriging_1.5.5  whisker_0.3-2      mco_1.0-15.1
## [31] withr_2.1.2        lazyeval_0.2.1     cli_1.0.0
## [34] survival_2.41-3    RJSONIO_1.3-0      crayon_1.3.4
## [40] nlme_3.1-131       xml2_1.2.0         tools_3.4.4
## [43] data.table_1.11.4  hms_0.4.2          plotly_4.8.0
## [46] munsell_0.5.0      compiler_3.4.4     rlang_0.2.2
## [49] plot3D_1.1.1       grid_3.4.4         rstudioapi_0.7
## [52] htmlwidgets_1.2    base64enc_0.1-3    gtable_0.2.0
## [55] R6_2.2.2           tfruns_1.4         lubridate_1.7.4
## [58] knitr_1.20         tensorflow_1.9     zeallot_0.1.0
## [61] bindr_0.1.1        fastmatch_1.1-0    rprojroot_1.3-2
## [64] stringi_1.2.4      parallel_3.4.4     Rcpp_0.12.18
## [67] tidyselect_0.2.4
Besides the R environment the version of keras and tensorflow might also be of interest. To retrieve this information we use the shell:
python3 -c 'import keras; print(keras.__version__)'
## Using TensorFlow backend.
## 2.2.2
and
python3 -c 'import tensorflow; print(tensorflow.__version__)'
## 1.11.0-rc0
Finally, I have to say that knitting rmarkdown with long-running keras model training did not work with RStudio for me (after a few hours it crashed every time). Instead, I have knitted this post directly in the shell with:
Rscript -e 'library(rmarkdown); rmarkdown::render("rmarkdown.Rmd", "html_document")'


text:
# Convex combinations of function values

Suppose $X \subset \mathbb{R}^n$ is a convex set and $f: X \rightarrow \mathbb{R}$ is a function which is bounded from below. Let $x_0 \in X$ be a point.

Define $$a = \inf \{ \alpha f(x_1) + (1 - \alpha) f(x_2) | x_1, x_2 \in X,\,\, \alpha \in [0,1], \text{ and } \alpha x_1 + (1 - \alpha) x_2 = x_0\}$$ and $$b = \inf \{ \sum_{i = 1}^k \alpha_i f(x_i)| k \geq 1, \,\, x_1, \ldots, x_k \in X,\,\, \alpha_i \in [0,1], \,\, \sum_i \alpha_i = 1, \text{ and } \sum_i \alpha_i x_i = x_0\}.$$ Of course, $a \geq b$. My question is whether $a = b$? If $n = 1$ (so $X \subset \mathbb{R}$, I can show that $a = b$. Since I extensively used properties of convex combinations in $\mathbb{R}$, I am not sure it is also true in higher dimensional spaces. I can neither construct a counter-example. Can someone help me out?

Thank you!

I have deleted my previous 1/2 baked answer to this question.

This question is covered by section 17 of Rockafellar's Convex Analysis, in particular Corollary 17.1.5 on p.157. What the OP calls $b$ is what Rockafellar calls $(\operatorname{conv}f)(x_0);$ the OP's formula for $b$ appears at the bottom of R's p.36. The correct formula for $a$ in the case $n>1$ is the conclusion of Cor. 17.1.5: $$(\operatorname{conv}f)(x) = \inf\left\{\sum_{i=1}^{n+1}\lambda_i f(x_i)\,\big| \sum_{i=1}^{n+1} \lambda_i x_i = x\right\}.\tag{*}$$

Here $\operatorname{conv}f$ is the pointwise supremum of all convex functions pointwise less than $f$.

The following counterexamples, based on Ashkan's example, show that the convex combinations in (*) need $n+1$ terms. Let $X\subset\mathbb R^n$ be a regular $n$-simplex, with vertex set $V=\{v_1,\dots,v_{n+1}\}$ and centroid $m=\sum_{i=1}^{n+1} v_i/ (n+1)$. Note that $m$ is at positive distance from each of the facets of $X$. Now define $f(x)=1$ for all $x\in X\backslash V$ and $f(v)=0$ for all $v\in V.$ Clearly $g = \operatorname{conv}f$ is the constant zero function. Now we try to approximate $g(m)$ with $n$-term sums of the form $S=\sum_{i=1}^n \lambda_i f(x_i),$ with $m=\sum_{i=1}^n \lambda_i x_i.$ If there is a $\lambda$ for which $S<\epsilon$, all but $\epsilon$ of the mass of $\lambda$ is concentrated on $V$. But $\lambda$ has only $n$ terms, so all but $\epsilon$ of the mass is concentrated on $n$ of the points of $V$. Hence $m$ is within distance $O(\epsilon)$ of one of the facets of $X$. Since $m$ is at positive distance from the facets, $S$ cannot be made arbitrarily close to $0.$

I found a counter example where $f$ defined on a compact convex set.

Let $f: \Bbb B \to \Bbb R$ and with

$$f(x)= \left\{\begin{matrix} -1 & x\in C \\ 1 & x \notin C \end{matrix}\right.$$

where $X= \Bbb B$ is the closed Unit ball in $\Bbb R^2$ and $C$ is the open left semicircle i.e. $$C= \{ (x,y) \in \Bbb R^2 ~ | ~ x^2 + y^2 = 1, ~~ x <0 \}$$

Then for $x_0 = 0$ observe that $a=0$ and $b < 0$

P.S: Note that

$$b \leq \frac{1}{3}f(\frac{-\sqrt2}{2}, \frac{+\sqrt2}{2} ) + \frac{1}{3} f (\frac{-\sqrt2}{2}, \frac{-\sqrt2}{2}, ) + \frac{1}{3} f(1, 0) = - \frac{1}{3},$$

Maybe you need assume continuity, however still I feel claim is not true!

format: plain


text:
New Bootstrap Method for Autoregressive Models

Title & Authors
New Bootstrap Method for Autoregressive Models
Hwang, Eunju; Shin, Dong Wan;

Abstract
A new bootstrap method combined with the stationary bootstrap of Politis and Romano (1994) and the classical residual-based bootstrap is applied to stationary autoregressive (AR) time series models. A stationary bootstrap procedure is implemented for the ordinary least squares estimator (OLSE), along with classical bootstrap residuals for estimated errors, and its large sample validity is proved. A finite sample study numerically compares the proposed bootstrap estimator with the estimator based on the classical residual-based bootstrapping. The study shows that the proposed bootstrapping is more effective in estimating the AR coefficients than the residual-based bootstrapping.
Keywords
Autoregressive model;stationary bootstrap;residual-based bootstrap;asymptotic property;
Language
English
Cited by
References
1.
Bose, A. (1988). Edgeworth correction by bootstrap in autoregressions, Annals of Statistics, 16, 1709-1722.

2.
Clements, M. P. and Kim, J. H. (2007). Bootstrap prediction intervals for autoregressive time series, Computational Statistics and Data Analysis, 51, 3580-3594.

3.
Freedman, D. A. (1981). Bootstrapping regression models, Annals of Statistics, 9, 1218-1228.

4.
Freedman, D. A. (1984). On bootstrapping two-stage least squares estimates in stationary linear models, Annals of Statistics, 12, 827-842.

5.
Goncalves, S. and Kilian, L. (2004). Bootstrapping autoregressions with conditional heteroskedasticity of unknown form, Journal of Econometrics, 123, 89-120.

6.
Goncalves, S. and Kilian, L. (2007). Asymptotic and bootstrap inference for AR(${\infty}$) processes with conditional heteroskedasticity, Econometric Reviews, 26, 609-641.

7.
Grigoletto, M. (1998). Bootstrap prediction intervals for autoregressions: Some alternatives, International Journal of Forecasting, 14, 447-456.

8.
Hwang, E. and Shin, D.W. (2011). Stationary bootstrapping for non-parametric estimator of nonlinear autoregressive model, Journal of Time Series Analysis, 32, 292-303.

9.
Hwang, E. and Shin, D. W. (2012a). Strong consistency of the stationary bootstrap under $\psi$-weak dependence, Statistics and Probability Letters, 82, 488-495.

10.
Hwang, E. and Shin, D. W. (2012b). Stationary bootstrap for kernel density estimators under $\psi$-weak dependence, Computational Statistics and Data Analysis, 56, 1581-1593.

11.
Kabaila, P. (1993). On bootstrap predictive inference for autoregressive processes, Journal of Time Series Analysis, 14, 473-484.

12.
Kim, J. H. (2004). Bias-correcting bootstrap prediction regions for vector autoregression, Journal of Forecasting, 23, 141-154.

13.
Lahiri, S. N. (1999). On second-order properties of the stationary bootstrap method for studentized statistics, In: Asymptotic, Nonparametrics, and Time Series, (Eds. Ghosh, S.), Marcel Dekker, New York, 683-711.

14.
Nordman, D. J. (2009). A note on the stationary bootstrap's variance, Annals of Statistics, 37, 359-370.

15.
Paparoditis, E. and Politis, D. N. (2005). Bootstrapping unit root tests for autoregressive time series, Journal of the American Statistical Association, 100, 545-553.

16.
Parker, C., Paparoditis, E. and Politis, D. N. (2006). Unit root testing via the stationary bootstrap, Journal of Econometrics, 133, 601-638.

17.
Pascual, L., Romo, J. and Ruiz, E. (2004). Bootstrap predictive inference for ARIMA processes, Journal of Time Series Analysis, 25, 449-465.

18.
Patton, A., Politis, D. N. and White, H. (2009). Correction to "Automatic block-length selection for the dependent bootstrap" by D. Politis and H. White, Econometric Reviews, 28, 372-375.

19.
Politis, D. N. (2003). The Impact of bootstrap methods on time series analysis, Statistical Science, 18, 219-230

20.
Politis, D. N. and Romano, J. P. (1994). The stationary bootstrap, Journal of the American Statistical Association, 89, 1303-1313.

21.
Politis, D. N. and White, H. (2004). Automatic block-length selection for the dependent bootstrap, Econometric Reviews, 23, 53-70.

22.
Swensen, A. R. (2003). Bootstrapping unit root tests for integrated processes, Journal of Time Series Analysis, 24, 99-126.

23.
Thombs, L. A. and Schucany,W. R. (1990). Bootstrap prediction intervals for autoregression, Journal of the American Statistical Association, 95, 486-492.


text:
# DS Interview Study Guide Part II: Software Engineering

This post continues my series on data science interviews. One of the major difficulty of doing data science interviews is that you must show expertise in a wide variety of skills. In particular, I see four key subject areas that you might be asked about during an interview:

1. Statistics
2. Software Engineering/Coding
3. Machine Learning
4. “Soft” Questions

This post focuses on software engineering & coding. It will be primarily a resource for aggregating content that I think you should be familiar with. I will mostly point to outside sources for technical exposition and practice questions.

I’ll link to these as appropriate throughout the post, but I thought it would be helpful to put up front a list of the primary resources that I’ve used when studying for interviews. Some of my favorites are:

• Data Structures and Algorithms in Python, for a good introduction to data structures such as linked lists, arrays, hashmaps, and so on. It also can give you good sense of how to write idiomatic Python code, for building fundamental classes.
• SQLZoo for studying SQL and doing practice questions. I particularly like the “assessments”.
• Cracking the Coding Interview for lots of practice questions organized by subject, and good general advice for the technical interviewing process.

I also use coding websites like LeetCode to practice various problems. I also look on Glassdoor to see what kinds of problems people have been asked.

As always, I’m working to improve this post, so please do leave comments with feedback.

# What Languages Should I Know?

In this section of data science interviews, your are generally asked to implement things in code. So, which language should you do it in? Generally, the best answer is (unsurprisingly) that you should work in Python. The next most popular choice is R; I’m not very familiar with R, so I can’t really speak to it’s capabilities.

There are a few reasons you should work in Python:

1. It’s widely adopted within industry.
2. It has high-quality, popular packages for working with data (see pandas, numpy, scipy, statsmodels, scikit-learn, matplotlib, etc).
3. It bridges the gap between academic work (e.g. using NumPy to build a fast solver for differential equations) and industrial work (e.g. using Django to build webservices).

This is far from an exhaustive list. Anyways, I mostly work in Python. I think it’s a nice language because it is clear and simple to write.

If you want to use another language, you should make sure that you can do everything you need to - this includes reading & writing data, cleaning/munging data, plotting, implementing statistical and machine learning models, and leveraging basic data types like hashmaps and arrays (more on those later).

I think if you wanted to do your interviews in R it would be fine, so long as you can do the above. I would strongly recommend against languages like MATLAB, which are proprietary and not open-source.

Languages like Java can be tricky since they might not have the data-oriented libraries that Python has. For example, I’ve worked profesionally in Scala, and am very comfortable manipulating data via the Spark API within it, but still wouldn’t want to have to use it in an interview; it just isn’t as friendly for general-purpose hacking as Python.

So is Python all you need? Well, not quite. You should also be familiar with SQL for querying databases; we’ll get into that later. I don’t think the dialect you use particularly matters. SQLZoo works with MySQL, which is fine. Familiarity with bash and shell-scripting is useful for a data scientist in their day-to-day work, but generally isn’t asked about in interviews. For the interviews, I’d say if you know one general-purpose language (preferably Python, or R if need be) and SQL, then you’ll be fine.

# General Tips for Coding Interviews

Coding interviews are notorious for being high-stress, so it’s important that you practice in a way that will maximize your comfort during the interview itself - you don’t want to add any unnecessary additional stress into an already difficult situation. There are a wide variety of philosophies and approaches to preparing yourself for and executing a successful interview. I’m going to talk about some points that resonate with me, but I’d also recommend reading Cracking the Coding Interview for a good discussion. Of course, this isn’t the final word on the topic - there are endless resources available online that address this.

## How to Prepare

When preparing for the interview, make sure to practice in an environment similar to the interview environment. There are a few aspects of this to keep in mind.

• Make sure that you replicate the writing environment of the interview. So, if you’ll be coding on a whiteboard, try to get access to a whiteboard to practice. At least practice on a pad of paper, so that you’re comfortable with handwriting code - it’s really quite different than using a text editor. If you’ll be coding in a Google Doc, practice doing that (protip: used a monospaced font). Most places I’ve interviewed at don’t let you evaluate your code to test it, so you have to be prepared for that.
• Time yourself! It’s important to make sure you can do these things in a reasonable amount of time. Generally, these things last 45 minutes per “round” (with multiple rounds for on-site interviews). Focus on being efficient at implementing simple ideas, so that you don’t waste a bunch of time with your syntax and things like that.
• Practice talking. If you practice by coding silently by yourself, then it might feel strange when you’re in the interview and have to talk through your process. The best is if you can have a friend who is familiar with interviewing play the interviewer, so that you can talk to them, get asked questions, etc. You can also record yourself and just talk to the recorder, so that you get practice externalizing your thoughts.

There are some services online that will do “practice” interviews for you. When I was practicing for a software engineer interview with Google, I used Gainlo for this - they were kind of expensive, but you interview with real Google software engineers, which I found helpful.

However, the interviews for a software engineering position at Google are very standardized in format. I haven’t used any of the services that do this for data science, and the interviews you’ll face are so varied. Therefore, I imagine it is harder to do helpful “mock interviews”. If you’ve used any of these services, I’d be very curious to hear about your experience.

## Tips for Interviewing

There are some things it’s important to keep in mind as you do the interview itself.

• Talk about your thought process. Don’t just sit sliently thinking, then go and write something on the board. Let the interviewer into your mind so that they can see how you are thinking about the problem. This is good advice at any point in a technical interview.
• Start with a simple solution you have confidence in. If you know that you can quickly write up a suboptimal solution (in this case, maybe insertion sort), then do that! You can discuss why that solution is sub-optimal, and they will often brainstorm with you about how to improve it. That said, if you are just as confident in writing up something more optimal (say, quicksort) then feel free to jump right to that.
• Sketch out your solution before doing real code. This is not necessary, but sometimes for complicated stuff it’s nice to write out your approach in pseudocode before jumping into real code. This can also help with exposing your thought process to the interviewer, and making sure they’re on board with how you’re thinking about it.
• Think about edge cases. Suppose they ask you to write a function that sorts a list. What if you’re given an empty list? What if you’re given a list of non-comparable things? (In Python, this might be a list of lists.) What does your function do in this case? Is that what you want it to do? There’s no right answer here, but you should definitely be thinking about this and asking the interview how they want the function to behave on these cases.
• Be sure to do a time complexity analysis on your solution. They want to know that you can think about efficiency, so unless they explicitly ask you not to do this, I’d recommend it. We’ll discuss more about what this means below.

For a more thorough discussion of preparation and day-of techniques, I’d recommend Cracking the Coding Interview.

## Tips for Coding

There are few things specifically in how the interviewee writes code that I think are worth mentioning. This kind of stuff usually isn’t a huge deal, but if you write good code, it can show professionalism and help leave a good impression.

• Name your variables well. If the variable is the average number of users per region, use num_users_per_region, or users_per_region, not avg_usr or num_usr. Unlike in mathematics, it’s good to have long, descriptive variables.
• Use built-ins when you can! Python already has functions for sorting, for building cartesian products of lists, for implementing various models (in statsmodels and scikit-learn), and endless other things. It also has some cool data structures already implemented, like the heap and queue. Get to know the itertools module; it has lots of usefull stuff. if you can use these built-ins effectively, it demonstrates skill and knowledge without adding much effort on your part.
• Break things into functions. If one step of your code is sorting a list, and you can’t use the built-in sorted() function, then write a separate function def sort() before you write your main function. This increases both readability and testability of code, and is essential for real-world software.
• Write idiomatic Python. This is a bit less important, but make sure to iterate directly over iterables, don’t do for i in range(len(my_iterable)). Also, familiarize yourself with enumerate and zip and know how to use them. Know how to use list compreshensions, and be aware that you can do a similar thing for dictionaries, sets, and even arguments of functions - for example, you can do max(item for item in l if item % 2 == 0) to find the maximum even number in l. Know how to do string formatting using either .format() for f-strings in Python 3.1

I’m only scratching the surface of how to write good code. It helps to read code that others have written to see what you don’t know. You can also look at code in large open-source libraries.

With all that said, let’s move on to some of the content that might be asked about in these interviews.

# Working with Data

One of the fundamental tasks of a data scientist is to load, manipulate, clean, and visualize data in various formats. I’ll go through some of the basic tasks that I think you should be able to do, and either include or link to Python implementations. If you work in R, or any other language, you should make sure that you can still do these things in your preferred language.

In Python, the key technologies are the packages pandas (for loading, cleaning, and manipulating data), numpy (for efficiently working with unlabeled numeric data), and matplotlib (for plotting and visualizing data).

This tutorial on DataCamp nicely deals with the basics of using pd.read_csv() to load data into Pandas. It is also possible to load from other formats, but in my experience writing to and from comma- or tab-separated plaintext is by far the most common approach for datasets that fit in memory.2

For example, suppose you had the following data in a csv file:

name,age,country,favorite color
steve,7,US,green
jennifer,14,UK,blue
franklin,,UK,black
calvin,22,US,

You can copy and paste this, into Notepad or whatever text editor you like3, and save it as data.csv.

You should be able to

• load in data from text, whether it is separated by commas, tabs, or some other arbitrary character (sometimes things are separated by the “pipe” character |). In this case, you can just do df = pd.read_csv('data.csv') to load it.
• Filter for missing data. If you wanted to find the row(s) where the age is missing, for example, you could do df[df['age'].isnull()]
• Filter for data values. For example, to find people from the US, do df[df['country'] == 'US']
• Replace missing data; use df.fillna(0) to replace missing data with zeros. Think for yourself about how you would want to handle missing data in this case - does it make sense to replace everything with zeros? What would make sense?

Dealing with missing data is, in particular, an important problem, and not one that has an easy answer. Towards Data Science has a decent post on this subject, but if you’re curious, there’s a lot to read about and learn here.

More advanced topics in pandas-fu include using groupby, joining dataframes (this is called a “merge” in pandas, but works the same as a SQL join), and reshaping data.

As I said before, loading and manipulating data is one of the fundamental tasks of a data scientist. You should probably be comfortable doing most or all of these tasks if asked. Pandas can be a bit unintuitive, so I’d recommend practicing if you aren’t already comfortable with it. Doing slicing and reshaping tasks in numpy is also an important skill, so make sure you are comfortable with that as well.

## Visualization

Another essential aspect of data work is visualization. Of course, this is an entire field unto itself; here, I’ll mostly be focusing on the practical aspects of making simple plots. If you want to start to learn more about the overarching principles of the visual representation of data, Tufte’s book is the classic in the field.

In Python, the fundamental tool used for data visualization is the library matplotlib. There exist many other libraries for more complicated visualization tasks, such as seaborn, bokeh, and plotly, but the only one that you really need to be comfortable with (in my opinion) is matplotlib.

You should be comfortable with:

• plotting two lists against one another
• changing the labels on the x- and y-axis of your plot, and adding a title
• changing the x- and y-limits of your plot
• plotting a bar graph
• plotting a histogram
• plotting two curves together, labelling them, and adding a legend

I won’t go through the details here - I’m sure you can find many good guides to each of these online. The matplotlib pyplot tutorial is a good place to start.4

It’s worth noting that you can plot directly from pandas, by doing df.plot(). This just calls out to matplotlib and plots your dataframe; I will often find myself both plotting from the pandas DataFrame.plot() method as well as directly using pyplot.plot(). They work on the same objects, and so you can use them together to make more complicated plots with multiple values plotted.

# Data Structures & Algorithms

Designing and building effective software is predicated on a solid understanding of the basic data structures that are available, and familiarity with the ways that they are employed in common algorithms. For me, learning this material opened up the world of software engineering - it illuminated the inner workings of computer languages. It also helped me understand the pros and cons of various approaches to problems, in ways that I wouldn’t have been able to before.

This subject is fundamental to software engineering interviews, but for data scientists, its importance can vary drastically from role to role. For engineering-heavy roles, this material can make up half or more of the interview, while for more statistician-oriented roles, it might only be very lightly touched upon. You will have to use your judgement to determine to what extent this material is important to you.

I learned this material when I was interviewing by reading the book Data Structures and Algorithms in Python.5 It’s really a great book - it has good, clear explanations of all the important topics, including complexity analysis and some of the basics of the Python language. I can’t recommend it highly enough if you want to get more familiar with this material.6 You can buy it, or look around online for the PDF - it shouldn’t be too hard to find.

## Time and Space Complexity Analysis

Before you begin writing algorithms, you need to know how to analyze their complexity. The “complexity” of an algorithm tells you how the amount of time (or space) that the algorithm takes depends on the size of the input data.

It is formalized using the so-called “big-O” notation. The precise mathematical definition of $\mathcal{O}(n)$ is somewhat confusing, so you can just think of it roughly as meaning that an algorithm that is $\mathcal{O}(n)$ “scales like $n$”; so, if you double the input size, you double the amount of time it takes. If an algorithm is $\mathcal{O}(n^3)$, then, doubling the input size means that you multiply the time it takes by $2^3 = 8$.7 You can see how even a $\mathcal{O}(n^2)$ algorithm wouldn’t work for large data; even if it runs in a reasonable amount of time (say, 5 seconds)for 10,000 points, it would take about 15,000 years to run on 1 billion data points. Obviously, this is no good.

So complexity analysis is critical. You don’t want to settle for a $\mathcal{O}(n^2)$ solution when a $\mathcal{O}(n)$ or $\mathcal{O}(n \log n)$ solution is available. I won’t get into how to do the analysis here, besides saying that I often like to annotate my loops with their complexity when I’m writing things. For example, here’s a (slow) approach to finding the largest k (unique) numbers in a list:

def get_top_k(k, input_list):
top_k = []
for _ in range(k):  # happens k times
remaining = [num for num in input_list if num not in top_k]  # O(n)
if remaining:
top_remaining = max(remaining)  # O(n)
top_k.append(top_remaining)  # O(1)

I know that the outer loop happend k times, and since finding the maximum of a list is $\mathcal{O}(n)$, the total task is $\mathcal{O}(nk)$.8 To learn more about how to do complexity analysis, I’d look at DS&A, Cracking the Coding Interview, or just look around online - I’m sure there are plenty of good resources out there.

You can also consider not just the time of computation, but the amount of memory (space) that your algorithm uses. This is not quite as common as time-complexity analysis, but is still important to be able to do.

A very useful resource for anyone studying for a coding interview is the big-O cheat sheet, which shows the complexity of access, search, insertion, and deletion for various data types, as well as the complexity of searching algorithms, and a lot more. I often use it as a reference, but of course it’s important that you understand why (for example) an array has $\mathcal{O}(n)$ insertion. Just memorizing complexities won’t help you much.

## Arrays & Hashmaps

In my opinion, the two essential data structures for a data scientist to know are the array and the hashmap. In Python, the list type is an array, while the dict type is a hashmap. Since both are used so commonly, you have to know their properties if you want to be able to design efficient algorithms and do your complexity analysis correctly.

Arrays are a data type where a piece of data (like a string) is linked to an index (in Python, this is an integer, starting with 0). I won’t go too deep into the details here, but for arrays, the important thing to know is that getting any element of an array is easy (i.e. doing mylist[5] is $\mathcal{O}(1)$, so it doesn’t depend on the size of the array) but adding elements (particularly in the beginning or middle of the array) is difficult; doing mylist.insert(k, 'foo') is $\mathcal{O}(n-k)$, where $k$ is the position you wish to insert at.9

Arrays are what we usually use when we’re building unordered, unlabelled collections of objects in Python. This is fine, since insertion at the end of an array is fast, and we’re often accessing slices of arrays in a complicated fashion (particularly in numpy). I generally use arrays by default, without thinking too much about it, and it generally works out alright.

Hashmaps also link values to keys, but in this case the key can be anything you want, rather than having to be an ordered set of integers. In Python, you build them by specifying the key and the value, like {'key': 'value'}. Hashmaps are magical in that accessing elements and adding elements are both $\mathcal{O}(1)$.10 Why is this cool? Well, say you wanted to store a bunch of people’s names and ages. You might think to do a list of tuples:

names_ages = [('Peter', 12), ('Kat', 25), ('Jeff', 41)]

Then, if you wanted to find out Jeff’s age, you would have to iterate through the list and find the correct tuple:

for name, age in name_ages:  # happens n times
if name == 'Jeff':
print(f"Jeff's age is {age}")

This is $\mathcal{O}(n)$ - not very efficient. With hashmaps, you can just do

name_ages = {'Peter': 12, 'Kat': 25, 'Jeff': 41}
print(f"Jeff's age is {name_ages['Jeff']}")  # O(1)! Wow!

It might not be obvious how cool this is until you see how to use it in problems. Cracking the Coding Interview has lots of good problems on hashmaps, but I’ll just reproduce some of the classics here. I think it’s worth knowing these, because they really can give you an intuitive sense of when and how hashmaps are valuable.

The first classic hashmap algorithm is counting frequencies of items in a list. That is, given a list, you want to know how many times each item appears. You can do this via the following:

def get_freqs(l):
freqs = {}
for item in l:  # happens O(n) times
if item not in freqs:  # This check is O(1)! Wow!
freqs[item] = 1
else:
freqs[item] += 1  # Also O(1)! Wow!
return freqs

Try and think of how you’d do this without hashmaps. Probably, you’d sort the list, and then look at adjacent values. But sorting is, at best $\mathcal{O}(\log n)$. This solution does it in $\mathcal{O}(n)$!

Another classic problem that is solved with hashmaps is to find all repeated elements in a list. This is really just a variant of the last, where you look for elements that have frequency greater than 1.

def get_repeated(l):
f = get_freqs(l)
return [item for item in f if f[item] > 1]

Now, if you only need one repeated element, you can be efficient and just terminate on the first one you find. For this, we’ll use a set, which is just a dict with values of None. That is to say, sets are also hashmaps. The important thing to know is that adding to them and checking if something is in them are both $\mathcal{O}(1)$.

def get_repeated(l):
items = set()
for item in l:  # happens O(n) times
if item not in items:  # This check is O(1)! Wow!
else:
return(item)
return None  # if this happens, all elements are unique

The last one we’ll do is a bit trickier. You’re given a list of numbers, and a “target”, and your task is to find a pair of numbers in the list that add up to the target. Try and think for yourself how you’d do this - the fact you use hashmaps is a big hint. You should be able to do it in $\mathcal{O}(n)$.

Have you thought about it? When I first encountered this one I had to look up the answer. But here’s how you do it in $\mathcal{O}(n)$:

def get_sum_pair(l, target):
nums_set = set()
for num in l:
other_num = target-num
if other_num in nums_set:
return (num, other_num)
return None

Note that other_num = target-num is the number that you would need to complete the sum pair; using a hashmap, you can check in $\mathcal{O}(1)$ if you’ve already seen it! Wow!

Hopefully you get it - hashmaps are cool. Go on LeetCode, or pop open your favorite data structures book, or even Cracking the Coding Interview, and get some practice with them.

## Sorting & Searching

Sorting and searching are two of the basic tasks you have to be familiar with for any coding interview. You can go into a lot of depth with these, but I’ll stick to the basics here, because that’s what I find most helpful.

### Sorting

Sorting is a nice problem in that the statement of the problem is fairly straightforward; given a list of numbers, reorder the list so that every element is less than or equal to the next. There are a number of approaches to sorting. The naive approach is called insertion sort; for example, it is what most people do when sorting a hand of cards. It has some advantages, but is $\mathcal{O}(n^2)$ in time, and so is not the most efficient available.

The two most common fast sorting algorithms are quicksort and mergesort. They are both $\mathcal{O}(n \log n)$ in time,11 and so scale close-to-linearly with the size of the list. I won’t go into the implementation details here; there are plenty of good discussions of them available on the internet.

When thinking about sorting, it’s also worth considering space complexity - can you sort without needing to carry around a second sorted copy of the list? If so, that’s a significant advantage, especially for larger lists. It’s also worth thinking about worst-case vs. average performance - how does the algorithm perform on a randomly shuffled list, and how does it perform on a list specifically designed to take the maximum number of steps for that algorithm to sort? Quicksort, for example, is actually $\mathcal{O}(n^2)$ in the worst case, but is $\mathcal{O}(n \log n)$ on average. Again, you can look to the big-O cheat sheet to make sure you’re remembering all your complexities correctly.

### Searching

The problem of searching is often stated as given a sorted list l and an object x, find the index at which an element x lives. (You should immediately ask: What should I return if x is not in l?)The name of the game here is binary search. You basically split the list, then if the number is greater than the split, search the top; otherwise, search the bottom. This is an example of a recursive algorithm, so the way it’s written can be a bit opaque to those not used to looking at recursive code. Once I can wrap my head around it, I find it quite elegant. The important thing to know is that this search is $\mathcal{O}(\log n)$, which means that you don’t touch every element in the list - it’s very fast, even for a large list. The key to this is that the list is already sorted - if it’s not sorted, then you’re out of luck; you’ve got to check every element to find x.

There are tons of examples of binary search in Python online, so I won’t put one here. That said, I have found it interesting to see how thinking in terms of binary search can help you in a variety of areas.

For example, suppose you had some eggs, and worked in a 40-story building, and wanted to know the highest floor you could drop the egg off of without it breaking (it’s kind of a dumb example cause the egg would probably break even on the first floor, but pretend it’s a super-tough egg.) You could drop it from the first floor, and see what happens. Say it doesn’t break. Then drop it from the 40th, and see what happens. Say it does break. Then, you bisect and use the midpoint - drop from the 20th floor. If it breaks here, you next try the 10th - if it doesn’t you next try the 30th. This allows you to find the correct floor much faster than trying each floor in succession.

Sorting and searching are fundamental algorithms, and have been well studied for decades. Having a basic fluency in them shows a familiarity with the field of computers science that many employers like to see. In my opinion, you should be able to quickly and easily implement the three sorting algorithms above, and binary search, in Python, or whatever your language of choice is.

# Working with SQL

Finally, let’s talk a bit about SQL. SQL is a tool used to interact with so-called “relational” databases, which just means that each row in a table has certain values (columns), and that those values have the same type for each row (that is, the schema is uniform throughout the table).12 It is not exactly a language, it’s more like a family of languages. There are many “dialects” which all have slight differences, but they behave the same with regards to core functionality; for example, you can do

SELECT column FROM table WHERE columns = 'value'

in any SQL-like language.13 Modern data-storage and -access solutions like Spark and Presto are very different from older databases in their underlying architecture, but still use a SQL dialect for accessing data.

Solving problems in SQL involves thinking in a quite different way than solving a similar problem on an array in Python. There is no real notion of iteration, or at least it’s not easily accessible, so most of the complicated action happens via table joins. I used SQLZoo, and particularly the “assessments”, to practice my SQL and get it up to snuff. LeetCode also has a SQL section (I think they call it “database”).

It’s essential to know SQL as a working data scientist. You’ll almost certainly use it in your day-to-day activities. That said, it’s not always asked in the interviews, so you might clarify with the company whether they will ask you SQL questions.

## A Note on Dialects

There are many dialects of SQL, and changing the dialect changes things like (for example) how you work with dates. It’s worth asking the company you’re interviewing with what dialect they want you to know, if they have one in mind. If you’re just writing SQL on a whiteboard, then I would be surprised if they were picky about this; I would just say something like “here I’d use DATE(table.dt_str) or whatever the string-to-date conversion function is in your dialect”. In this case it’s just details that move around, but the big picture is generally the same for different dialects.

# Conclusion

Coding interviews are stressful. From what I can tell, that’s just the way it is. For me, the best antidote to that is being well-prepared. I think companies are moving more towards constructive, cooperative interview formats, and away from the classic Google brain-teaser kind of questions, which helps with this, but you can still expect to be challenged during these interviews.

Remember to be kind to yourself. You’ll probably fail many times before you succeed. That’s fine, and is what happens to almost everyone. Just keep practicing, and keep learning from your mistakes. Good luck!

1. You should be using Python 3 at this point, but also be familiar with the differences between 2 and 3, and be able to write code in Python 2 if need be.

2. For “big data” stored in the cloud, an efficient format called Parquet is the standard. In my experience, however, it’s uncommon to work with parquet files directly in Pandas; you often read them into a distributed framework like Spark and work with them in that context.

3. The correct answer is, of course, emacs.

4. pyplot is an API within matplotlib that was designed in order to mimic the MATLAB plotting API. It is generally what I use; I begin most of my matplotlib work with from matplotlib import pyplot as plt. I only rarely need to import matplotlib direct, and that’s generally for configuration work.

5. I read the book when preparing for a software engineer interview at Google, so I picked up a lot more than was necessary for a data science interview. I still find the material helpful, however, and it’s nice to be able to demonstrate that you have gone above and beyond in a realm that data scientists sometimes neglect (efficient software design).

6. It goes well beyond what you’ll need for a data science interview, however - it gets into tree structures, graphs (and graph traversal algorithms), and other more advanced topics. I’d recommend focusing on complexity analysis, arrays, and hashmaps as the most important data structures that a data scientist will use day-to-day.

7. This is only approximately true, or rather it is is asymptotically true; this scaling law holds in the limit as $n\rightarrow\infty$

8. It’s a bit weird to use both $n$ and $k$ in your complexity - mathematically, what this means is that we consider them separate variables , and we can take the limit of either one independently from the other. If, for example, you knew that $k = n/4$, so you always wanted the top quarter of the list, then this would be $\mathcal{O}(n^2)$, since $n/4 = \mathcal{O}(n)$

9. I’m glossing over some details here - the numbers I quote above are for a fixed-size array. So, if you build up an array by adding elements at the end, it may seem like you get to just do a bunch of $\mathcal{O}(1)$ .appends, but in reality, you have to occasionally resize the array to make more space, which slows things down to an average append time of $\mathcal{O}(n)$. If you want a list-like type where inserting elements is easy ($\mathcal{O}(1)$) but accessing elements is difficult ($\mathcal{O}(n)$), then you want a linked list. Linked lists aren’t as important for data scientists to use, so I won’t get into them much here.

10. You might wonder why we would ever use an array over a hashmap if hashmaps are strictly superior with respect to their complexity. It’s a good question. The answer is that arrays take up less space (they don’t have to store the keys, only the values) and they are much easier to work with in code (they look cleaner, and are more intuitive for unordered data). Furthermore, if you had a hashmap that linked integers 0 through 10 to strings, and you wanted to change the element at key 5, then you’d have to go through what is currently at keys 5 through 10, and increment their keys by one, so you would end up back at an inefficient insertion algorithm like you have with arrays.

11. This is true on average; see the section below for a discussion of average vs. worst-case complexity.

12. Non-relational database formats, like HBase and NoSQL, basically function like giant hashmaps; they have a single “key”, and then the “value” can contain arbitrary data - you don’t have to have certain columns in there. The advantage of this is flexibility, but the disadvantage is that sorting and filtering are slower because the database doesn’t have a pre-defined schema.

13. Technically, SQL is an ANSI Standard that many different dialects implement - so, to call yourself a SQL dialect, you must have features defined by this standard, like the SELECT, FROM, and WHERE clauses shown above.

Updated:


text:
# Question c7570

Nov 9, 2017

See a solution process below:

#### Explanation:

First, let's write $1 \frac{1}{8} \text{ cups}$ as an improper fraction:

$1 \frac{1}{8} = 1 + \frac{1}{8} = \left(\frac{8}{8} \times 1\right) + \frac{1}{8} = \frac{8}{8} + \frac{1}{8} = \frac{8 + 1}{8} = \frac{9}{8}$

Let's call the amount of fiber we are looking for: $f$

We can then write the relationship and solve for $f$:

$\frac{f}{\frac{9}{8} \text{cups") = (3"g")/(1/2"cups}}$

$\textcolor{red}{\frac{9}{8}} \textcolor{red}{\text{cups") xx f/(9/8"cups") = color(red)(9/8)color(red)("cups") xx(3"g")/(1/2"cups}}$

cancel(color(red)(9/8)color(red)("cups")) xx f/color(red)(cancel(color(black)(9/8"cups"))) = color(red)(9/8)cancel(color(red)("cups")) xx(3"g")/(1/2color(red)(cancel(color(black)("cups"))))#

$f = \textcolor{red}{\frac{9}{8}} \times \frac{3 \text{g}}{\frac{1}{2}}$

$f = \frac{\frac{27 \text{g}}{8}}{\frac{1}{2}}$

$f = \frac{\left(27 \text{g}\right) \times 2}{8 \times 1}$

$f = \frac{54 \text{g}}{8}$

$f = 6.75 \text{g}$

There are close to 7 grams of fiber for the oatmeal left in the box.


text:
# Predicting a cyclic target

I'm familiar with using trigonometric functions to transform cyclic variables for use as features in training a model (most commonly hour of the day or month of the year); I'm now trying to figure out the best way of doing this for using these types of variables as the target for my model. (Imagine a model predicting the month when a particular event is most likely to occur, for instance). Neither using a strictly-increasing representation (so that January, 2018 is close to December, 2017 but very far from January, 2017) nor treating the month as a categorical variable is ideal, but the trigonometric encoding done for features requires both sin and cos parts to have a unique representation of each month, so using one of the two isn't a workable approach either. Is there a better option that I'm missing?

• You can still do the sin/cos trick, but use a bivariate model. Oct 9 '18 at 22:59