1 R & Python: как их совместно использовать посредством библиотеки reticulate

С недавних пор особая версия Rstudio (версия 1.2.xxx) позволяет выполнять код на R и Python совместно использую оперативную память в программе Rstudio. При этом пакет rmarkdown должен быть не старше версии 1.10. Учитывайте, что “дневная” версия RStudio пока не стабильно выполняет код Python и ее нужно перегружать всякий раз перез запуском. Кроме того, в коде Python он не любит кириллические символы.

Мы применяем пакеты предустановленные в поставке Machine Learning Server 9.х, которая содержит наряду с R также и Python - конкретно сборку Anaconda 4.2 with Python 3.5.

## R Chunk

# Loading packages of R

library('tidyverse')        # An opinionated collection of R Packages designed for Data Science. 

## -- Attaching packages ------------------------------------------------------------------------------- tidyverse 1.2.1 --

## v ggplot2 3.1.0       v purrr   0.3.2  
## v tibble  2.1.1       v dplyr   0.8.0.1
## v tidyr   0.8.3       v stringr 1.4.0  
## v readr   1.3.1       v forcats 0.4.0

## -- Conflicts ---------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library('magrittr')         # A Forward-Pipe Operator for R

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:purrr':
## 
##     set_names

## The following object is masked from 'package:tidyr':
## 
##     extract

# install.packages('reticulate')
library('reticulate')       # Interface to 'Python'           

# reticulate::conda_list()
reticulate::use_python(python = 'C:\\Program Files\\Microsoft\\MLS\\PYTHON_SERVER', required = TRUE)
# reticulate::use_python(python = 'C:\\R\\MLS\\PYTHON_SERVER', required = TRUE)

reticulate::py_config()

## python:         C:\Program Files\Microsoft\MLS\PYTHON_SERVER/python.exe
## libpython:      C:/Program Files/Microsoft/MLS/PYTHON_SERVER/python35.dll
## pythonhome:     C:\PROGRA~1\MIE74D~1\MLS\PYTHON~1
## version:        3.5.2 |Anaconda 4.2.0 (64-bit)| (default, Jul  5 2016, 11:41:13) [MSC v.1900 64 bit (AMD64)]
## Architecture:   64bit
## numpy:          C:\PROGRA~1\MIE74D~1\MLS\PYTHON~1\lib\site-packages\numpy
## numpy_version:  1.14.5
## 
## NOTE: Python version was forced by use_python function

Так же обратите внимание на то, что различные версии **Python* библиотек могут иметь различный синтаксис и из-за объектно-ориентированности этого языка программирования это будет вносить существеннуе изменения в код, что можеит приводить к сбоям уже отлажанных программ. Например, в версии 0.18 библиотеки sklearn, занимающейся машинным обучением у объекта GridSearchCV существовал атрибут grid_scores_. Однако в версии 0.20 он уже безвозвратно изчез и стало необходимым обращаться к атрибуту cv_results_. При этом, составные части последнего также будут изменяться от версии к версии в последующем.

Было бы хорошо, если такие изменения касались только прикладных библиотек. Однако это не так. Существенные изменения в синтаксисе проиходят даже в таких базовых библиотеках как pandas. Например, в версии 0.21 прекратил функционирование метод DataFrame.reindex_axis(). Так, что работая с такой “лапочкой” как Python будьте осторожны, не забывайте, что “сетчатый” - коварная рептилия.

2 Чтение и предобработка данных

Будем работать с данными за 2012-2019 гг., которые были получены из разнообразных баз данных в Интернете.

• Валовой внутренний продукт Казахстана с учетом сезонных колебаний,

• Счет текущих операций (платежный баланс) Казахстана,

• другие макроэкономические показатели Казахстана,

• индексы цен от Всемирного банка по биржевым товарам, составляющих основу экспорта в платежном баланс Казахстана.

## Python Chunk

# Importing Dataframe from MS Excel file

import numpy as np
import pandas as pd
import sklearn
import scipy

print('The used python packages', '\n') 

## The used python packages

print('The numpy version is {}.'.format(np.__version__))

## The numpy version is 1.14.5.

print('The pandas version is {}.'.format(pd.__version__))

## The pandas version is 0.19.2.

print('The scikit-learn version is {}.'.format(sklearn.__version__))

## The scikit-learn version is 0.18.1.

print('The scipy version is {}.'.format(scipy.__version__), '\n')

## The scipy version is 0.19.0.

xl = pd.ExcelFile("ML2019_rodionov_HW8_Data.xlsx")
# print xl.sheet_names
df = xl.parse("Data", skiprows=2)

df.drop(df.tail(1).index, inplace=True) # drop last n rows for "January 2019"
df["Period"] = pd.to_datetime(df["Period"])

df.head()
# df.info()

##       Period  Badrate  LN_Badrate         PoG_CP         PS_CP          NT_CP  \
## 0 2012-01-01      0.0   -6.907755  786052.700000  1.059344e+06  187179.833333   
## 1 2012-02-01      0.0   -6.907755  786052.700000  1.059344e+06  187179.833333   
## 2 2012-03-01      0.0   -6.907755  786052.700000  1.059344e+06  187179.833333   
## 3 2012-04-01      0.0   -6.907755  915496.433333  1.122793e+06  212134.500000   
## 4 2012-05-01      0.0   -6.907755  915496.433333  1.122793e+06  212134.500000   
## 
##          GDP_CP         PoG_05     PS_05         NT_05      ...       \
## 0  2.032577e+06  381646.600000  481124.1  43321.066667      ...        
## 1  2.032577e+06  381646.600000  481124.1  43321.066667      ...        
## 2  2.032577e+06  381646.600000  481124.1  43321.066667      ...        
## 3  2.250424e+06  384440.966667  522228.9  42671.666667      ...        
## 4  2.250424e+06  384440.966667  522228.9  42671.666667      ...        
## 
##    iFATS_OILS     iGRAINS  iOTHERFOOD  iRAW_MATERIAL     iTIMBER  \
## 0  113.562990  129.395382  106.955920     105.724914  111.329464   
## 1  117.047836  130.962494  110.966636     107.847977  111.058334   
## 2  122.580266  132.792130  113.217004     106.463610  109.771439   
## 3  129.173677  129.661085  109.323425     106.255568  110.410051   
## 4  125.242031  131.711678  104.370695     104.597615  111.078986   
## 
##    iOTHERRAWMAT  iFERTILIZERS     iMETMIN    iBASEMET  iPRECIOUSMET  
## 0     99.595886    138.557763  100.502475  101.512249    136.910716  
## 1    104.337191    139.326660  104.005008  105.806667    146.277188  
## 2    102.846232    136.905741  103.543268  104.558683    140.782530  
## 3    101.712306    154.337038  100.953867  100.892434    137.600443  
## 4     97.509712    161.103661   96.628640   97.319598    130.908813  
## 
## [5 rows x 68 columns]

Зафиксируем целевую переменную y, представляющую отношение суммы активов банков, которые оказались в дефолте к общей сумме активов банков второго уровня, выраженное в процентных пунктах (Badrate), и удалим ее из данных - независимых переменных. Отдельно вынесем обозначения периодов Period. Удалим признаки Badrate, LN_Badrate и Period из исходных данных.

## Python Chunk

y = df["Badrate"]
Period = df["Period"]
print(y.describe(include='all'))

## count    84.000000
## mean      0.015251
## std       0.022695
## min       0.000000
## 25%       0.000000
## 50%       0.000000
## 75%       0.025810
## max       0.079678
## Name: Badrate, dtype: float64

df.drop(["Period", "Badrate", "LN_Badrate"], axis=1, inplace=True)

Разделим набор данных на две части: обучающую train_df и тестовую выборки test.df. Последняя содержит 30% периодов наблюдений. Посмотрим размерность получившихся массивов.

