Tidymodels による機械学習モデル構築

俺様コードによる tidymodels を用いた機械学習モデル構築のサンプルを記載する。
※基礎的な集計の段階は終了しているという前提で話を進めていく

使用したパッケージの一覧

• rsample
• recipes
• parsnip
• dials
• yardstick

データの分割(訓練/テスト)

データを訓練用とテスト用に分割する。
strata 引数に目的変数を指定する事で目的変数の分布を変えずにデータを分割する事が可能。

• rsample::initial_split
• rsample::analysis / rsample::assesment

library(tidyverse)
library(tidymodels)

# 元データを訓練用とテスト用に分割
# strata を指定する事でクラスの分布を保持したまま分割
lst.splitted <- rsample::initial_split(iris, prop = 0.75, strata = "Species") %>% {
  list(
    train = rsample::analysis(.),
    test  = rsample::assessment(.)
  )
}

クロスバリデーション用データの分割

クロスバリデーション用にデータを分割する。
strata に目的変数を指定しているのは前述した rsample::initial_split と同じ理由。

• rsample::vfold_cv

# クロスバリデーション用の分割を定義
df.cv <- rsample::vfold_cv(
  lst.splitted$train,
  v = 5,
  strata = "Species"
)

前処理の定義

事前に実施した基礎集計の結果やドメイン知識を元にデータの前処理を行う。
いわゆる特徴量エンジニアリングを行う箇所。

※ここで定義している処理はサンプルとして適当に行っているものであり、実施によって精度は下がっていると思われるw

• recipes::recipe
• recipes::step_rm
• recipes::step_log
• recipes::step_center
• recipes::step_scale

# 前処理レシピの定義
recipe <- recipes::recipe(Species ~ ., lst.splitted$train) %>%

  # Sepal.Width を削除
  recipes::step_rm(Sepal.Width) %>%

  # Sepal.Length を対数変換
  recipes::step_log(Sepal.Length) %>%

  # 説明変数を基準化
  recipes::step_center(all_predictors()) %>%
  recipes::step_scale(all_predictors())

モデルの定義

今回は RandomForest でモデルを定義する。
各ハイパーパラメータに parsnip::varying を指定する事で後から複数の値を設定可能。今回はグリッドサーチでパラメータを指定する(後述)。

parsnip::set_engine では指定した engine パッケージ特有のパラメータをセットする事ができる。
下記の例では ranger::ranger 関数の num.threads パラメータ(並列処理に用いるコア数)を指定している。

• parsnip::rand_forest
• parsnip::varying
• parsnip::set_engine

# モデルの定義
model <- parsnip::rand_forest(
  mode = "classification",
  mtry = parsnip::varying(),
  min_n = parsnip::varying(),
  trees = parsnip::varying()
) %>%
  parsnip::set_engine(engine = "ranger", num.threads = 4)

グリッドサーチの定義

3 x 3 x 3 = 27 パターンのハイパーパラメータの組み合わせを定義。
今回はグリッドサーチを行っているが dials::grid_random を用いてランダムサーチを行う事も可能。

• dials::grid_regular

# グリッドサーチ用の組み合わせパターンを定義
df.grid.params <- dials::grid_regular(
  dials::mtry  %>% dials::range_set(c(1, 3)),
  dials::min_n %>% dials::range_set(c(2, 6)),
  dials::trees %>% dials::range_set(c(500, 1500)),
  levels = 3
)

モデルの学習と評価

処理の構造としては下記の 2 重ループとなっている。

1. ハイパーパラメータ一覧のループ
2. クロスバリデーションによる分割のループ

上記 2 重ループの内部で下記の各ステップを実施。

1. モデルの学習
2. 学習済モデルによる予測
3. モデルの評価

上記ループが終了したら 1 段目のループ(ハイパラ一覧のループ)の内部にて下記を実施。

• CV 毎の評価スコアを平均してハイパーパラメータ毎の評価スコアを算出

上記によって得られる結果を評価スコアでソートする事で最も性能の良いハイパーパラメータを特定する。

• recipes::prep
• recipes::bake
• parsnip::fit
• yardstick::metric_set
• yardstick::accuracy
• yardstick::precision
• yardstick::recall
• yardstick::f_meas

