2019年5月27日月曜日

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〜 を参考に高速なアルゴリズムへの置き換えが必要だと思われる。

2019年5月18日土曜日

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

カウントデータの末尾に一定の値(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)されてしまう。

2019年5月12日日曜日

喫煙率の推定

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% に近づいていく(大数の法則)。