## Python Chunk

# Splitting Dataframe into Train & Testing Sets

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    df,
    y,
    test_size=0.3, # 30% of the data will go to the control (test) set
    random_state=17
)

print('X_train.shape = ', X_train.shape)

## X_train.shape =  (58, 65)

print('X_test.shape  =', X_test.shape)

## X_test.shape  = (26, 65)

## R Chunk

# Chart Function for `ggplot2` & `plotly` packages

ShowCaptionsofChart <- function(p, is_ggplot2) {
# Modify the layout of a plotly visualization & Add Mathjax 
  
  # Load required packages
  if (!requireNamespace('ggplot2', quietly = TRUE)) 
    stop('ggplot2 is needed for this function to work. Install it via install.packages("ggplot2")',
         call. = FALSE)
  if (!requireNamespace('plotly', quietly = TRUE)) 
    stop('plotly is needed for this function to work. Install it via install.packages("plotly")',
         call. = FALSE)

  p %>% 
    plotly::layout(annotations =
       list(x = 0, y = 1, text = paste('Graph is created by', ifelse(is_ggplot2 == TRUE, '`ggplot2` & `plotly` packages',
                                                                     '`plotly` package'), 'of R'),
            showarrow = FALSE, xref='paper', yref='paper',
            # xanchor='right', yanchor='auto', xshift=0, yshift=0,
            font=list(size=15, color='gray'))
       ) %>%
        plotly::layout(annotations =
           list(x = 1, y = 0, 
                text = 'Совместная работа языков программирования R и Python в программе RStudio 1.2.xxx',
                showarrow = FALSE, xref = 'paper', yref = 'paper',
                # xanchor='right', yanchor='auto', xshift=0, yshift=0,
                font=list(size=10, color='gray'))
           ) %>%
         plotly::config(mathjax = 'cdn')

}

ShowPlotOfModel <- function(NameOfModel, NameOfColor, MAE) {
# We cannot have R function with signature like function(a = 1, b) since Python function requires that arguments without default values appear before arguments with default values.
    
  # Load required packages
  if (!requireNamespace('reticulate', quietly = TRUE)) 
    stop('ggplot2 is needed for this function to work. Install it via install.packages("reticulate")',
         call. = FALSE)
  if (!requireNamespace('ggplot2', quietly = TRUE)) 
    stop('ggplot2 is needed for this function to work. Install it via install.packages("ggplot2")',
         call. = FALSE)
  if (!requireNamespace('plotly', quietly = TRUE)) 
    stop('plotly is needed for this function to work. Install it via install.packages("plotly")',
         call. = FALSE)

  my_shapes = c(20, 15)
  
  preds <- 
    py$y_test %>%
      data.frame %>%
        rownames_to_column %>% 
          setNames(c('index', 'Badrate')) %>%
            dplyr::mutate( Period = py$Period[index]
                         , Badrate = (exp(py$y_pred) - 0.001) * 100 %>% round(2) 
                         , Legend = 'Predicted' ) %>%
              dplyr::bind_rows(yt %>% dplyr::mutate( Legend = 'Observed' )) %>% 
                dplyr::select(-index)

  errorbars <-
    dplyr::left_join( dplyr::filter(preds, Legend == 'Predicted'), dplyr:: filter(preds, Legend == 'Observed'),
                      by = c('Period' = 'Period') ) %>% 
      dplyr::mutate( lower = if_else(Badrate.x < Badrate.y, Badrate.x, Badrate.y)
                   , upper = if_else(Badrate.x < Badrate.y, Badrate.y, Badrate.x)
                   , Legend = 'Predicted') %>% 
        dplyr::select(Period, lower, upper, Legend)
  
  preds %<>% left_join(errorbars, by = c('Period' = 'Period', 'Legend' = 'Legend'))
  
  MAE = sprintf('MAE = %.2f p.p.', MAE)
  
  R2 = sprintf('R2 = %.2f', py$r2_score(py$y_test, (exp(py$y_pred)-0.001) ))

  # R2 = plotly::TeX( paste('\\displaystyle \\large R^2 =',
  #                         sprintf('%.2f', py$r2_score(py$y_test, (exp(py$y_pred)-0.001) ))) )
  
  p <- 
    ggplot(data = preds) +
      geom_line(data = filter(preds, Legend == 'Observed'), aes(x = Period, y = Badrate), color = 'dodgerblue') +
      geom_point(aes(x = Period, y = Badrate, colour = Legend, shape = Legend
                             , text = paste0('Period:  ', Period,
                                         '<br>Badrate(', Legend, '): ', sprintf('%.2f', Badrate))
                             ), size = 2) +
      geom_linerange(data = filter(preds, Legend == 'Predicted'), aes(x = Period, ymin = lower, ymax = upper)
                     , linetype = 'dashed', color = NameOfColor) +
      scale_colour_manual(values =  c('dodgerblue3', NameOfColor),
                          guide = guide_legend(override.aes = list(linetype = c('solid', 'blank'),
                                                                   shape = my_shapes))) +
      scale_shape_manual(values = my_shapes) +
      scale_x_datetime(date_labels = '%b-%y', date_breaks = '4 months') +
      annotate(geom = 'text', x = as.POSIXct('2013-01-01'), y = max(yt$Badrate) * 0.95, label = MAE, color = 'red') +
      annotate(geom = 'text', x = as.POSIXct('2013-01-01'), y = max(yt$Badrate) * 0.9, label = R2, color = 'red') +
      theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
      labs( title = paste('Предсказание значений Badrate по модели', NameOfModel)
          # , subtitle = 'Graph is created by `ggplot2` & `plotly` packages of R'
          , x = 'Периоды наблюдений', y = 'Badrate, %'
          # , caption = 'Совместная работа языков программирования R и Python в программе RStudio 1.2.xxx'
          )

  plotly::ggplotly(p, width = 900, height = 500, tooltip = c('text') ) %>% 
    ShowCaptionsofChart(is_ggplot2 = TRUE)
  
}

Проведем разведочный анализ зависимой переменной y. Поскольку применение Python библиотеки matplotlib на бета-версии Rstudio 1.2.xxx пока встречает определенные трудности, да и сами ее графические возможности пока не столь изощрены, то мы будем обращаться для каждого графического представления к одной из самых мощных графических библиотек R ggplot2, которая входит в комплект tidyverse. При помощи R библиотеки интерактивной графике plotly (пакет написан на Java) можно получить самую детальную информацию на схеме.

## R Chunk

# EDA of Train Set

library('plotly')           # Create Interactive Web Graphics via 'plotly.js' dc;dc;

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

yt = data.frame(Period = py$Period                             # Преобразование к R типу 'POSIXt' / 'POSIXct'
              , Badrate = round(py$y * 100, 2) %>% as.vector   # Преобразование к R типу 'vector', содержащих 'double'
              )

# Time Series Plot of `Badrate`
p <-
  ggplot(data = yt, aes(x = Period, y = Badrate, colour = Badrate)) +
    geom_line(colour = 'grey60') +
    geom_point(shape = 15, size = 4) +
    scale_color_gradient(low = 'blueviolet', high ='brown1') +
    scale_x_datetime(date_labels = '%b-%y', date_breaks = '4 months') +
    geom_hline(aes(yintercept = mean(yt$Badrate)), color = 'blue') +                                  # Mean
    annotate(geom = 'text', x = as.POSIXct('2013-01-01'), y = mean(yt$Badrate) %>% ceiling,
             label = plotly::TeX(paste0('\\displaystyle \\large \\mu = ', round(mean(yt$Badrate), 2))),
             color = 'red') +
    stat_smooth( color = 'lightcoral', fill = 'lightcoral', method = 'loess' ) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    labs( title = 'Отношение суммы активов банков, которые оказались в дефолте
         к общей сумме активов банков второго уровня'
        , subtitle = 'Graph is created by `ggplot` & `plotly` packages of R'
        , x = 'Периоды наблюдений', y = 'Badrate, %'
        , caption = 'Совместная работа языков программирования R и Python в программе RStudio 1.2.xxx')