df.results <-

  # ハイパーパラメータをモデルに適用
  # ※2020.05 時点で下記のコードは動かない
  merge(df.grid.params, model) %>%
  # 下記のコードに改修(2020.05.19)
  purrr::pmap(df.grid.params, set_args, object = model) %>%

  # ハイパーパラメータの組み合わせごとにループ
  purrr::map(function(model.applied) {

    # クロスバリデーションの分割ごとにループ
    purrr::map(df.cv$splits, model = model.applied, function(df.split, model) {

      # 前処理済データの作成
      df.train <- recipe %>%
        recipes::prep() %>%
        recipes::bake(rsample::analysis(df.split))
      df.test <- recipe %>%
        recipes::prep() %>%
        recipes::bake(rsample::assessment(df.split))

      model %>%

        # モデルの学習
        {
          model <- (.)

          parsnip::fit(model, Species ~ ., df.train)
        } %>%

        # 学習済モデルによる予測
        {
          fit <- (.)

          list(
            train = predict(fit, df.train, type = "class")[[1]],
            test  = predict(fit, df.test,  type = "class")[[1]]
          )
        } %>%

        # 評価
        {
          lst.predicted <- (.)

          # 評価指標の一覧を定義
          metrics <- yardstick::metric_set(
            yardstick::accuracy,
            yardstick::precision,
            yardstick::recall,
            yardstick::f_meas
          )

          # train データでモデルを評価
          df.result.train <- df.train %>%
            dplyr::mutate(
              predicted = lst.predicted$train
            ) %>%
            metrics(truth = Species, estimate = predicted) %>%
            dplyr::select(-.estimator) %>%
            dplyr::mutate(
              .metric = stringr::str_c("train", .metric, sep = "_")
            ) %>%
            tidyr::spread(key = .metric, value = .estimate)

          # test データでモデルを評価
          df.result.test <- df.test %>%
            dplyr::mutate(
              predicted = lst.predicted$test
            ) %>%
            metrics(truth = Species, estimate = predicted) %>%
            dplyr::select(-.estimator) %>%
            dplyr::mutate(
              .metric = stringr::str_c("test", .metric, sep = "_")
            ) %>%
            tidyr::spread(key = .metric, value = .estimate)

          dplyr::bind_cols(
            df.result.train,
            df.result.test
          )
        }
    }) %>%

      # CV 分割全体の平均値を評価スコアとする
      purrr::reduce(dplyr::bind_rows) %>%
      dplyr::summarise_all(mean)
  }) %>%

  # 評価結果とパラメータを結合
  purrr::reduce(dplyr::bind_rows) %>%
  dplyr::bind_cols(df.grid.params) %>%

  # 評価スコアの順にソート(昇順)
  dplyr::arrange(
    desc(test_accuracy)
  ) %>%

  dplyr::select(
    mtry,
    min_n,
    trees,

    train_accuracy,
    train_precision,
    train_recall,
    train_f_meas,

    test_accuracy,
    test_precision,
    test_recall,
    test_f_meas
  )

上記で得られる df.results のサンプル。

mtry	min_n	trees	train_accuracy	train_precision	train_recall	train_f_meas	test_accuracy	test_precision	test_recall	test_f_meas
2	2	500	1.000	1.000	1.000	1.000	0.946	0.949	0.946	0.946
2	4	500	0.991	0.991	0.991	0.991	0.946	0.949	0.946	0.946
3	4	500	0.987	0.987	0.987	0.987	0.946	0.949	0.946	0.946
2	6	500	0.985	0.985	0.985	0.985	0.946	0.949	0.946	0.946
3	6	500	0.982	0.983	0.982	0.982	0.946	0.949	0.946	0.946

ベストモデルの構築

最も性能の良い(=評価スコア最大)モデルを構築する。
この段階では訓練データの全体を用いて学習を行う事に注意。

# 訓練データに前処理を適用
df.train.baked <- recipe %>%
  recipes::prep() %>%
  recipes::bake(lst.splitted$train)

# 最も性能の良いハイパーパラメータを用いたモデルを構築
best_model <- update(
  model,
  mtry  = df.results[1,]$mtry,
  min_n = df.results[1,]$min_n,
  trees = df.results[1,]$trees
) %>%

  # 訓練データ全体を用いてモデルの学習を行う
  parsnip::fit(Species ~ ., df.train.baked)

テストデータによる検証

テストデータを用いた評価を行い、構築したモデルの汎化性能を検証する。

# テストデータに前処理を適用
df.test.baked <- recipe %>%
  recipes::prep() %>%
  recipes::bake(lst.splitted$test)  

# 汎化性能を検証
df.test.baked %>%

  # ベストモデルを用いて予測
  dplyr::mutate(
    predicted = predict(best_model, df.test.baked)[[1]]
  ) %>%

  # 精度(Accuracy)を算出
  yardstick::accuracy(Species, predicted)

まとめ

必要になる度に前に書いたコードを思い出したりコピペしたりでノウハウをまとめられていなかったので、今回は良い機会になった。
機械学習周りは scikit-learn 一択かなと思ってた頃もあったけど tidymodels がいい感じなのでぜひ使う人が増えてくれると嬉し。
他の人がどんな感じでやってるのかも知りたいところ。

時系列のクロスバリデーション

tidymodels に含まれる rsample パッケージの rolling_origin という関数が便利そうだったので試してみた内容を記載してみる。表題の通り時系列データに対するクロスバリデーションを行う。

この関数を知る切っ掛けになったのがこちら。ありがとうございます。

https://blog.hoxo-m.com/entry/2019/06/08/220307

rolling_origin 関数