plotly::ggplotly(p, width = 900, height = 500, tooltip = c('x', 'y', 'color') ) %>%
  ShowCaptionsofChart(is_ggplot2 = TRUE)

# Combine Histogram and Density plots  of `Badrate`

p <- 
  ggplot(data = yt, aes(x = Badrate, colour = Badrate
)) +
    geom_histogram(aes(y = ..density..), colour = 'black', fill = 'azure', bins = 10)+
    geom_density(alpha = .2, fill = 'lightcoral') +
    geom_vline(aes(xintercept=mean(Badrate)-sd(Badrate)), color = 'blue', linetype = 'dashed') + # Stand. Deviation
    geom_vline(aes(xintercept=mean(Badrate)*1.1), color = 'blue') +                                  # Mean
    geom_vline(aes(xintercept=mean(Badrate)+sd(Badrate)), color = 'blue', linetype = 'dotted') + # Stand. Deviation
    annotate(geom = 'text', x = mean(yt$Badrate)-sd(yt$Badrate), y=0.6,
             label = plotly::TeX('\\displaystyle \\large \\mu - sdt'), color = 'red') +
    annotate(geom = 'text', x = mean(yt$Badrate), y = 0.6,
             label = plotly::TeX(paste0('\\displaystyle \\large \\mu --', round(mean(yt$Badrate), 2))), color = 'red') +
    annotate(geom = 'text', x = mean(yt$Badrate)+sd(yt$Badrate), y=0.6,
             label = plotly::TeX('\\displaystyle \\large \\mu + sdt'), color = 'red') +
    labs( title = 'Кривая плотности распределения переменной Badrate'
        , subtitle = 'Graph is created by `ggplot` & `plotly` packages of R'
        , x = 'Badrate, %', y = 'Плотность распределения'
        , caption = 'Совместная работа языков программирования R и Python в программе RStudio 1.2.xxx')

plotly::ggplotly(p, width = 900, height = 600, tooltip = c('y', 'x', 'color') ) %>% 
  ShowCaptionsofChart(is_ggplot2 = TRUE)

3 Корреляционный анализ признаков и отбор предикторов

Очевидно, что количество переменных для обучения заметно больше количества точек. Необходимо проверить на обучающей выборке X_train в какой мере эти переменные независимы друг от друга.

## Python Chunk

# https://github.com/TheLoneNut/CorrelationMatrixClustering/blob/master/CorrelationMatrixClustering.ipynb
#
# Two pass clustering
# 1-We cluster the corr matrix
#   We sort the survey data according to this clustering
# 2-For cluster bigger than a threshold we cluster those sub-clusters
#   We sort the survey data according to these clustering

import scipy.cluster.hierarchy as sch

cluster_th = 4

X = X_train.corr().values
d = sch.distance.pdist(X)
L = sch.linkage(d, method='complete')
ind = sch.fcluster(L, 0.5*d.max(), 'distance')

columns = [X_train.columns.tolist()[i] for i in list(np.argsort(ind))]
df = X_train.reindex_axis(columns, axis=1)

unique, counts = np.unique(ind, return_counts=True)
counts = dict(zip(unique, counts))

i = 0
j = 0
columns = []
for cluster_l1 in set(sorted(ind)):
    j += counts[cluster_l1]
    sub = df[df.columns.values[i:j]]
    if counts[cluster_l1]>cluster_th:        
        X = sub.corr().values
        d = sch.distance.pdist(X)
        L = sch.linkage(d, method='complete')
        ind = sch.fcluster(L, 0.5*d.max(), 'distance')
        col = [sub.columns.tolist()[i] for i in list((np.argsort(ind)))]
        sub = sub.reindex_axis(col, axis=1)
    cols = sub.columns.tolist()
    columns.extend(cols)
    i = j
df = df.reindex_axis(columns, axis=1)

Для признаков обучающей выборки построим упорядоченную корреляционную матрицу, которая позволит выявить группы взаимосвязанных признаков и выявит структуру данных.

После проведения двухшагового кластерного анализа на признаковом пространстве, попробуем выделить из этих нескольких базовых групп признаков те, которые могли бы стать предикторами в нашей регрессионой модели.

## R Chunk

# Plotting Heat Chart of Correlation Matrix

p <- 
  cor(py$df, method = c('spearman')) %>% 
    round(digits = 4) %>% 
      as.data.frame %>% 
        tibble::rownames_to_column(var = 'Var1') %>%
          tidyr::gather(key = Var2, value = Spearman, -Var1) %>% 
            ggplot2::ggplot(data = ., aes(x = Var1, y = Var2
                                        , text = paste0(Var2, ' ~ ', Var1, ':', Spearman)
                                          )) +
              geom_tile(aes(fill = Spearman)) + 
              scale_fill_gradient2(low = 'firebrick2', 
                                   mid = 'floralwhite', 
                                   high = 'forestgreen', 
                                   midpoint = 0) +
              theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1, size = 9),
                    axis.text.y = element_text(size = 9),
                    plot.title = element_text(size = 13.5, face = 'bold'),
                    legend.title = element_text(face = 'bold', size = rel(1))) +
              labs( title = 'Корреляционная матрица обучающей выборки `X_train`'
                  , subtitle = 'Graph is created by `ggplot` & `plotly` packages of R'
                  , x = NULL, y = NULL, fill = 'Spearman`s rho'
                  , caption = 'Совместная работа языков программирования R и Python в программе RStudio 1.2.xxx')

plotly::ggplotly(p, width = 1000, height = 900, tooltip = c('text') ) %>% 
  ShowCaptionsofChart(is_ggplot2 = TRUE)

Из выделенных кластеров признаков на обучающей выборке отобрано несколько переменных, которые могут стать предикторами для регрессионной модели. Мы создали сокращенный набор предикторов на обучающей выборке X_train_short.

## Python Chunk

X_train_short=X_train[[
'NT_CP',
'GDP_CP',
'PoG_05', 
'PS_05',
'NT_05',
'GDP_05',
'PoG_CP_1YI',
'LoansOI',
'LnLoans',
'IoT'
]] # by ln(Badrate) from ElasticNet

В результате этого шага признаковое пространство сократилось почти в семь раз и мы удалили все мультиколлинеарные признаки, оставив только десять предикторов.

## Python Chunk

# Describing short Train & Testin Sets

print('X_train_short\n', type(X_train_short))

## X_train_short
##  <class 'pandas.core.frame.DataFrame'>

print(X_train_short.shape)

## (58, 10)

X_train_short.describe(include="all")

##                NT_CP        GDP_CP         PoG_05          PS_05  \
## count      58.000000  5.800000e+01      58.000000      58.000000   
## mean   248835.996552  3.619994e+06  413142.893434  590007.948498   
## std     65530.987161  1.148603e+06   50148.063739   72051.666761   
## min    133800.600000  2.032577e+06  353718.600000  481124.100000   
## 25%    212134.500000  2.688893e+06  381597.400000  551808.858333   
## 50%    236454.566667  3.519345e+06  395557.416667  574488.166667   
## 75%    272897.633333  4.356984e+06  441336.200000  622482.058333   
## max    424002.533333  6.339544e+06  567931.284175  766652.300000   
## 
##                NT_05        GDP_05  PoG_CP_1YI       LoansOI    LnLoans  \
## count      58.000000  5.800000e+01   58.000000  5.800000e+01  58.000000   
## mean    67169.929421  1.082419e+06  109.147753  3.818230e+06  16.256230   
## std     37085.410995  1.502854e+05    7.310526  7.825052e+05   0.124005   
## min      9704.266667  9.060918e+05   96.602738  2.341709e+06  15.989490   
## 25%     43321.066667  1.001092e+06  103.923487  3.389576e+06  16.181804   
## 50%     52113.750000  1.042549e+06  110.765305  4.002058e+06  16.314717   
## 75%    108129.633333  1.134661e+06  112.716098  4.182302e+06  16.354485   
## max    133142.313248  1.593780e+06  129.617739  5.302620e+06  16.397433   
## 
##                IoT  
## count    58.000000  
## mean   3281.731994  
## std     759.848043  
## min    1932.366465  
## 25%    2737.798260  
## 50%    3001.441088  
## 75%    3842.402013  
## max    4907.970945