時系列データのクロスバリデーションに関する考え方はこちらのサイト Rでのナウなデータ分割のやり方: rsampleパッケージによる交差検証 や書籍 前処理大全 に任せるとして、rolling_origin 関数の skip 引数の意味が分かりづらかったので少し試してみる事にする。

まず、今回使用するデータを準備。
"S4248SM144NCEN" というカラム名が何かは知らないけどデータの名称からして恐らく飲酒量とかその辺りだと思われるのでとりあえず "amount" という項目名に変更しておく。

library(tidyverse)
library(tidymodels)

df.drinks <- rsample::drinks %>%
  dplyr::rename(amount = S4248SM144NCEN)

#	date	amount
1	1992-01-01	3,459
2	1992-02-01	3,458
3	1992-03-01	4,002

このデータは 1992/01/01 から 2017/09/01 まで続く月次データであり、全部で 309 件のレコードが存在する。
このデータを下記の要件にて分割する事を考えてみる。

• 学習期間: 24ヶ月
• 検証期間: 12ヶ月
• スライド: 1ヶ月

上記を実現したのがこちら

df.rolling <- rsample::rolling_origin(
  df.drinks,
  initial = 24,
  assess = 12,
  skip = 0,
  cumulative = F
)

df.rolling
# Rolling origin forecast resampling 
# A tibble: 274 x 2
   splits          id      
                
 1 <split [24/12]> Slice001
 2 <split [24/12]> Slice002
 3 <split [24/12]> Slice003
 4 <split [24/12]> Slice004
 5 <split [24/12]> Slice005
 6 <split [24/12]> Slice006
 7 <split [24/12]> Slice007
 8 <split [24/12]> Slice008
 9 <split [24/12]> Slice009
10 <split [24/12]> Slice010
# … with 264 more rows

学習データの 1 番目を取り出してみた。確かに 1992/1/1 を含む 24 ヶ月間のデータとなっている。

df.rolling$splits[[1]] %>%
  rsample::analysis()

# A tibble: 24 x 2
   date       amount
         
 1 1992-01-01   3459
 2 1992-02-01   3458
 3 1992-03-01   4002
 4 1992-04-01   4564
 5 1992-05-01   4221
 6 1992-06-01   4529
 7 1992-07-01   4466
 8 1992-08-01   4137
 9 1992-09-01   4126
10 1992-10-01   4259
# … with 14 more rows

検証用データの 1 番目を取り出してみる。こちらも想定通り 1994 年の 12 ヶ月分のデータとなっている。

df.rolling$splits[[1]] %>%
  rsample::assessment()

# A tibble: 12 x 2
   date       amount
         
 1 1994-01-01   3075
 2 1994-02-01   3377
 3 1994-03-01   4443
 4 1994-04-01   4261
 5 1994-05-01   4460
 6 1994-06-01   4985
 7 1994-07-01   4324
 8 1994-08-01   4719
 9 1994-09-01   4374
10 1994-10-01   4248
11 1994-11-01   4784
12 1994-12-01   4971

2 番目の学習用/検証用のデータも取り出してみる。想定通り 1 ヶ月スライドしたデータとなっている。

df.rolling$splits[[2]] %>%
  rsample::analysis()
# A tibble: 24 x 2
   date       amount
         
 1 1992-02-01   3458
 2 1992-03-01   4002
 3 1992-04-01   4564
 4 1992-05-01   4221
 5 1992-06-01   4529
 6 1992-07-01   4466
 7 1992-08-01   4137
 8 1992-09-01   4126
 9 1992-10-01   4259
10 1992-11-01   4240
# … with 14 more rows

df.rolling$splits[[2]] %>%
  rsample::assessment()
# A tibble: 12 x 2
   date       amount
         
 1 1994-02-01   3377
 2 1994-03-01   4443
 3 1994-04-01   4261
 4 1994-05-01   4460
 5 1994-06-01   4985
 6 1994-07-01   4324
 7 1994-08-01   4719
 8 1994-09-01   4374
 9 1994-10-01   4248
10 1994-11-01   4784
11 1994-12-01   4971
12 1995-01-01   3370

ここで最初に疑問に思ったのが rolling_origin 関数の skip 引数の仕様である。例えば下記のようなデータ仕様を考えてみる。

• 学習期間: 24ヶ月
• 検証期間: 12ヶ月
• スライド: 12ヶ月

つまり 1 年おきにデータを学習/予測(ex. 年の始めに 2 年分のデータを使って直近 1 年間の予測を行う)するようなケースとなる。 この時に直感的には下記のようなコードを想定していた。

df.rolling2 <- rsample::rolling_origin(
  df.drinks,
  initial = 24,
  assess = 12,
  skip = 12,
  cumulative = F
)

12 ヶ月データをスキップするので skip 引数に 12 を指定しているが、これだと下記のように意図した通りにはならない。 2 番目の学習データが 1993/1/1 ではなく 1993/2/1 からのスタートになってしまっており、実際には 13 ヶ月のスキップとなってしまっている。