X_test_short =  X_test[X_train_short.columns]
print('\nX_test_short\n', type(X_test_short))

## 
## X_test_short
##  <class 'pandas.core.frame.DataFrame'>

print(X_test_short.shape)

## (26, 10)

Затем отберем из контрольной выборки X_test набор применяемых предикторов в выборке X_test_short. В дальнейшем будем применять только сокращенное признаковое пространство этих предикторов.

4 Определение набора перспективных регрессионый моделей

Таковы основные характеристики оставшихся признаков на обучающей выборке. Оценку будем производить с помощью коэффициента детерминации (R²).

Отберем методы, на которых построим регрессионые модели. Они будут прелдставлять различные типы моделей:

• Обобщенная линейная модель

  • ElasticNet - Регуляризация L₁ и L₂ линейной регрессионной модели (Частотный подход),

  • ARDRegression - Automatic Relevance Determination Regression (ARDRegression) - Регуляризация (L₂) линейной регрессии (Байесовский подход);

• Кластеризация

  • KNeighborsRegressor - Метод ближайших соседей;

• Деревья решений

  • DecisionTreeRegressor - “Деревья решений”,

  • ExtraTreeRegressor - Метод Extremely randomized tree regressor;

• Ансамблевые методы

  • BaggingRegressor - Баггинговый регрессор,

  • GradientBoostingRegressor - Градиентный Бустинговый регрессор.

При этом предсказывать будем не только напрямую значения Badrate - Y1, но также его натуральный логарифм ln(Badrate) - Y2, увеличенный на 0.001. Такой подход позволит,

• с одной стороны, избежать проблем с получение натурального логарифма от нуля,

• а с другой стороны, не допустить появления в прогнозе отрицательных значений Badrate.

Ведь зависимая переменная по сути представляют собой лево-цензуированный вектор неотрицательных значений (не меньше нуля), т.е. с точки зрения эконометрики нам следует построить Tobit модель.

Теперь попробуем настроить эти регрессионные модели с гипер-параметрами по умолчению оценивая их по R² или средняя абсолютная ошибка - mae (от англ. Mean Absolute Error) на контрольной выборке.

## Python Chunk

# Pre-Modeling by `sklearn` package

from pandas import DataFrame
from sklearn.linear_model import ElasticNet
from sklearn.linear_model import ARDRegression
from sklearn.neighbors import KNeighborsRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.tree import ExtraTreeRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.ensemble import BaggingRegressor
from sklearn.metrics import r2_score

# https://code-examples.net/ru/docs/scikit_learn/modules/generated/sklearn.ensemble.extratreesregressor
models = [ElasticNet( max_iter=5000, random_state=17 ), 
          ARDRegression(n_iter=5000),
          KNeighborsRegressor( n_neighbors=3 ), 
          DecisionTreeRegressor( random_state=17 ),
          ExtraTreeRegressor( random_state=17, criterion='mae' ),
          GradientBoostingRegressor( random_state=17 ),
          BaggingRegressor(DecisionTreeRegressor(), random_state=17)
          ]

Y_train = pd.DataFrame({'Y1':y_train, 'Y2':np.log(y_train+0.001)})
Y_test = pd.DataFrame({'Y1':y_test, 'Y2':np.log(y_test+0.001)})

TestModels = pd.DataFrame()
tmp = {}

for model in models:

    m = str(model)
    tmp['Model'] = m[:m.index('(')]    
    
    for i in  [1, 2]:
        np.random.seed(17)
        model.fit( X_train_short, Y_train[['Y%s'%str(i)]].values.ravel() )

        tmp['R2_Y%s'%str(i)] = r2_score(Y_test[['Y%s'%str(i)]], model.predict(X_test_short))

    TestModels = TestModels.append([tmp])

## ElasticNet(alpha=1.0, copy_X=True, fit_intercept=True, l1_ratio=0.5,
##       max_iter=5000, normalize=False, positive=False, precompute=False,
##       random_state=17, selection='cyclic', tol=0.0001, warm_start=False)
## ElasticNet(alpha=1.0, copy_X=True, fit_intercept=True, l1_ratio=0.5,
##       max_iter=5000, normalize=False, positive=False, precompute=False,
##       random_state=17, selection='cyclic', tol=0.0001, warm_start=False)
## ARDRegression(alpha_1=1e-06, alpha_2=1e-06, compute_score=False, copy_X=True,
##        fit_intercept=True, lambda_1=1e-06, lambda_2=1e-06, n_iter=5000,
##        normalize=False, threshold_lambda=10000.0, tol=0.001, verbose=False)
## ARDRegression(alpha_1=1e-06, alpha_2=1e-06, compute_score=False, copy_X=True,
##        fit_intercept=True, lambda_1=1e-06, lambda_2=1e-06, n_iter=5000,
##        normalize=False, threshold_lambda=10000.0, tol=0.001, verbose=False)
## KNeighborsRegressor(algorithm='auto', leaf_size=30, metric='minkowski',
##           metric_params=None, n_jobs=1, n_neighbors=3, p=2,
##           weights='uniform')
## KNeighborsRegressor(algorithm='auto', leaf_size=30, metric='minkowski',
##           metric_params=None, n_jobs=1, n_neighbors=3, p=2,
##           weights='uniform')
## DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
##            max_leaf_nodes=None, min_impurity_split=1e-07,
##            min_samples_leaf=1, min_samples_split=2,
##            min_weight_fraction_leaf=0.0, presort=False, random_state=17,
##            splitter='best')
## DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
##            max_leaf_nodes=None, min_impurity_split=1e-07,
##            min_samples_leaf=1, min_samples_split=2,
##            min_weight_fraction_leaf=0.0, presort=False, random_state=17,
##            splitter='best')
## ExtraTreeRegressor(criterion='mae', max_depth=None, max_features='auto',
##           max_leaf_nodes=None, min_impurity_split=1e-07,
##           min_samples_leaf=1, min_samples_split=2,
##           min_weight_fraction_leaf=0.0, random_state=17, splitter='random')
## ExtraTreeRegressor(criterion='mae', max_depth=None, max_features='auto',
##           max_leaf_nodes=None, min_impurity_split=1e-07,
##           min_samples_leaf=1, min_samples_split=2,
##           min_weight_fraction_leaf=0.0, random_state=17, splitter='random')
## GradientBoostingRegressor(alpha=0.9, criterion='friedman_mse', init=None,
##              learning_rate=0.1, loss='ls', max_depth=3, max_features=None,
##              max_leaf_nodes=None, min_impurity_split=1e-07,
##              min_samples_leaf=1, min_samples_split=2,
##              min_weight_fraction_leaf=0.0, n_estimators=100,
##              presort='auto', random_state=17, subsample=1.0, verbose=0,
##              warm_start=False)
## GradientBoostingRegressor(alpha=0.9, criterion='friedman_mse', init=None,
##              learning_rate=0.1, loss='ls', max_depth=3, max_features=None,
##              max_leaf_nodes=None, min_impurity_split=1e-07,
##              min_samples_leaf=1, min_samples_split=2,
##              min_weight_fraction_leaf=0.0, n_estimators=100,
##              presort='auto', random_state=17, subsample=1.0, verbose=0,
##              warm_start=False)
## BaggingRegressor(base_estimator=DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
##            max_leaf_nodes=None, min_impurity_split=1e-07,
##            min_samples_leaf=1, min_samples_split=2,
##            min_weight_fraction_leaf=0.0, presort=False, random_state=None,
##            splitter='best'),
##          bootstrap=True, bootstrap_features=False, max_features=1.0,
##          max_samples=1.0, n_estimators=10, n_jobs=1, oob_score=False,
##          random_state=17, verbose=0, warm_start=False)
## BaggingRegressor(base_estimator=DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
##            max_leaf_nodes=None, min_impurity_split=1e-07,
##            min_samples_leaf=1, min_samples_split=2,
##            min_weight_fraction_leaf=0.0, presort=False, random_state=None,
##            splitter='best'),
##          bootstrap=True, bootstrap_features=False, max_features=1.0,
##          max_samples=1.0, n_estimators=10, n_jobs=1, oob_score=False,
##          random_state=17, verbose=0, warm_start=False)