df.rolling2$splits[[2]] %>%
  rsample::analysis()
# A tibble: 24 x 2
   date       amount
         
 1 1993-02-01   3261
 2 1993-03-01   4160
 3 1993-04-01   4377
 4 1993-05-01   4307
 5 1993-06-01   4696
 6 1993-07-01   4458
 7 1993-08-01   4457
 8 1993-09-01   4364
 9 1993-10-01   4236
10 1993-11-01   4500
# … with 14 more rows

これはそもそも rolling_origin 関数における skip 引数のデフォルト値が 0 であり、この指定によって 1 ヶ月おきのスキップを意図している事による。 つまりそこから更にずらす幅を追加したい時に skip 引数を使いなさいという事なのだろうなと。

これはドキュメントをちゃんと読むとそれっぽい事が記載されていた。ちゃんと読みなさいよと

When skip = 0, the resampling data sets will increment by one position.

あとから考えてみると skip = 0 の指定でスライドしないのであれば無限ループになってしまう orz

なので、上記の仕様であれば下記のコードが正解となる。

df.rolling3 <- rsample::rolling_origin(
  df.drinks,
  initial = 24,
  assess = 12,
  skip = 11, # 11 = 12 - 1
  cumulative = F
)

df.rolling3$splits[[2]] %>%
  rsample::analysis()
# A tibble: 24 x 2
   date       amount
         
 1 1993-01-01   3031
 2 1993-02-01   3261
 3 1993-03-01   4160
 4 1993-04-01   4377
 5 1993-05-01   4307
 6 1993-06-01   4696
 7 1993-07-01   4458
 8 1993-08-01   4457
 9 1993-09-01   4364
10 1993-10-01   4236
# … with 14 more rows

最初から複数月のスライドを意図したデータを作ろうとしていたのが失敗で、スライド幅 1 のケースから順番に試していれば混乱しなかったんだろうなと。難しい

クロスバリデーションによるモデル比較

せっかく rolling_origin の使い方がだいたい分かったので、これを用いて時系列モデルの比較を行ってみる事にする。

まず比較対象となるモデルの一覧を list で作成する。実際には当該リストに含まれるのは学習済みモデルを返す関数である事に注意。

lst.models <- list(
  # ローカルレベルモデル + 確定的季節成分
  d1_ss = function(df.data) {
    library(KFAS)
    model <- SSModel(
      amount ~
        SSMtrend(degree = 1, Q = NA) +
        SSMseasonal(period = 12, Q = 0),
      data = df.data,
      H = NA
    )
    fitSSM(
      model,
      inits = c(1, 1),
      method = "BFGS"
    )
  },

  # ローカルレベルモデル + 確率的季節成分
  d1_sv = function(df.data) {
    library(KFAS)
    model <- SSModel(
      amount ~
        SSMtrend(degree = 1, Q = NA) +
        SSMseasonal(period = 12, Q = NA),
      data = df.data,
      H = NA
    )
    fitSSM(
      model,
      inits = c(1, 1, 1),
      method = "BFGS"
    )
  },

  # ローカル線形トレンドモデル + 確定的季節成分
  d2_ss = function(df.data) {
    library(KFAS)
    model <- SSModel(
      amount ~
        SSMtrend(degree = 2, Q = list(matrix(NA), matrix(NA))) +
        SSMseasonal(period = 12, Q = 0),
      data = df.data,
      H = NA
    )
    fitSSM(
      model,
      inits = c(1, 1, 1),
      method = "BFGS"
    )
  },

  # ローカル線形トレンドモデル + 確率的季節成分
  d2_sv = function(df.data) {
    library(KFAS)
    model <- SSModel(
      amount ~
        SSMtrend(degree = 2, Q = list(matrix(NA), matrix(NA))) +
        SSMseasonal(period = 12, Q = NA),
      data = df.data,
      H = NA
    )
    fitSSM(
      model,
      inits = c(1, 1, 1, 1),
      method = "BFGS"
    )
  }
)

リストの要素を指定してデータを渡せば KFAS による学習済みモデルが返却される。

lst.models$d1_ss(df.drinks)

### ↓以下結果 ###

$optim.out
$optim.out$par
[1]  9.355875 12.207159

$optim.out$value
[1] 2290.374

$optim.out$counts
function gradient 
      29       10 

$optim.out$convergence
[1] 0

$optim.out$message
NULL


$model
Call:
SSModel(formula = amount ~ SSMtrend(degree = 1, Q = NA) + SSMseasonal(period = 12, 
    Q = 0), data = df.data, H = NA)

State space model object of class SSModel

Dimensions:
[1] Number of time points: 309
[1] Number of time series: 1
[1] Number of disturbances: 2
[1] Number of states: 12
Names of the states:
 [1]  level        sea_dummy1   sea_dummy2   sea_dummy3   sea_dummy4   sea_dummy5   sea_dummy6   sea_dummy7   sea_dummy8 
[10]  sea_dummy9   sea_dummy10  sea_dummy11
Distributions of the time series:
[1]  gaussian

Object is a valid object of class SSModel.

こんな感じに予測も可能。

predict(lst.models$d1_ss(df.drinks)$model)

Time Series:
Start = 1 
End = 309 
Frequency = 1 
             fit
  [1,]  2719.484
  [2,]  3162.725
  [3,]  4213.908
  [4,]  4128.718
  [5,]  4758.920 ...

上記を用いて下記のようにモデル毎に評価結果を返す関数を作って実験を行う。
評価指標は CV データ毎に算出した RMSE の平均を用いる事にする。

# h: 予測する範囲 ex. h=1 で 1 期先予測
# skip: スキップ数
rmse_per_model_cv <- function(h, skip = 0) {
  # 時系列 CV の作成
  df.rolling <- rsample::rolling_origin(df.drinks, initial = 24, assess = h,  cumulative = F, skip = skip)

  # 各モデルの適用
  purrr::map(lst.models, function(model) {

    # クロスバリデーション
    purrr::map_dbl(df.rolling$splits, function(split) {
      df.train <- rsample::analysis(split)
      df.test  <- rsample::assessment(split)

      # モデルの学習
      fit <- model(df.train)

      # 学習済みモデルによる予測と RMSE の算出
      # columns: date, amount, predicted
      df.test %>%
        dplyr::mutate(
          predicted = predict(fit$model, n.ahead = h)[, "fit"] %>% as.numeric()
        ) %>%
        yardstick::rmse(amount, predicted) %>%
        .$.estimate
    })
  })
}

モデル毎に 24 ヶ月の学習データを用いて 12 期先までの予測を行う試行を 23 回行い、23 回分の RMSE の平均値をモデルの評価値とする。

# 予測結果の一覧を取得
rmses <- rmse_per_model_cv(h = 12, skip = 11)

# モデル毎に RMSE の平均値を算出
purrr::map_dbl(rmses, ~ mean(.))

   d1_ss    d1_sv    d2_ss    d2_sv 
561.9146 472.5309 442.3523 449.9628

今回のデータだとモデル d2_ss(ローカル線形トレンド+確定的季節成分) の予測精度が最も高いという結果になった。

Iris で k-NN Feature Extraction

前回 は k-NN Feature Extraction を自分で作成したテストデータに適用して試してみた。 そこで今回は Iris データに対して適用してみたいと思う。

関数のコードは https://github.com/you1025/knn_feature_extraction/blob/master/knn_feature_extraction.R

Iris

library(tidyverse)

df.iris <- iris %>%
  tibble::as_tibble()

head(df.iris)

#	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
1	5.1	3.5	1.4	0.2	setosa
2	4.9	3.0	1.4	0.2	setosa
3	4.7	3.2	1.3	0.2	setosa
4	4.6	3.1	1.5	0.2	setosa
5	5.0	3.6	1.4	0.2	setosa
6	5.4	3.9	1.7	0.4	setosa

散布図で可視化

df.iris %>%
  ggplot(aes(Sepal.Length, Sepal.Width)) +
    geom_point(aes(colour = Species), alpha = 2/3)

k-NN Feature Extraction の適用

df.knn_d_columns <- df.iris %>%
  add_knn_d_columns(Species)

head(df.knn_d_columns)

#	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species	setosa_1	versicolor_1	virginica_1
1	5.1	3.5	1.4	0.2	setosa	0.1000000	2.0904545	3.5916570
2	4.9	3.0	1.4	0.2	setosa	0.1414214	1.9131126	3.4799425
3	4.7	3.2	1.3	0.2	setosa	0.1414214	2.0856654	3.6083237
4	4.6	3.1	1.5	0.2	setosa	0.1414214	1.9157244	3.4205263
5	5.0	3.6	1.4	0.2	setosa	0.1414214	2.1424285	3.6166283
6	5.4	3.9	1.7	0.4	setosa	0.3316625	2.0566964	3.4263683

2 軸を変更しながら可視化

# setosa & versicolor
g1 <- df.knn_d_columns %>%
  ggplot(aes(setosa_1, versicolor_1)) +
    geom_point(aes(colour = Species), alpha = 2/3)

# versicolor & virginica
g2 <- df.knn_d_columns %>%
  ggplot(aes(versicolor_1, virginica_1)) +
    geom_point(aes(colour = Species), alpha = 2/3)

# virginica & setosa
g3 <- df.knn_d_columns %>%
  ggplot(aes(virginica_1, setosa_1)) +
    geom_point(aes(colour = Species), alpha = 2/3)

gridExtra::grid.arrange(g1, g2, g3)

どの組み合わせもいい感じに分離できそう。

3 次元プロットも試してみた

plotly::plot_ly(
  x = df.knn_d_columns$setosa_1,
  y = df.knn_d_columns$versicolor_1,
  z = df.knn_d_columns$virginica_1,
  type = "scatter3d",
  mode = "markers",
  color = df.knn_d_columns$Species,
  size = 0
)