TestModels.set_index('Model', inplace=True)

print('Validation on Testing Set', '\n', TestModels)

## Validation on Testing Set 
##                                R2_Y1     R2_Y2
## Model                                        
## ElasticNet                 0.547305  0.656433
## ARDRegression              0.238050  0.370327
## KNeighborsRegressor        0.379940  0.353826
## DecisionTreeRegressor      0.484155  0.944444
## ExtraTreeRegressor         0.889875  0.869585
## GradientBoostingRegressor  0.821687  0.989186
## BaggingRegressor           0.714205  0.946951

print('\n', 'Mean R² Y1 = Badrate    : %.4f' % TestModels.R2_Y1.mean(),
      '\n', 'Mean R² Y2 = ln(Badrate): %.4f' % TestModels.R2_Y2.mean())

## 
##  Mean R² Y1 = Badrate    : 0.5822 
##  Mean R² Y2 = ln(Badrate): 0.7330

Выведем на график качество прогнозирования этими методами по коэффициенту детерминации R² как по собственным значениям Badrate, так и по его натуральному логарифму ln(Badrate).

## R Chunk

# Plotting Graphs with R2 scores by Methods & Outcomes
library('latex2exp')        # Use LaTeX Expressions in Plots

## 
## Attaching package: 'latex2exp'

## The following object is masked from 'package:plotly':
## 
##     TeX

py$TestModels %>%
  round( digits = 4) %>% 
    rownames_to_column %>%
      setNames(c('Methods', 'Badrate', 'ln(Badrate)')) %>% 
        tidyr::gather(key = Transformation, value = R2, -Methods) %>% 
        ggplot(data = ., aes(x = Methods, y = R2, fill = Transformation)) +
          geom_bar(stat = 'identity', show.legend = FALSE) +
          geom_text(aes(label = round(R2, 2)), vjust = 1.6, color = 'azure', size = 3.5) +
          facet_grid(. ~ Transformation, scales = 'free', space = 'free') +
          theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1, size = 9),
                legend.position = 'none') +
          labs( title = latex2exp::TeX('Коэффициент детерминации $R^{2}$ различных моделей')
              , subtitle = 'Graph is created by `ggplot` package of R'
              , x = 'Методы', y = latex2exp::TeX('$R^{2}$')
              , caption = 'Совместная работа языков программирования R и Python в программе RStudio 1.2.xxx')

Наилучшие результаты по коэффициенту детерминации R² на контрольной выборке оказались у среди линейных методов у ElasticNet, по “Деревьям решений” - ExtraTreeRegressor, а для ансамблей отличились оба - BaggingRegressor и GradientBoostingRegressor.

Поскольку значения R² выше по Y2, то предсказывать будем не напрямую значения Badrate, а её превращенную форму - натуральный логарифм ln(Badrate) - Y2, увеличенный на 0.001.

Произведем настройки гипер-параметров у первой пары этих методов при помощи GridSearchCV.

5 Настройка гипер-параметров у модели ElasticNet

Построим наилучшую линейную модель - ElasticNet с регуляризацией L₁ (LASSO) и L₂ (Ridge или по акад. Тихонову). Для этого подберем сетку гипер-параметров α и ρ для этой разновидности регуляризации.

Рассмотрим значения R² в пространстве гипер-параметров ρ - L₁ и α - L₂ при помощи R библиотеки интерактивной графике plotly.

## R Chunk

# Plotting 3D Graph of  Grid Scores

best_score = round(py$grid$best_score_, 4)
alpha = py$grid$best_estimator_$alpha
rho = py$grid$best_estimator_$l1_ratio
  
Scores <- py$scores3D 

plotly::plot_ly(width = 1000, height = 600) %>%
  plotly::add_surface(x = py$alphas3D, y = py$rhos3D, z = Scores) %>%
    plotly::add_trace(name = 'Максимальное значение Score', 
              x = c(alpha, alpha), y = c(rho, rho), z = c(0, best_score),
              type='scatter3d', mode = 'lines+markers',
              marker = list(size = 7, opacity = 0.8), line = list(width = 5, color = 'red2')) %>% 
      plotly::layout(title = plotly::TeX('R^2 \\ \\text{в сетке гипер-параметров линейной модели} \\ ElasticNet'),
               # LaTeX symbols aren’t supported on 3D graphs at the moment: github.com/plotly/plotly.js/issues/608
               scene = list(xaxis = list(type = 'log', title = 'alphas'),
                            yaxis = list(title = 'rhos'),
                            zaxis = list(title = 'Scores'),
                            camera = list(eye = list(x = -1.25, y = -1.25, z = 1.25),
                                        center = list(x = 0, y = 0, z = 0))
                            ) ) %>%
        ShowCaptionsofChart(is_ggplot2 = FALSE)

На трехмерном графике видно, что максимальное значения R², которое выступает в качестве score в этой модели едва достигает 0.7554, при этом значение α - 0.0001, а ρ - 0.7. Теперь настроим с полученными параметрами линейную регрессионную модель на зависимую переменную по обучающей выборке.

## Python Chunk

# Fitting ElasticNet Model

print('Best score (R²) = %.4f' % grid.best_score_, '\n',
      'Best alpha      = ', grid.best_estimator_.alpha, '\n',
      'Best l1_ratio   = ', grid.best_estimator_.l1_ratio)

## Best score (R²) = 0.7554 
##  Best alpha      =  0.0001 
##  Best l1_ratio   =  0.7

regressor = ElasticNet(alpha = grid.best_estimator_.alpha, l1_ratio = grid.best_estimator_.l1_ratio, max_iter=5000)

regressor.fit(X_train_short, Y_train[['Y2']].values.ravel())

## ElasticNet(alpha=0.0001, copy_X=True, fit_intercept=True, l1_ratio=0.7,
##       max_iter=5000, normalize=False, positive=False, precompute=False,
##       random_state=None, selection='cyclic', tol=0.0001, warm_start=False)

Спрогнозируем на контрольной выборке X_test_short по модели ElasticNet зависимую переменную Badrate.

## Python Chunk

# Evaluating ElasticNet Model on the Testing Set

from sklearn.metrics import mean_absolute_error

# Make predictions using the testing set
y_pred = regressor.predict(X_test_short)  # Y2 = ln(Badrate)

# result = pd.DataFrame(data=regressor.coef_)
# result.to_csv('coef.csv', index=False)

# result = pd.DataFrame({'Y1': Y_test[["Y1"]], 'Y2': Y_test[["Y2"]]})
# print(result)

# result = pd.DataFrame(data=y_pred)
# result = pd.DataFrame(data=Y_test)
# print(result, '\n')

# result.to_csv('predictions.csv', index=True)