いいね

k-NN Feature Extraction を試す

機械学習について色々と調べていたら特徴量エンジニアリングの 1 つとして k-NN Feature Extraction という手法がある事を知ったので試してみた。考え方としては Python: k-NN Feature Extraction について が分かり易く説明してくれていたので参考にさせてもらった。

簡易版の実装は https://github.com/you1025/knn_feature_extraction/blob/master/knn_feature_extraction.R で参照できるとして、以下ではいくつかのサンプルデータに対して当該処理で作成される特徴量がどのようなものとなるのかを確認していく。

サンプル 1: 4 つの正方形(2 クラス)

例としてネットによく挙がってるデータで試してみる。
元データは線形分離が出来ない典型的なデータであり、何らかの非線形変換なしには線形分離は不可能である。そこで k-NN Feature Extraction を実施して新しい特徴量で分布を可視化してみる。

library(tidyverse)

# データの定義
N <- 500
df.data <- tibble(
  x = runif(N, min = -1, max = 1),
  y = runif(N, min = -1, max = 1),
  class = dplyr::if_else(x * y > 0, "a", "b") %>% factor()
)

# 可視化(変換前)
df.data %>%
  ggplot(aes(x, y)) +
    geom_point(aes(colour = class), show.legend = F)

# 特徴量の追加
df.data.knn_d <- df.data %>%
  add_knn_d_columns(col_class = class, k = 1)

# 可視化(変換後)
df.data.knn_d %>%
  ggplot(aes(class_a_1, class_b_1)) +
    geom_point(aes(colour = class), show.legend = F)

変換前

変換後

かなり良い感じに分離できそうなデータに変換されている事が分かる。

サンプル 2: 単位円(2 クラス)

単位円盤上に分布する　π / 3 ごとに 2 つのクラスが切り替わるデータで試してみる。
こちらも単純な直線では分離できないデータである。

library(tidyverse)

# データの定義
N <- 500
df.data <- tibble(
  r = runif(N, 0, 1),
  theta = runif(N, 0, 2 * pi),
  x = r * cos(theta),
  y = r * sin(theta),
  class = cut(theta, breaks = seq(0, 2 * pi, length.out = 7), labels = rep(letters[1:2], 3))
) %>%
  dplyr::select(x, y, class)

# 可視化(変換前)
df.data %>%
  ggplot(aes(x, y)) +
    geom_point(aes(colour = class), show.legend = F)

# 特徴量の追加
df.data.knn_d <- df.data %>%
  add_knn_d_columns(class)

# 可視化(変換後)
df.data.knn_d %>%
  ggplot(aes(class_a_1, class_b_1)) +
  geom_point(aes(colour = class), show.legend = F)

変換前

変換後

こちらもサンプル 1 と同様にかなりきれいに分離可能なデータに変換された。

サンプル 3: 単位円(3 クラス)

単位円盤上に分布する π / 3 ごとに 3 つのクラスが切り替わるデータで試してみる。
こちらはサンプル 2 のデータの 3 クラス版となる。

library(tidyverse)

# データの定義
N <- 500
df.data <- tibble(
  r = runif(N, 0, 1),
  theta = runif(N, 0, 2 * pi),
  x = r * cos(theta),
  y = r * sin(theta),
  class = cut(theta, breaks = seq(0, 2 * pi, length.out = 7), labels = rep(letters[1:3], 2))
) %>%
  dplyr::select(x, y, class)

# 可視化(変換前)
df.data %>%
  ggplot(aes(x, y)) +
  geom_point(aes(colour = class), show.legend = F)

# 特徴量の追加
df.data.knn_d <- df.data %>%
  add_knn_d_columns(col_class = class)

# 可視化(変換後)
df.data.knn_d %>%
  ggplot(aes(class_a_1, class_b_1)) +
  geom_point(aes(colour = class), show.legend = F)

変換前

変換後

クラス a までの距離およびクラス b までの距離の 2 次元で表示したところ、クラス c が平面全体に散らばってしまいクラス a と b の分離には良さそうだがクラス c を含めた線形分離は厳しい事が分かる。そこで 3 次元でのマッピングを試してみる。

plotly::plot_ly(
  x = df.data.knn_d$class_a_1,
  y = df.data.knn_d$class_b_1,
  z = df.data.knn_d$class_c_1,
  type = "scatter3d",
  mode = "markers",
  color = df.data.knn_d$class,
  size = 0
)

画像だと少し分かりづらいかもしれないがクラス a / b / c の各点がそれぞれ y-z 平面 / z-x 平面 / x-y 平面 の上に張り付いており、線形分離が可能な分布となっている。

サンプル 4: 長方形(2 クラス)

比率が 1:1000 の長方形で試してみる。
こちらはサンプル 1 と類似しているデータだが x と y のスケールが 大きく異なっており、この事がどのように影響するのかを確認していく。

library(tidyverse)

# データの定義
N <- 500
df.data <- tibble(
  x = runif(N, -1, 1),
  y = runif(N, -1000, 1000),
  class = dplyr::if_else(x * y > 0, "a", "b") %>% factor()
)

# 可視化(変換前)
df.data %>%
  ggplot(aes(x, y)) +
    geom_point(aes(colour = class), show.legend = F) +
    scale_y_continuous(labels = scales::comma)

# 特徴量の追加
df.data.knn_d <- df.data %>%
  add_knn_d_columns(class)

# 可視化(変換後)
df.data.knn_d %>%
  ggplot(aes(class_a_1, class_b_1)) +
    geom_point(aes(colour = class), show.legend = F)

変換前

変換後

上記は全く分離できそうな雰囲気が無い。
k-means などのユークリッド距離(L2距離)を用いるアルゴリズムと同様に、あまりにもスケールの異なるデータでは数値の大きいデータに引っ張られてまともに処理が出来ない事が原因だと思われる。
仮にサンプル 1 において x, y が共に m(メートル) 単位のデータでありサンプル 4 においては y の単位を cm(センチメートル) に変更したと考えてみると、本質的に同じデータであっても採用する単位によって簡単に結果がおかしくなってしまう事が分かる。この例からも実際の分析に用いる際にはスケーリングが必須となる事が伺える。

まとめ

試してみて実装も比較的容易で使い勝手の良さそうな手法だと感じられる。
問題は処理時間であり、サンプルに用いた 500 件程度のデータでも 10 秒ほどの処理時間が必要となった。データ数 N に対し O(N^2) のアルゴリズムとなるはずなのでこれだとまともなサイズのデータ(万オーダー)には使えない可能性が高い。この点の解消のためには 特徴量エンジニアリング〜fastknn〜 を参考に高速なアルゴリズムへの置き換えが必要だと思われる。

ノイズの付加による見せかけの自己相関

カウントデータの末尾に一定の値(0)の系列が追加される事で独立(だと考えられる)な系列に自己相関が検出されるという事例に遭遇した。
そこでこの事例が一般的な事なのかを確認してみる。

ポアソン過程

まずポアソン分布に従う乱数の系列を考えてみる。このときこの乱数列は自己相関を持たない。

seq.pois <- rpois(100, 2)
forecast::tsdisplay(seq.pois, main = "poisson process")

ポアソン過程 + 0 x10

上記データの末尾に 0 を 10 個追加。ちょっと変化してるっぽい。

c(seq.pois, rep(0, 10)) %>%
  forecast::tsdisplay(main = "poisson process")

Ljung-Box 検定でも自己相関は認められず。

c(seq.pois, rep(0, 10)) %>%
  Box.test(lag = 7, type = "Ljung-Box")

 Box-Ljung test

data:  .
X-squared = 7.3861, df = 7, p-value = 0.3898

ポアソン過程 + 0 x30

さらに 20 個追加。一歩前進。

c(seq.pois, rep(0, 30)) %>%
  forecast::tsdisplay(main = "poisson process + 0x30")

検定では自己相関あり。

c(seq.pois, rep(0, 30)) %>%
  Box.test(lag = 7, type = "Ljung-Box")

 Box-Ljung test

data:  .
X-squared = 44.469, df = 7, p-value = 1.734e-07

ポアソン過程 + 0 x50

さらに 20 個追加して 50 個。ACF プロットを見る限りガッツリ自己相関あり。

c(seq.pois, rep(0, 50)) %>%
  forecast::tsdisplay(main = "poisson process + 0x50")

もはや疑いようのないレベル

c(seq.pois, rep(0, 50)) %>%
  Box.test(lag = 7, type = "Ljung-Box")

 Box-Ljung test

data:  .
X-squared = 96.077, df = 7, p-value < 2.2e-16

シミュレーション

上記のチェックでは偶然の可能性もあるので、何回も繰り返して検定で棄却される(自己相関が存在)確率を求めてみる。

# シミュレーション関数
pois_and_zeros <- function(k, lambda = 2) {
  purrr::map_dbl(1:1000, k = k, lambda = lambda, function(round, k, lambda) {
    seq.pois <- rpois(n = 100, lambda = lambda)
    c(seq.pois, rep(0, k)) %>%
      Box.test(lag = 7, type = "Ljung-Box") %>%
      .$p.value
  })
}

mean(pois_and_zeros( 0) < 0.05) # 0.052
mean(pois_and_zeros(10) < 0.05) # 0.332
mean(pois_and_zeros(20) < 0.05) # 0.994
mean(pois_and_zeros(50) < 0.05) # 1

# 可視化
# 0 の付加を 0〜30 個の 1 刻みで算出
tibble(
  k = seq(0, 30, 1),
  ratio_signifficant = purrr::map_dbl(seq(0, 30, 1), ~ mean(pois_and_zeros(.) < 0.05))
) %>%
  ggplot(aes(k, ratio_signifficant)) +
    geom_point() +
    geom_line() +
    scale_y_continuous(limits = c(0, 1), labels = scales::percent)