mae = (mean_absolute_error(Y_test[['Y1']], np.exp(y_pred)-0.001) * 100)
R2 = r2_score(Y_test[['Y1']], np.exp(y_pred)-0.001)

print('Parametric ElasticNet Model - \n',
      'Mean Absolute Error (MAE) on Testing Set (%%): %.4f' % mae, '\n',
      'R² on Testing Set Badrate through Ln(Badrate): %.4f' % R2)

# The coefficients

## Parametric ElasticNet Model - 
##  Mean Absolute Error (MAE) on Testing Set (%): 1.2025 
##  R² on Testing Set Badrate through Ln(Badrate): 0.4373

print('Coefficients for Y2 = ln(Badrate): \n', regressor.coef_)

## Coefficients for Y2 = ln(Badrate): 
##  [-3.49419814e-06  1.98213941e-07  9.90893136e-05  3.93057781e-05
##   8.28432558e-05 -4.20826802e-05 -8.81616015e-03 -3.79760332e-06
##   1.56679368e+01 -1.48363783e-03]

Полученная величина коэффициента детерминации R² на контрольной выборке заметно хуже, чем при обучении.

Посмотрим на интерактивном графике, созданном R библиотеками ggplot2 и plotly спрогнозированные значения по настроенной модели ElasticNet.

## R Chunk

# Plotting Forecast by ElasticNet Model

ShowPlotOfModel('ElasticNet', 'coral', MAE = py$mae)

6 Настройка гипер-параметров у модели ExtraTreeRegressor

Следом обучим другую регрессионную модель - разновидность “деревьев решений” ExtraTreesRegressor. Используем параметр random_state=17 для воспроизводимости результатов.

## Python Chunk

# Tunning GridSearchCV ExtraTreeRegressor

extra_tree = ExtraTreeRegressor(criterion='mae', random_state=17, max_leaf_nodes=None,
                                min_weight_fraction_leaf=0.0, splitter='random')

max_features_ = np.arange(1, X_train_short.shape[1], 1)
max_depth_ = np.arange(1, X_train_short.shape[1], 1)

# Create the ExtraTreeRegressor grid
ET_grid = {'max_features': max_features_,
           'max_depth': max_depth_
          }

# Fit the ExtraTreeRegressor search model
extra_tree_grid = GridSearchCV(estimator=extra_tree, param_grid=ET_grid, cv=kfold, return_train_score=True)

extra_tree_grid.fit( X_train_short, Y_train[['Y2']].values.ravel() )

# The Best Parametres

## GridSearchCV(cv=KFold(n_splits=5, random_state=17, shuffle=True),
##        error_score='raise',
##        estimator=ExtraTreeRegressor(criterion='mae', max_depth=None, max_features='auto',
##           max_leaf_nodes=None, min_impurity_split=1e-07,
##           min_samples_leaf=1, min_samples_split=2,
##           min_weight_fraction_leaf=0.0, random_state=17, splitter='random'),
##        fit_params={}, iid=True, n_jobs=1,
##        param_grid={'max_depth': array([1, 2, 3, 4, 5, 6, 7, 8, 9]), 'max_features': array([1, 2, 3, 4, 5, 6, 7, 8, 9])},
##        pre_dispatch='2*n_jobs', refit=True, return_train_score=True,
##        scoring=None, verbose=0)

print(extra_tree_grid.best_params_, '\n')

# The Best score

## {'max_depth': 9, 'max_features': 5}

print('The best of score: ', extra_tree_grid.best_score_, '\n')

## The best of score:  0.9325221116412763

extra_tree = ExtraTreeRegressor(criterion='mae', random_state=17, max_leaf_nodes=None,
                                max_features=extra_tree_grid.best_estimator_.max_features,
                                max_depth=extra_tree_grid.best_estimator_.max_depth,
                                min_weight_fraction_leaf=0.0, splitter='random')

extra_tree.fit( X_train_short, Y_train[['Y2']].values.ravel() )

# print(pd.DataFrame(extra_tree_grid.cv_results_)[['mean_test_score']])

# Get Test Scores Mean from each Grid Search for 3D Interactive Plotting

## ExtraTreeRegressor(criterion='mae', max_depth=9, max_features=5,
##           max_leaf_nodes=None, min_impurity_split=1e-07,
##           min_samples_leaf=1, min_samples_split=2,
##           min_weight_fraction_leaf=0.0, random_state=17, splitter='random')

max_features3D, max_depth3D = np.meshgrid(max_features_, max_depth_)
# scores3D = pd.DataFrame(extra_tree_grid.grid_scores_).mean_validation_score   # sklearn ver. 0.18 !!!
scores3D = pd.DataFrame(extra_tree_grid.cv_results_)[['mean_test_score']]       # sklearn ver. 0.20 !!!
scores3D = np.array(scores3D).reshape(len(max_features_), len(max_depth_))

Рассмотрим детальнее значения R² в пространстве двух настраиваемых гипер-параметров max_features и max_depth на трехмерном графике.

## R Chunk

# Plotting 3D Graph of  Grid Scores

best_score = round(py$extra_tree_grid$best_score_, 4)
max_features = py$extra_tree_grid$best_estimator_$max_features
max_depth = py$extra_tree_grid$best_estimator_$max_depth
  
Scores <- py$scores3D 

plotly::plot_ly(width = 1000, height = 600) %>%
  plotly::add_surface(x = py$max_features3D, y = py$max_depth3D, z = Scores) %>%
    plotly::add_trace(name = 'Максимальное значение Score', 
              x = c(max_features, max_features), y = c(max_depth, max_depth), z = c(0, best_score),
              type='scatter3d', mode = 'lines+markers',
              marker = list(size = 7, opacity = 0.8), line = list(width = 5, color = 'red2')) %>% 
      plotly::layout(title = plotly::TeX('R^2 \\ \\text{в сетке гипер-параметров линейной модели} \\ ElasticNet'),
               # LaTeX symbols aren’t supported on 3D graphs at the moment: github.com/plotly/plotly.js/issues/608
               scene = list(xaxis = list(# type = 'log', 
                                         title = 'max_features'),
                            yaxis = list(title = 'max_depth'),
                            zaxis = list(title = 'Scores'),
                            camera = list(eye = list(x = 1, y = 1, z = 1.25),
                                        center = list(x = 0, y = 0, z = 0))
                            ) ) %>%
        ShowCaptionsofChart(is_ggplot2 = FALSE)

remove(max_features, max_depth)

Как мы видим выбор двух наилучших значений отбираемых гипер-параметров представляется весьма разумными.

Теперь сделаем прогноз для контрольной выборки и оценим на ней качество модели ExtraTreeRegressor.

## Python Chunk

# Evaluation ExtraTree on the Testing set

y_pred = extra_tree.predict(X_test_short)  # Y2 = ln(Badrate)

mae = (mean_absolute_error(Y_test[["Y1"]], np.exp(y_pred)-0.001) * 100)
R2 = r2_score(Y_test[['Y1']], np.exp(y_pred)-0.001)

print('Non-Parametric ExtraTree Model - \n',
      'Mean Absolute Error (MAE) on Testing Set (%%): %.4f' % mae, '\n',
      'R² on Testing Set Badrate through Ln(Badrate): %.4f' % R2)

## Non-Parametric ExtraTree Model - 
##  Mean Absolute Error (MAE) on Testing Set (%): 0.4503 
##  R² on Testing Set Badrate through Ln(Badrate): 0.8694

Приведем визуальное отображение прогнозных значени, полученных методом ExtraTreeRegressor. Красными точками на графике отмечены прогнозы. К сожалению, не возможно провести границы доверительных интервалов, поскольку предсказания были получены непараметрическим методами.

## Python Chunk

# Plotting Forecast by ExtraTree Model

ShowPlotOfModel('ExtraTree', 'red2', MAE = py$mae)