完全な乱数列であっても 5% ちょっとの確率で見せかけの自己相関が検出されてしまう。
また k (0 の付加個数)の増加に伴いだいたい k=20 辺りにおいてほぼ 100% の確率で自己相関が検出(Lag=7)されてしまう。

喫煙率の推定

7 年ぶりくらいに新しい記事を書いてみる。
適当に疑問に思った事などを計算して確かめてはすぐ捨てるみたいな感じを続けてたんだけどたまには記録に残していこうかなと。
Qiita にしようか迷ったけど何となくこっちに。

浜田 宏「その問題、数理モデルが解決します の第1章で喫煙率の推定が題材になっていて面白い内容だったので考えた事や試した事などを書いてみる。

書籍では喫煙やドラッグなどのように社会的に望ましくない質問への解答はバイアスが掛かるという事で下記のルールによるアンケートを提案している(回答のランダム化)。

1. コインを投げる(回答者しか分からない)
1. 表が出たら質問へ(正直に)回答
2. 裏が出たら質問の内容に依らず「はい」と回答

上記を行う事で回答者単位での質問への回答は保護されるが、全体としての「はい」を回答する比率は算出可能というのはとても上手いと思う。以前に何かの本で読んでた内容だが今回ちゃんと考えてみて改めて上手いなと感じられる。

その上でもう少し踏み込んで考えてみるという事で、得られた回答から喫煙率の範囲を推定する事を試してみたいと思う。

まず下記の設定の下で考えてみる。

• アンケート回答数: N
• コイン投げで表の出る確率: p
• アンケートで「はい(=喫煙者)」と回答された数: K
• 推定される「はい」の比率(=喫煙率): q

このとき
\begin{align} N \times p \times q + N \times (1 - p) = K \end{align}
により
\begin{align} q = \frac{1}{p} \left\{ \frac{K}{N} - \left( 1 - p \right) \right\} \end{align} となるので N, p, K にそれぞれ 1000, 0.5, 600 を代入する事で喫煙率 q は 20% と算出される。

前提として質問への回答に嘘がないというのは良いとして、コイン投げによる試行にはブレがあり p = 1/2 の下でも表の出る回数は 500 とは限らずその結果によって推定される q は範囲を持つ事になると考えられる。そこで以下では二項分布に従う変数 H を導入し、それを用いて喫煙率 q の範囲を考えてみる。

• アンケート回答数: N
• コイン投げで表の出る確率: p
• コイン投げで表の出た回数(推定値): H
• アンケートで「はい(=喫煙者)」と回答された数: K
• 推定される「はい」の比率(=喫煙率): q

このとき
\begin{align} H \sim Binom(N, p) \end{align} を満たしつつ
\begin{align} q = \frac{K - (N - H)}{H} \end{align} であり、H を 500 に設定すると N, K がそれぞれ 1000, 600 の下でやはり q は 20% という結果を得る。

以下に実際に H の分布を記載してみると 500 を中心に主に 400〜600 付近に分布している事が分かる。

H の 2.5%〜97.5% タイルを算出してみるとその範囲は 469〜531 であり、q が H の単調増加な関数である事からq の範囲(95%)は式(4) の H に 469 と 531 を代入した 14.7%〜24.7% となる。

# H の 2.5%〜97.5% タイル
H.lower <- qbinom(p = 0.025, size = 1000, prob = 1/2)
H.upper <- qbinom(p = 0.975, size = 1000, prob = 1/2)
H <- c(H.lower, H.upper)

# q
(600 - (1000 - H)) / H # 0.1471215 0.2467043

上記が理論上の q の範囲(95%)であり、以下ではシミュレーションによる q の推定を行ってみる。

N <- 1000
K <- 600

# 二項分布に従う乱数の生成
H <- rbinom(n = 1000, size = N, prob = 1/2)

# q の 2.5%〜97.5% タイル
q <- (K - (N - H)) / H
quantile(q, probs = c(0.025, 0.975)) # 0.1471215 0.2467043

# 可視化
tibble::enframe(q) %>%
  ggplot(aes(value)) +
    geom_histogram(aes(y=..density..), binwidth = 0.005) +
    geom_density(fill = "red", alpha = 1/3) +
    geom_vline(xintercept = quantile(q, probs = c(0.025, 0.975)), linetype = 2, colour = "blue", alpha = 1/2)

理論値とそんなに離れておらず精度はそんなに悪くない。 さらにアンケート数を増やす(Nを大きくする)方向で q を算出してみると対数の法則に従って q の 95% 区間が狭まりながら 20% に収束していく様子が見て取れる。

上記で見てきたようにコイン投げの回数が二項分布に従ってばらつく事を加味すると喫煙率 q もばらつく事となり、その分布はアンケート数 N を増加させる事で 20% に近づいていく(大数の法則)。