Полученные в результате прогнозы лежат ближе к наблюдаемым точкам. Качество предсказания по модели “Дерева решений” ExtraTree заметно лучше, чем по линейной модели ElasticNet: MAE = 0.45 p.p. и
R² = 0.87.

7 Обучение регрессионой модели BaggingRegressor

Теперь натренируем ансамблевый метод Баггинга “Деревья решений” на обучающей выборке с параметром random_state=17 для воспроизводимости результатов.

## Python Chunk

# Evaluating BaggingRegressor

max_features = 5

br_model = BaggingRegressor(DecisionTreeRegressor(), random_state=17, max_features=max_features)

br_model.fit( X_train_short, Y_train[['Y2']].values.ravel() )

# Make predictions using the testing set

## BaggingRegressor(base_estimator=DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
##            max_leaf_nodes=None, min_impurity_split=1e-07,
##            min_samples_leaf=1, min_samples_split=2,
##            min_weight_fraction_leaf=0.0, presort=False, random_state=None,
##            splitter='best'),
##          bootstrap=True, bootstrap_features=False, max_features=5,
##          max_samples=1.0, n_estimators=10, n_jobs=1, oob_score=False,
##          random_state=17, verbose=0, warm_start=False)

y_pred = br_model.predict(X_test_short)  # Y2 = ln(Badrate)

mae = (mean_absolute_error(Y_test[['Y1']], np.exp(y_pred)-0.001) * 100)
R2 = r2_score(Y_test[['Y1']], np.exp(y_pred)-0.001)

print('Non-Parametric BaggingRegressor Model - \n',
      'Mean Absolute Error (MAE) on Testing Set (%%): %.4f' % mae, '\n',
      'R² on Testing Set Badrate through Ln(Badrate): %.4f' % R2)

## Non-Parametric BaggingRegressor Model - 
##  Mean Absolute Error (MAE) on Testing Set (%): 0.6453 
##  R² on Testing Set Badrate through Ln(Badrate): 0.7466

Мы получили методом BaggingRegressor на выбранных гипер-параметрах результат лучше, чем по по обобщенной линейной модели с регуляризация L₁ и L₂, но заметно хуже, чем по методу “деревьев решений” ExtraTreeRegressor. Посмотрим на интерактивном графике спрогнозированные значения по настроенной модели BaggingRegressor.

## R Chunk

# Plotting Forecast by BaggingRegressor Model

ShowPlotOfModel('BaggingRegressor', 'violet', MAE = py$mae)

8 Обучение ансамбля GradientBoostingRegressor

Теперь обучим ансамблевый метод “Градиентного бустинга” на обучающей выборке с параметром random_state=17 для воспроизводимости результатов.

## Python Chunk

# Evaluating GradientBoostingRegressor

min_samples_split = 3

gb_model = GradientBoostingRegressor(random_state=17, max_features = max_features,
                                     min_samples_split = min_samples_split)

gb_model.fit( X_train_short, Y_train[['Y2']].values.ravel() )

# Make predictions using the testing set

## GradientBoostingRegressor(alpha=0.9, criterion='friedman_mse', init=None,
##              learning_rate=0.1, loss='ls', max_depth=3, max_features=5,
##              max_leaf_nodes=None, min_impurity_split=1e-07,
##              min_samples_leaf=1, min_samples_split=3,
##              min_weight_fraction_leaf=0.0, n_estimators=100,
##              presort='auto', random_state=17, subsample=1.0, verbose=0,
##              warm_start=False)

y_pred = gb_model.predict(X_test_short)  # Y2 = ln(Badrate)

mae = (mean_absolute_error(Y_test[['Y1']], np.exp(y_pred)-0.001) * 100)
R2 = r2_score(Y_test[['Y1']], np.exp(y_pred)-0.001)

print('Non-Parametric GradientBoostingRegressor Model - \n',
      'Mean Absolute Error (MAE) on Testing Set (%%): %.4f' % mae, '\n',
      'R² on Testing Set Badrate through Ln(Badrate): %.4f' % R2)

## Non-Parametric GradientBoostingRegressor Model - 
##  Mean Absolute Error (MAE) on Testing Set (%): 0.4122 
##  R² on Testing Set Badrate through Ln(Badrate): 0.9024

Мы получили методом “Градиентного бустинга” на установленных гипер-параметрах результат выше, чем по по обобщенной линейной модели с регуляризация L₁ и L₂. Посмотрим на интерактивном графике спрогнозированные значения по настроенной модели GradientBoostingRegressor.

## R Chunk

# Plotting Forecast by GradientBoostingRegressor Model

ShowPlotOfModel('GradientBoostingRegressor', 'gold', MAE = py$mae)

Полученный результат по модели GradientBoostingRegressor следует счесть наилучшим среди полученных в рамках Частотного подхода, так как у нее имеются наименьшая средняя абсолютная ошибка (англ. Mean Absolute Error), равная 0.41 p.p..

Из графика видно, что значительные по абсолютной величине ошибки (около 2 процентных пункта) приходятся на вторую половину 2016 г. Однако уже в начале 2017 г. таких крупных расхождений не наблюдается. В остальном прогнозные точки располагаются на графике вблизи фактических значений.

Особый интерес представляет прогноз для января 2016 г., когда отмечались самые низкие цены на нефть и другие биржевые товары, которые оставляют основу экспорта Казахстана, что неминуемо повлияло как на размер Валового Внутреннего Продукта, так и на состояние Платежного баланса страны. Однако цены довольно быстро стали расти и банковская система относительно безболезненно перенесла эту стрессовую нагрузку.

9 Заключение

Мы построили устойчивую модель прогнозирования уровней отношение суммы активов банков, которые оказались в дефолте к общей сумме активов банков второго уровня - Badrate. Правильное её применение позволит избежать штрафных санкций со стороны Регулятора и возможных репутационных потерь в результате этого.

## Python Chunk

# The End of Session

print(pd.show_versions())

## 
## INSTALLED VERSIONS
## ------------------
## commit: None
## python: 3.5.2.final.0
## python-bits: 64
## OS: Windows
## OS-release: 8.1
## machine: AMD64
## processor: Intel64 Family 6 Model 58 Stepping 9, GenuineIntel
## byteorder: little
## LC_ALL: None
## LANG: None
## LOCALE: None.None
## 
## pandas: 0.19.2
## nose: 1.3.7
## pip: 19.0.3
## setuptools: 27.2.0
## Cython: 0.25.2
## numpy: 1.14.5
## scipy: 0.19.0
## statsmodels: 0.8.0
## xarray: None
## IPython: 5.3.0
## sphinx: 1.5.4
## patsy: 0.4.1
## dateutil: 2.6.0
## pytz: 2017.2
## blosc: None
## bottleneck: 1.2.0
## tables: 3.2.2
## numexpr: 2.6.2
## matplotlib: 2.0.0
## openpyxl: 2.4.1
## xlrd: 1.0.0
## xlwt: 1.2.0
## xlsxwriter: 0.9.6
## lxml: 3.7.3
## bs4: 4.5.1
## html5lib: None
## httplib2: None
## apiclient: None
## sqlalchemy: 1.1.9
## pymysql: None
## psycopg2: None
## jinja2: 2.8
## boto: 2.42.0
## pandas_datareader: 0.2.1
## None

Наконец, посмотрим какие библиотеки применял R.

# The End of Session

devtools::session_info()

## - Session info ----------------------------------------------------------
##  setting  value                       
##  version  R version 3.5.3 (2019-03-11)
##  os       Windows 8.1 x64             
##  system   x86_64, mingw32             
##  ui       RTerm                       
##  language (EN)                        
##  collate  Russian_Russia.1251         
##  ctype    Russian_Russia.1251         
##  tz       Asia/Dhaka                  
##  date     2019-04-01                  
## 
## - Packages --------------------------------------------------------------
##  package     * version date       lib
##  assertthat    0.2.1   2019-03-21 [1]
##  backports     1.1.3   2018-12-14 [1]
##  broom         0.5.1   2018-12-05 [1]
##  callr         3.2.0   2019-03-15 [1]
##  cellranger    1.1.0   2016-07-27 [1]
##  cli           1.1.0   2019-03-19 [1]
##  colorspace    1.4-1   2019-03-18 [1]
##  crayon        1.3.4   2017-09-16 [1]
##  crosstalk     1.0.0   2016-12-21 [1]
##  data.table    1.12.3  2019-03-29 [1]
##  desc          1.2.0   2018-05-01 [1]
##  devtools      2.0.1   2018-10-26 [1]
##  digest        0.6.18  2018-10-10 [1]
##  dplyr       * 0.8.0.1 2019-02-15 [1]
##  evaluate      0.13    2019-02-12 [1]
##  forcats     * 0.4.0   2019-02-17 [1]
##  fs            1.2.7   2019-03-19 [1]
##  generics      0.0.2   2018-11-29 [1]
##  ggplot2     * 3.1.0   2018-10-25 [1]
##  glue          1.3.1   2019-03-12 [1]
##  gtable        0.3.0   2019-03-25 [1]
##  haven         2.1.0   2019-02-19 [1]
##  hms           0.4.2   2018-03-10 [1]
##  htmltools     0.3.6   2017-04-28 [1]
##  htmlwidgets   1.3     2018-09-30 [1]
##  httpuv        1.5.0   2019-03-15 [1]
##  httr          1.4.0   2018-12-11 [1]
##  jsonlite      1.6     2018-12-07 [1]
##  knitr         1.22    2019-03-08 [1]
##  labeling      0.3     2014-08-23 [1]
##  later         0.8.0   2019-02-11 [1]
##  latex2exp   * 0.4.0   2015-11-30 [1]
##  lattice       0.20-38 2018-11-04 [2]
##  lazyeval      0.2.2   2019-03-15 [1]
##  lubridate     1.7.4   2018-04-11 [1]
##  magrittr    * 1.5     2014-11-22 [1]
##  Matrix        1.2-15  2018-11-01 [2]
##  memoise       1.1.0   2017-04-21 [1]
##  mime          0.6     2018-10-05 [1]
##  modelr        0.1.4   2019-02-18 [1]
##  munsell       0.5.0   2018-06-12 [1]
##  nlme          3.1-137 2018-04-07 [2]
##  pillar        1.3.1   2018-12-15 [1]
##  pkgbuild      1.0.3   2019-03-20 [1]
##  pkgconfig     2.0.2   2018-08-16 [1]
##  pkgload       1.0.2   2018-10-29 [1]
##  plotly      * 4.8.0   2018-07-20 [1]
##  plyr          1.8.4   2016-06-08 [1]
##  prettyunits   1.0.2   2015-07-13 [1]
##  processx      3.3.0   2019-03-10 [1]
##  promises      1.0.1   2018-04-13 [1]
##  ps            1.3.0   2018-12-21 [1]
##  purrr       * 0.3.2   2019-03-15 [1]
##  R6            2.4.0   2019-02-14 [1]
##  Rcpp          1.0.1   2019-03-17 [1]
##  readr       * 1.3.1   2018-12-21 [1]
##  readxl        1.3.1   2019-03-13 [1]
##  remotes       2.0.2   2018-10-30 [1]
##  reshape2      1.4.3   2017-12-11 [1]
##  reticulate  * 1.11.1  2019-03-06 [1]
##  rlang         0.3.2   2019-03-21 [1]
##  rmarkdown     1.12    2019-03-14 [1]
##  rprojroot     1.3-2   2018-01-03 [1]
##  rstudioapi    0.10    2019-03-19 [1]
##  rvest         0.3.2   2016-06-17 [1]
##  scales        1.0.0   2018-08-09 [1]
##  sessioninfo   1.1.1   2018-11-05 [1]
##  shiny         1.2.0   2018-11-02 [1]
##  stringi       1.4.3   2019-03-12 [1]
##  stringr     * 1.4.0   2019-02-10 [1]
##  tibble      * 2.1.1   2019-03-16 [1]
##  tidyr       * 0.8.3   2019-03-01 [1]
##  tidyselect    0.2.5   2018-10-11 [1]
##  tidyverse   * 1.2.1   2017-11-14 [1]
##  usethis       1.4.0   2018-08-14 [1]
##  viridisLite   0.3.0   2018-02-01 [1]
##  withr         2.1.2   2018-03-15 [1]
##  xfun          0.5     2019-02-20 [1]
##  xml2          1.2.0   2018-01-24 [1]
##  xtable        1.8-3   2018-08-29 [1]
##  yaml          2.2.0   2018-07-25 [1]
##  source                                
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  Github (Rdatatable/data.table@8603485)
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.0)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.3)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.2)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
##  CRAN (R 3.5.1)                        
## 
## [1] C:/Users/ARodionov/Documents/R/win-library/3.5
## [2] C:/Program Files/R/R-3.5.3/library

sessionInfo()

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Russian_Russia.1251  LC_CTYPE=Russian_Russia.1251   
## [3] LC_MONETARY=Russian_Russia.1251 LC_NUMERIC=C                   
## [5] LC_TIME=Russian_Russia.1251    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] latex2exp_0.4.0   plotly_4.8.0      reticulate_1.11.1
##  [4] magrittr_1.5      forcats_0.4.0     stringr_1.4.0    
##  [7] dplyr_0.8.0.1     purrr_0.3.2       readr_1.3.1      
## [10] tidyr_0.8.3       tibble_2.1.1      ggplot2_3.1.0    
## [13] tidyverse_1.2.1  
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.0        pkgload_1.0.2     jsonlite_1.6     
##  [4] viridisLite_0.3.0 modelr_0.1.4      shiny_1.2.0      
##  [7] assertthat_0.2.1  cellranger_1.1.0  remotes_2.0.2    
## [10] yaml_2.2.0        sessioninfo_1.1.1 pillar_1.3.1     
## [13] backports_1.1.3   lattice_0.20-38   glue_1.3.1       
## [16] digest_0.6.18     promises_1.0.1    rvest_0.3.2      
## [19] colorspace_1.4-1  htmltools_0.3.6   httpuv_1.5.0     
## [22] Matrix_1.2-15     plyr_1.8.4        pkgconfig_2.0.2  
## [25] devtools_2.0.1    broom_0.5.1       haven_2.1.0      
## [28] xtable_1.8-3      scales_1.0.0      processx_3.3.0   
## [31] later_0.8.0       generics_0.0.2    usethis_1.4.0    
## [34] withr_2.1.2       lazyeval_0.2.2    cli_1.1.0        
## [37] crayon_1.3.4      readxl_1.3.1      mime_0.6         
## [40] memoise_1.1.0     evaluate_0.13     ps_1.3.0         
## [43] fs_1.2.7          nlme_3.1-137      xml2_1.2.0       
## [46] pkgbuild_1.0.3    tools_3.5.3       data.table_1.12.3
## [49] prettyunits_1.0.2 hms_0.4.2         munsell_0.5.0    
## [52] callr_3.2.0       compiler_3.5.3    rlang_0.3.2      
## [55] grid_3.5.3        rstudioapi_0.10   htmlwidgets_1.3  
## [58] crosstalk_1.0.0   labeling_0.3      rmarkdown_1.12   
## [61] gtable_0.3.0      reshape2_1.4.3    R6_2.4.0         
## [64] lubridate_1.7.4   knitr_1.22        rprojroot_1.3-2  
## [67] desc_1.2.0        stringi_1.4.3     Rcpp_1.0.1       
## [70] tidyselect_0.2.5  xfun_0.5