2020年9月18日金曜日

さけのわデータで遊んでみる

日本酒に関するデータを Web API で提供してくれている さけのわデータ【生配信】さけのわデータを見て楽しむ で紹介されていたので試してみる。


1. さけのわデータ


利用規約を遵守しつつデータの取得を行う。

1.1 利用規約


サイト上 の利用規約によると 『さけのわデータにて提供されるデータを利用するのに必要なのは 帰属表示 (後述) だけです』 かつ

  • 無料です
  • 商用、非商用を問わず利用できます
  • データを加工して利用できます

との事なので規約に従って帰属表示をちゃんと表示してみる。

〜この記事では さけのわデータ を利用しています〜

ありがとうございます。

1.2 データ取得


httr パッケージで Web API からデータを取得し jsonlite パッケージで json から data.frame への変換を行う。
※データフォーマットに関しては さけのわデータプロジェクト を参照
library(tidyverse)
library(httr)
library(jsonlite)

# 銘柄マスタ
df.brands <-

  # API を叩いて json データを取得
  httr::GET("https://muro.sakenowa.com/sakenowa-data/api/brands") %>%
  httr::content(as = "text", encoding = "utf-8") %>%

  # json を data.frame に変換
  jsonlite::fromJSON() %>%
  .$brands %>%

  # data.frame を tibble に変換
  tibble::as_tibble()

# 醸造所マスタ
df.breweries <-

  # API を叩いて json データを取得
  httr::GET("https://muro.sakenowa.com/sakenowa-data/api/breweries") %>%
  httr::content(as = "text", encoding = "utf-8") %>%

  # json を data.frame に変換
  jsonlite::fromJSON() %>%
  .$breweries %>%

  # data.frame を tibble に変換
  tibble::as_tibble()

# 地域マスタ
df.areas <-

  # API を叩いて json データを取得
  httr::GET("https://muro.sakenowa.com/sakenowa-data/api/areas") %>%
  httr::content(as = "text", encoding = "utf-8") %>%

  # json を data.frame に変換
  jsonlite::fromJSON() %>%
  .$areas %>%

  # data.frame を tibble に変換
  tibble::as_tibble()


# 銘柄ごとのフレーバーデータ
df.flavor_charts <-

  # API を叩いて json データを取得
  httr::GET("https://muro.sakenowa.com/sakenowa-data/api/flavor-charts") %>%
  httr::content(as = "text", encoding = "utf-8") %>%

  # json を data.frame に変換
  jsonlite::fromJSON() %>%
  .$flavorCharts %>%

  # data.frame を tibble に変換
  tibble::as_tibble()

# ランキングデータ
df.rankings <-

  # API を叩いて json データを取得
  httr::GET("https://muro.sakenowa.com/sakenowa-data/api/rankings") %>%
  httr::content(as = "text", encoding = "utf-8") %>%

  # json を data.frame に変換
  jsonlite::fromJSON() %>%

  # 総合ランキングを対象(地域ランキングは今回は対象外)
  .$overall %>%

  # data.frame を tibble に変換
  tibble::as_tibble()


2. クラスタリング


銘柄ごとのフレーバー値(f1〜f6)で k-means によるクラスタリングを行う。

2.1 クラスタ数の推定


cluster パッケージを用いて良さげなクラスタ数を推定する。
※参考: R K-means法のクラスタ数を機械的に決定する方法
set.seed(1025)

clsgap.flavor_charts <- df.flavor_charts %>%

  # フレーバーデータ以外を除去
  dplyr::select(-brandId) %>%

  # 標準化の実施
  scale() %>%

  # 推定に用いる統計量を算出
  # K.max: 推定するクラスタ数の上限
  # B: モンテカルロ法のブートストラップ回数(らしい)
  # ※k-means の収束を安定させるために kmeans 関数の引数 iter.max を大きめに指定している事に注意
  cluster::clusGap(kmeans, K.max = 10, B = 100, iter.max = 100)

# 結果の確認
# 今回はクラスタ数 5 が推定されている
clsgap.flavor_charts
# Clustering Gap statistic ["clusGap"] from call:
# cluster::clusGap(x = ., FUNcluster = kmeans, K.max = 10, B = 100, iter.max = 100)
# B=100 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
#  --> Number of clusters (method 'firstSEmax', SE.factor=1): 5  # 推定されたクラスタ数

2.2 クラスタリングの実施


上記で推定されたクラスタ数 5 を用いて k-means によるクラスタリングを実施。
km.flavor_charts <- df.flavor_charts %>%

  # フレーバーデータ以外を除去
  dplyr::select(-brandId) %>%

  # 標準化の実施
  scale() %>%

  # matrix から tibble へ変換
  tibble::as_tibble() %>%

  # k-means の実施
  # nstart: ランダムにクラスタリングを実施する回数(デフォルトの 1 だと結果が安定しない事が多いので 5〜10 程度を指定する方が良い(と思う))
  kmeans(centers = 5, iter.max = 100, nstart = 5)
km.flavor_charts$centers で取得できる 5 つのクラスタ中心の特徴を確認する。
km.flavor_charts$centers %>%

  # matrix を tibble に変換
  # rownames を指定する事で matrix の表側に指定された内容を列として追加可能
  tibble::as_tibble(rownames = "cluster") %>%
  dplyr::mutate(cluster = factor(cluster)) %>%

  # => long-form
  tidyr::pivot_longer(cols = -cluster, names_to = "feature", values_to = "score") %>%

  # 可視化
  ggplot(aes(feature, cluster)) +
    geom_tile(aes(fill = score), colour = "black", show.legend = F) +
    geom_text(aes(label = round(score, 2)), size = 5, alpha = 1/2) +
    scale_fill_gradient(low = "white", high = "tomato") +

    # 余白を除去
    coord_cartesian(expand = F)

横軸に各フレーバー、縦軸に各クラスタを配置。セル内の数値は標準化処理済みのフレーバー平均値。
各クラスタごとにうまく特徴を抽出できているような感じはある。
※各フレーバー f1〜f6 の中身に関しては さけのわデータプロジェクト を参照

  • クラスタ 1
    • f4(穏やか)が高め
  • クラスタ 2
    • f3(重厚)が高め
    • f6(軽快)が低め
  • クラスタ 3
    • f2(芳醇)が若干高めだが概ね平均的
  • クラスタ 4
    • f1(華やか)が高め
    • f6(軽快)が高め
  • クラスタ 5
    • f5(ドライ)が高め


さらにクラスタごとに各フレーバーの分布を確認する事を考える。
df.flavor_charts %>%

  # クラスタ番号を追加
  dplyr::mutate(cluster = forcats::as_factor(km.flavor_charts$cluster)) %>%

  # => long-form
  tidyr::pivot_longer(cols = dplyr::starts_with("f"), names_to = "feature", values_to = "score") %>%

  # 可視化
  ggplot(aes(feature, score)) +
    geom_boxplot() +
    facet_grid(cluster ~ .)
横軸に各フレーバー、縦軸にオリジナルのフレーバー値(0〜1)を表示。
概ね前述したクラスタ中心に沿って各クラスタごとに特徴的な分布をしている事が見て取れる。
クラスタ 3 は際立った特徴が無く平均的な銘柄を中心とするクラスタである事がよく分かる。



3. 銘柄全体の可視化


6 項目(次元)のフレーバーデータを平面上で表現するために第 1&2 主成分スコアを用いて各銘柄を配置する事を考える。

3.1 主成分分析


銘柄ごとのフレーバー値を用いて主成分分析を実施する。
第 1&2 主成分で約 70% の寄与率であり、それなりに特徴を捉えた表現が出来る可能性あり。
pca.flavor_charts <- df.flavor_charts %>%

  # フレーバーデータ以外を除去
  dplyr::select(-brandId) %>%

  # 主成分分析を実施
  # scale = T で標準化の実施を指定
  prcomp(scale = T)

# 結果を確認
summary(pca.flavor_charts)
# Importance of components:
#                           PC1    PC2    PC3     PC4     PC5     PC6
# Standard deviation     1.6518 1.2178 0.9635 0.71732 0.48521 0.33189
# Proportion of Variance 0.4548 0.2472 0.1547 0.08576 0.03924 0.01836
# Cumulative Proportion  0.4548 0.7019 0.8567 0.94240 0.98164 1.00000
pca.flavor_charts$rotation で取得できる各種成分ごとの固有ベクトルを用いて第 1&2 主成分の特徴を確認する。
pca.flavor_charts$rotation %>%

  # matrix を tibble に変換
  # rownames を指定する事で matrix の表側に指定された内容を列として追加可能
  tibble::as_tibble(rownames = "flavor") %>%

  # => long-form
  tidyr::pivot_longer(cols = -flavor, names_to = "PC", values_to = "score") %>%

  # 可視化
  ggplot(aes(flavor, score)) +

    # 第 3〜6 主成分をグレーで表示
    geom_line(aes(group = PC), colour = "gray", data = function(df) { dplyr::filter(df, PC %in% stringr::str_c("PC", 3:6, sep = "")) }, alpha = 1/2) +
    geom_point(colour = "gray", data = function(df) { dplyr::filter(df, PC %in% stringr::str_c("PC", 3:6, sep = "")) }, alpha = 1/2) +

    # 第 1&2 主成分のみを強調表示
    geom_line(aes(group = PC, colour = PC), data = function(df) { dplyr::filter(df, PC %in% stringr::str_c("PC", 1:2, sep = "")) }, size = 1, alpha = 2/3) +
    geom_point(aes(colour = PC), data = function(df) { dplyr::filter(df, PC %in% stringr::str_c("PC", 1:2, sep = "")) }, size = 3)
ざっくり第 1 主成分が軽さ、第 2 主成分が辛口の度合いという感じ?(適当)
※各フレーバー f1〜f6 の中身に関しては さけのわデータプロジェクト を参照


3.2 銘柄の配置


算出した第 1&2 主成分を用いて銘柄を配置していく。
df.flavor_charts %>%

  # 各主成分スコアを追加
  dplyr::bind_cols(tibble::as_tibble(pca.flavor_charts$x)) %>%

  # クラスタ番号を追加
  dplyr::mutate(cluster = forcats::as_factor(km.flavor_charts$cluster)) %>%

  # 銘柄の名称を追加
  dplyr::left_join(df.brands, by = c("brandId" = "id")) %>%

  # ランキング(トップ 50)を追加
  dplyr::left_join(df.rankings, by = "brandId") %>%

  dplyr::select(
    brand_name = name,
    PC1,
    PC2,
    cluster,
    rank
  ) %>%

  # 可視化
  ggplot(aes(PC1, PC2)) +

    # 銘柄の名称を追加
    geom_text(aes(label = brand_name, colour = cluster), family = "Osaka", size = 3, alpha = 2/3) +

    # クラスタごとの中心(平均値)を追加
    geom_point(
      aes(x = PC1, y = PC2, fill = cluster),
      data = function(df) {
        dplyr::group_by(df, cluster) %>%
          dplyr::summarise(dplyr::across(dplyr::starts_with("PC"), mean))
      },
      size = 5,
      shape = 24 # 表示を▲に指定
    ) +

    # ランキング(トップ 50)を追加
    geom_text(aes(label = rank), data = function(df) { dplyr::filter(df, dplyr::between(rank, 1, 50)) }, size = 5, alpha = 2/3)

横軸を第 1 主成分、縦軸を第 2 主成分として各銘柄の位置を表現した。
5 つの三角(▲)は各クラスタの重心(平均値)の位置を表しており、黒文字の数字はランキング対象銘柄の位置を表す。
概ね各クラスタごとに配置が分離されており、第 1&2 主成分でそれなりに特徴を捉える事が出来ているものと思われる。
ランキングのトップ 50 に入るような人気銘柄の多くがクラスタ 3 と 4 の領域に存在しており、 辛口を避け(=下半分)軽い飲み心地(=やや右寄り)の尖りすぎない(=原点に近い)銘柄が好まれる傾向にあるのではないかと考えられる。


ランキングをトップ 10 に絞ってみると上記の傾向がより顕著となり、仮説としてはそれなりに悪くないのかも。



4. 決定木


ランクを数量と見做して回帰を行う事の是非は一旦見ないふりをして各ランクを分類するための決定木を構築する。
今回は ggparty パッケージを使用して ggplot 上でグラフを作成する(決定木をキレイに可視化できるようになって嬉しい)。
※ggparty の詳細に関しては ggparty: Graphic Partying を参照

各エッジ上に分割の条件となる数値が表示されるが、現状では簡単に丸め処理を行う方法が無い(恐らく)ので ggparty 作者が作成した非公式の add_splitvar_breaks_index_new 関数をダウンロードして用いている。
※詳細は Number of decimal places on edges of a decision-tree plot with ggparty 参照
library(ggparty)

# github 上の add_splitvar_breaks_index_new 関数をロードする
# エッジ上のラベルに表示する数値を丸めるために用いる
fnc <- readr::read_file("https://raw.githubusercontent.com/martin-borkovec/ggparty/martin/R/add_splitvar_breaks_index_new.R") %>%
  parse(text = .)
eval(fnc)


df.flavor_charts %>%

  # 各主成分スコアを追加
  dplyr::bind_cols(tibble::as_tibble(pca.flavor_charts$x)) %>%

  # クラスタ番号を追加
  dplyr::mutate(cluster = forcats::as_factor(km.flavor_charts$cluster)) %>%

  # ランキングの追加
  dplyr::inner_join(df.rankings, by = "brandId") %>%

  dplyr::select(-c(
    brandId,
    score
  )) %>%

  # 決定木を構築
  rpart::rpart(rank ~ ., data = ., minbucket = 5)%>%
  partykit::as.party() %>%

  {
    pt <- (.)

    # ノード分割に用いる数値の表示を丸める
    rounded_labels <- add_splitvar_breaks_index_new(
      party_object = pt,
      plot_data = ggparty:::get_plot_data(pt),
      round_digits = 3
    )


    # 可視化
    ggparty(pt) +

      geom_edge() +
      geom_edge_label(aes(label = unlist(rounded_labels)), size = 3) +

      # ノード分割に用いる項目を追加
      geom_node_label(
        aes(colour = splitvar),
        line_list = list(
          aes(label = splitvar),
          aes(label = stringr::str_c("N = ", nodesize, sep = ""))
        ),
        line_gpar = list(
          list(size = 12, fontface = "bold"),
          list(size = 9, colour = "black")
        ),
        show.legend = F,
        ids = "inner"
      ) +

      # 終端ノードに箱ひげ図を指定
      geom_node_plot(
        gglist = list(
          geom_boxplot(aes(x = "", y = rank)),
          labs(
            x = NULL
          )
        ),
        nudge_x = -0.01,
        shared_axis_labels = T
      ) +

      # 終端ノードの件数を追加
      geom_node_label(
        line_list = list(
          aes(label = stringr::str_c("N = ", nodesize, sep = ""))
        ),
        line_gpar = list(
          list(size = 9, colour = "black")
        ),
        nudge_y = -0.4,
        ids = "terminal"
      )
  }

  • f3(重厚)が 0.256(下位25%)〜0.276(下位35%)に位置する銘柄群が最も人気の高いクラス(一番左)として識別されており、これは第 1 主成分 PC1 の f3 が負の方向に大きく特徴づけられている事とも整合性がとれている
  • f3(重厚)が 0.276(下位35%)以上であり、PC1 のスコアがそれなりに低い(-0.738以下)銘柄群が最も人気の低いクラス(一番右)として識別されている事もこれまで見てきた特徴と整合性がとれている
※上記の人気の高い/低いは人気 Top 50 の中での相対的なものである事に注意



5. まとめ

  • さけのわデータは利用規約が緩く使い勝手が良い
  • json ベースの Web API は R でも簡単に扱える
  • 人気銘柄の傾向をそれなりに捉えられて楽しい
  • 今回は使用していないフレーバータグも試してみたい

参考

2020年9月9日水曜日

ベルヌーイ施行における最尤推定と事後分布

複数回のベルヌーイ施行においてパラメータ $\theta = p$ を最尤推定および事後分布(からの点推定)で推定する際に
  • どの程度の差異が発生するのか
  • どのように差異が収斂していくのか
をシミュレーションで検証する。

この記事は以前に記載した ベルヌーイ分布と事後分布 の 「1.1.3 無情報事前分布における MAP と MLE の一致」 で簡潔に述べた内容を含み、より広い条件で検証するものである。


1. 導出された結果(再掲)


以前の記事で導出した結果のみを下記に記載する。
※詳細な定義および導出の過程は ベルヌーイ分布と事後分布 を参照

  • $p(\theta | x^{n})$: 事後分布
  • $x_{j}$: サンプルデータ
  • $n$: サンプル数
  • $a, b$: 事前分布 $Beta(a,b)$ のパラメータ

\begin{eqnarray} p(\theta | x^{n}) &=& Beta(a + \sum_{j=1}^{n} x_{j}, n + b - \sum_{j=1}^{n} x_{j}) \\ \hat{\theta}_{MLE} &=& \frac{1}{n} \sum_{j=1}^{n} x_{j} \\ \hat{\theta}_{MAP} &=& \frac{a + \sum x_{j} - 1}{n + a + b - 2} \\ \hat{\theta}_{EAP} &=& \frac{ a + \sum x_{j} }{ n + a + b } \end{eqnarray}

2. 事前分布


上記ではベルヌーイ施行における事前分布としてベータ分布を仮定しており、以下ではパラメータ $a, b$ を 2 パターン設定して検証を行う。
  1. (a, b) = (1, 1)
  2. (a, b) = (3, 6)
上記 1. は無情報事前分布であり下記の $Beta(1,1)$ に示すように定義域が $[0, 1]$ の一様分布となる。
library(tidyverse)

# Beta 分布のパラメータ
lst.beta <- list(
  a = 3,
  b = 6
)

tibble(
  p = seq(0, 1, 0.01),
  d_1_1 = dbeta(0, 1, 1),                  # Beta(1,1)
  d_a_b = dbeta(p, lst.beta$a, lst.beta$b) # Beta(a,b)
) %>%

  # to wide-form
  tidyr::pivot_longer(cols = -p, names_to = "param", values_to = "density") %>%

  # 可視化
  ggplot(aes(p, density)) +
    geom_line() +
    labs(y = NULL) +
    facet_grid(
      . ~ param,
      labeller = labeller(
        param = as_labeller(c(
          "d_1_1" = "B(1,1)",
          "d_a_b" = stringr::str_c("Beta(", lst.beta$a, ",", lst.beta$b, ")")
        ))
      )
    )



2. データの生成


真のパラメータ $p = 0.3$ を指定して 1,000 個の TRUE/FALSE からなる乱数 $x$ を生成する。
set.seed(1025)

n <- 1000
p <- 0.3

# データ生成
x <- rbernoulli(n, p)

# 先頭 10 件の確認
x[1:10]
# FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE FALSE

# x の平均値
# 真のパラメータ p=0.3 にかなり近い
mean(x)
# 0.29

上記 $x$ を用いて各累積施行ごとにパラメータ $p$ を推定していく。
※コメント内の番号 (1)〜(4) は 「1. 導出された結果(再掲)」 における数式の番号に紐づく
df.estimations <- tibble(
  # 延べ試行回数
  n = 1:length(x),
  
  # 成功: 1
  # 失敗: 0
  x = as.integer(x),

  # 延べ成功回数
  cumsum_x = cumsum(x),

  # 最尤推定による点推定値 - (2)
  MLE = cumsum_x / n,

  # MAP - (3)
  MAP_1_1 = (1 + cumsum_x - 1) / (n + 1 + 1 -2),                            # 事前分布: Beta(1,1)
  MAP_a_b = (lst.beta$a + cumsum_x - 1) / (n + lst.beta$a + lst.beta$b -2), # 事前分布: Beta(3,6)

  # EAP - (4)
  EAP_1_1 = (1 + cumsum_x) / (1 + 1 + n),                                   # 事前分布: Beta(1,1)
  EAP_a_b = (lst.beta$a + cumsum_x) / (lst.beta$a + lst.beta$b + n),        # 事前分布: Beta(3,6)

  # 事後分布の算出 - (1)
  post_dists = purrr::map2(n, cumsum_x, function(n, cumsum_x, a, b) {
    tibble::tibble(
      p = seq(0, 1, 0.01),
      post_d_1_1 = dbeta(p, 1 + cumsum_x, n + 1 - cumsum_x),                # 事前分布: Beta(1,1)
      post_d_a_b = dbeta(p, a + cumsum_x, n + b - cumsum_x)                 # 事前分布: Beta(3,6)
    )
  }, a = lst.beta$a, b = lst.beta$b)
)

df.estimations
      n             x       cumsum_x MLE MAP_1_1 MAP_a_b EAP_1_1 EAP_a_b post_dists
1 0 0 0 0 0.25 0.333 0.3 <tibble [101 × 3]>
2 0 0 0 0 0.222 0.25 0.273 <tibble [101 × 3]>
3 0 0 0 0 0.2 0.2 0.25 <tibble [101 × 3]>
4 0 0 0 0 0.182 0.167 0.231 <tibble [101 × 3]>
5 0 0 0 0 0.167 0.143 0.214 <tibble [101 × 3]>


上記 post_dists 項目には各行(累積施行)ごとに事後分布のデータが tibble 形式で格納されており、項目を指定しないと確認する事はできない。
下記は 1 回目の施行によりデータを 1 つ獲得した後の事後分布を 5 行だけ表示したものであり、post_d_xxx は確率密度を表す。

df.estimations[1,]$post_dists[[1]]
      p       post_d_1_1 post_d_a_b
0 2 0
0.01 1.98 0.0237
0.02 1.96 0.0893
0.03 1.94 0.189
0.04 1.92 0.316


3. 事後分布の確認


各施行ごとの事後分布を確認していく。

3.1 少ない試行回数での挙動


試行回数が少ない時の挙動を確認する。主に 1〜10 回をメインに観察する。
df.estimations%>%

  # 試行回数が特定の時点のみを対象
  dplyr::filter(n %in% c(1, 2, 3, 4, 5, 10, 25, 50, 100)) %>%

  # 事後分布の抽出
  tidyr::unnest(post_dists) %>%
  dplyr::select(n, MLE, p, post_d_1_1, post_d_a_b) %>%

  # => long-form
  tidyr::pivot_longer(cols = c(post_d_1_1, post_d_a_b), names_to = "type", values_to = "density") %>%

  # 可視化
  ggplot(aes(p, density)) +
    geom_line(aes(colour = type), alpha = 1, show.legend = F) +
    geom_vline(aes(xintercept = MLE), size = 1.0, linetype = "solid", colour = "gray") +
    labs(
      x = NULL,
      y = NULL
    ) +
    facet_grid(
      n ~ type,
      scale = "free_y",
      labeller = labeller(
        type = as_labeller(c(
          "post_d_1_1" = "Beta(1,1)",
          "post_d_a_b" = stringr::str_c("Beta(", lst.beta$a, ",", lst.beta$b, ")")
        ))
      )
    )

左列には $Beta(1,1)$ を事前分布(無情報事前分布)とする事後分布が、右列には $Beta(3,6)$ を事前分布とする事後分布が表示されている。
右側の番号は試行回数を表しており、下に進むほど試行回数が大きくなる。
各セルにおける横軸上の三角は MLE を表しており、各行ごとに左右の列で同じ位置に表示されている事に注意。


上記より把握できる事を下記にまとめる。
  • 試行回数が少ない場合に MLE はデータに過度に反応し、試行回数 1〜5 では p = 0 という極端な値を提示してしまっている
  • 左列である Beta(1,1) (無情報事前分布)側におけるピークの位置(MAP) と MLE は一致している(理論上完全に一致する)
  • 最初の 5 回で FALSE が続いている事を受けて Beta(1,1) 側はデータに過度に反応した事後分布を生成してしまっている
  • 最初の 5 回の FALSE で Beta(1,1) 側は試行回数が増える程に p = 0 の確信が高まっている様子が伺える
  • 試行回数が少ない時点では Beta(1,1) 側と Beta(3,6) 側では差異が大きい
  • 試行回数が大きくなるにつれ Beta(1,1) 側と Beta(3,6) 側の差異は縮まる傾向にある
  • Beta(1,1) 側および Beta(3,6) 側の両方で、試行回数が増えるにつれてピーク付近でのばらつきが減少している(横が狭くなる)事およびピークの位置が高くなっている事から段々と推定の確信が高まっている様子が伺える

3.2 事前分布による差異


$Beta(1,1)$ と $Beta(3,6)$ の異なる事前分布における挙動の差異を確認する。
df.estimations %>%

  dplyr::select(n, MAP_1_1, MAP_a_b, EAP_1_1, EAP_a_b) %>%

  # => long-form
  tidyr::pivot_longer(cols = c(MAP_1_1, MAP_a_b, EAP_1_1, EAP_a_b), names_pattern = "(MAP|EAP)_(.+_.+)", names_to = c("type", "param"), values_to = "prob") %>%

  # 可視化
  ggplot(aes(n, prob)) +
    geom_line(aes(color = param), alpha = 1/2) +
    geom_point(aes(colour = param), size = 0.75, alpha = 1/5) +

    # 真のパラメータ 0.3 を表示
    geom_hline(yintercept = 0.3, alpha = 1/3) +

    scale_y_continuous(limits = c(0, 0.5)) +
    scale_colour_hue(
      name = "Prior Dist",
      labels = c(
        "1_1" = "Beta(1,1)",
        "a_b" = stringr::str_c("Beta(", lst.beta$a, ",", lst.beta$b, ")")
      )
    ) +
    labs(
      x = "Trials",
      y = "Estimated Parameter"
    ) +
    facet_grid(type ~ .)

横軸に試行回数、縦軸にパラメータ p の推定値を表示する。
水平に伸びた直線はパラメータ p の真の値 0.3 を示している。
上下で推定値が異なる事に注意。


  • Beta(1,1) (無情報事前分布)を用いた推定値の方が上下の振れ幅が大きい傾向にある
  • MAP/EAP ともに試行回数が少ない範囲では事前分布の違いにより挙動が大きく異る
  • 試行回数が大きくなると事前分布による差異はほぼ見られくなる

3.3 点推定値による差異


MAP と EAP における挙動の差異を確認する。
df.estimations %>%

  dplyr::select(n, MAP_1_1, MAP_a_b, EAP_1_1, EAP_a_b) %>%

  # => long-form
  tidyr::pivot_longer(cols = c(MAP_1_1, MAP_a_b, EAP_1_1, EAP_a_b), names_pattern = "(MAP|EAP)_(.+_.+)", names_to = c("type", "param"), values_to = "prob") %>%

  # 可視化
  ggplot(aes(n, prob)) +
    geom_line(aes(color = type), alpha = 1/2) +
    geom_point(aes(colour = type), size = 0.75, alpha = 1/5) +

    # 真のパラメータ 0.3 を表示
    geom_hline(yintercept = 0.3, alpha = 1/3) +

    scale_y_continuous(limits = c(0, 0.5)) +
    labs(
      colour = NULL,
      x = "Trials",
      y = "Estimated Parameter"
    ) +
    facet_grid(
      param ~ .,
      labeller = labeller(
        param = as_labeller(c(
          "1_1" = "Beta(1,1)",
          "a_b" = stringr::str_c("Beta(", lst.beta$a, ",", lst.beta$b, ")")
        ))
      )
    )

横軸に試行回数、縦軸にパラメータ p の推定値を表示する。
水平に伸びた直線はパラメータ p の真の値 0.3 を示している。
上下で事前分布が異なる事に注意。


  • $Beta(1,1)$/$Beta(3,6)$ ともに試行回数が少ない範囲では推定量ごとの挙動に差異が見られる
  • 試行回数が大きくなると推定量による差異はほぼ見られくなる
  • EAP(赤色) > MAP(青色) の傾向がある(3.2.1 で後述)

3.2.1 EAP > MAP について


特定の状況においてほぼ EAP > MAP が成立する事を示す。
以下では EAP と MAP の差分を考える事で上記を検証していく。

\begin{eqnarray} \hat{\theta}_{EAP} - \hat{\theta}_{MAP} = \frac{ a + \sum x_{j} }{ n + a + b } - \frac{a + \sum x_{j} - 1}{n + a + b - 2}. \nonumber \end{eqnarray} 上記は分母が基本的に正となるので、計算を楽にするため以下では分子だけを考える。 \begin{eqnarray} 分子 &=& \{ (a + \sum x_{j})(n + a + b - 2) - (a + \sum x_{j} - 1)(n + a + b) \} \nonumber \\ &=& \Bigl[ \bigl\{ (a + \sum x_{j} - 1)(n + a + b - 2) + (n + a + b - 2) \bigr\} - (a + \sum x_{j} - 1)(n + a + b) \Bigr] \nonumber \\ &=& \{ (a + \sum x_{j} - 1)(- 2) + (n + a + b - 2) \} \nonumber \\ &=& \{ (b - a) + 2n (\frac{1}{2} - \frac{1}{n} \sum x_{j}) \} \nonumber \end{eqnarray} 上記より、パラメータ $p$ の不偏推定量である標本平均 $\frac{1}{n} \sum x_{j}$ が $\frac{1}{n} \sum x_{j} < \frac{1}{2}$ を満たせばほぼ $\hat{\theta}_{EAP}>\hat{\theta}_{MAP}$ が成立する。
今回はベルヌーイ分布の真のパラメータが $p = 0.3$ であり $\frac{1}{2}$ より小さい事からほぼ \begin{eqnarray} \therefore \hat{\theta}_{EAP} > \hat{\theta}_{MAP} \nonumber \end{eqnarray} が成立するものと考えられる(適当)。

3.3 試行回数による分布の変化


試行回数の増加に伴う分布の変化をマクロな視点で確認する。
df.estimations %>%

  dplyr::select(n, post_dists) %>%

  # 事後分布の抽出
  tidyr::unnest(post_dists) %>%

  # => long-form
  tidyr::pivot_longer(cols = c(post_d_1_1, post_d_a_b), names_to = "param", values_to = "density") %>%

  # 可視化
  ggplot(aes(n, p)) +
    geom_tile(aes(fill = density), show.legend = F) +
    scale_fill_gradient(low = "#132B43", high = "tomato") +
    labs(
      x = "Trials",
      y = "p"
    ) +
    facet_grid(
      param ~ .,
      labeller = labeller(
        param = as_labeller(c(
          "post_d_1_1" = "Beta(1,1)",
          "post_d_a_b" = stringr::str_c("Beta(", lst.beta$a, ",", lst.beta$b, ")")
        ))
      )
    ) +

    # 不要な余白を除去
    coord_cartesian(expand = F)

横軸に試行回数、縦軸にパラメータ p を表示し、各点における事後分布の確率密度を色で表現する。
黒に近いほど確率密度が低く、赤く明るいほど確率密度が高い事を示す。
上下で事前分布が異なる事に注意。


  • 全体を外観すると事前分布が異なる事による差異はほぼ見られない
  • 左端の方では明確な赤色が見られず全体としてぼやけた感じとなっている
  • 右側に進むにつれて赤と黒が明確に分離していく様子が見て取れる
  • 試行回数の増加と共にパラメータ推定の確信度が上がっていると考えられる
  • 試行回数が大きくなると分布に殆ど変化が見られなくなり、多少のデータ傾向の揺れに振り回されなくなってくる(左端の挙動と対称的)

上記から、ある程度の試行回数が存在しない限りパラメータ p の推定はあまり意味を持たない可能性があると思われる。


まとめ

  • サンプルが少ない場合に、MLE はデータに過度に反応し、データによってはかなり極端な値を提示してしまう
  • サンプルが少ない場合に、推定された事後分布はデータに過度に反応する可能性がある
  • サンプルが少ない場合に、推定された事後分布は事前分布の違いに大きく影響を受ける可能性がある
  • サンプルが少ない場合に、点推定量(MLE/MAP/EAP)による違いが現れる可能性がある
  • サンプルが増加するにつれ、推定される事後分布が多少のデータ傾向の揺れに影響を受けづらくなり安定する
  • サンプルが増加するにつれ、推定される事後分布は事前分布の違いにあまり影響を受けなくなる
  • サンプルが増加するにつれ、点推定量の違いによる差異はほぼ見られなくなる
  • ある程度のサンプルが存在しないとパラメータの推定はかなり厳しい

2020年9月3日木曜日

JuMP でビンパッキング問題おためし

Julia の数理最適化パッケージ JuMP を用いて混合整数線型計画問題を試してみる。
JuMP で通常の線形計画問題を解くのは JuMP で線型計画問題おためし を参照。


1. 環境の準備


今回用いるパッケージは下記。
  • JuMP
  • CBC (MILP:混合整数計画問題 のソルバ)
※環境構築に関しては JuMP で線型計画問題おためし を参照
※各 Solver に関しては Getting Solvers を参照


2. ビンパッキング問題


ビンパッキング問題の詳細は Wikipedia を参照してもらうとして、 今回は下記の状況を扱う。

  • それぞれ重さの異なる荷物を大きめのビンに詰める事を考える
  • 使用するビンの数を最も少なくする格納方法を考えたい
  • 荷物 10 個
    • 重さ 2 kg: 4 個
    • 重さ 3 kg: 3 個
    • 重さ 4 kg: 2 個
    • 重さ 5 kg: 1 個
  • ビンごとの最大容量 5 kg
※ちなみに 「ビン」 は 「瓶」 ではなく 「Bin(容器)」 の事


3. 定式化


荷物の数を n とすると必要なビンの数は最大でも n 個となるので、使われないビンの存在も許容してビンの数を(多めに) n と設定する。
このとき各ビンに荷物が 1 つでも格納されたかどうかを判定する変数 y[1〜n] を、各成分に 0 か 1 が格納されるベクトルで表現する。

\begin{eqnarray} \forall j \in \{ 1, 2, \dots, n \}, \ \ y_{j} \in \{ 0, 1 \} \nonumber . \end{eqnarray} 次に行列 $X = \{ x_{i, j} \}$ を下記で定義する。
  • 行(i): 各荷物(1〜n)に対応
  • 列(j): 各ビン(1〜n)に対応
  • 各 $x_{i, j}$ には 0 か 1 が格納される

最後にベクトル $s$ 及びスカラ $B$ を下記で定義する。
  • $s$: 各荷物の重さ
    • 今回の場合は $s = [ 2, 2, 2, 2, 3, 3, 3, 4, 4, 5 ]$
  • $B$: 各ビンごとの最大容量
    • 今回の場合は $B = 5$

3.1 目的関数


出来るだけビンの数を少なくしたいので、各 $y_{j}$ の和を最小とする事を考える。 \begin{eqnarray} minimize \sum_{j} y_{j} \end{eqnarray}

3.2 制約条件(1)


まず、各荷物 1〜n がそれぞれ 1 つのビンに必ず格納される事を表現する。
\begin{eqnarray} \forall i, \ \ \sum_{j} x_{i, j} = 1 \nonumber \end{eqnarray}
この条件は n 個の成分が全て 1 であるベクトル $e (\forall i, \ e_{i} = 1)$ を用いて下記のように書き換え可能。 \begin{eqnarray} X \cdot e = 1 \end{eqnarray}

3.3 制約条件(2)


次に、各ビンに格納された荷物の重さの合計が $B$ を以下となる事を表現する。 \begin{eqnarray} \forall j, \ \ \sum_{i} s_{i} \cdot x_{i, j} \leqq B \cdot y_{j} \nonumber \end{eqnarray}
この条件は行列の転置を用いて下記のように書き換え可能。 \begin{eqnarray} {}^t\!X \cdot s \leqq B \cdot y \end{eqnarray}

3.4 制約条件(3)


最後に、行列 $X$ およびベクトル $y$ の各成分に関する条件を表現する。
\begin{eqnarray} \forall i, j, \ \ x_{i, j} \in \{ 0, 1 \} \\ \forall j, \ \ y_{j} \in \{ 0, 1 \} \end{eqnarray}

3.5 まとめ


上記 (1)〜(5) を用いて下記のようにまとめられる。
\begin{eqnarray} minimize & & \sum_{j} y_{j} \nonumber \\ subject \ to & & X \cdot e = 1 \nonumber \\ & & {}^t\!X \cdot s \leqq B \cdot y \nonumber \\ & & \forall i, j, \ \ x_{i, j} \in \{ 0, 1 \} \nonumber \\ & & \forall j, \ \ y_{j} \in \{ 0, 1 \} \nonumber \end{eqnarray}

4. JuMP


JuMP を用いて最適解を求める。
MILP(混合整数線型問題) を扱う事のできる CBC をソルバーに用いる。
using JuMP
using Cbc
using LinearAlgebra

# 荷物の数 => 最大のビンの数
n = 10

# ビンごとの最大容量
B = 5

# 荷物ごとの重さ
s = [2, 2, 2, 2, 3, 3, 3, 4, 4, 5]


# CBC を用いたモデルを作成
model = Model(Cbc.Optimizer)

# 変数を追加
# X: 荷物とビンの関係
# y: ビン使用の有無
# binary = true で 2 値(0 or 1)変数である事を指定
@variable(model, X[i=1:n, j=1:n], binary = true) # (4)
@variable(model, y[1:n], binary = true)          # (5)

# 目的関数を最小化
@objective(model, Min, sum(y))                   # (1)

# 制約条件を追加
@constraint(model, X * ones(n) .== 1)            # (2)
@constraint(model, transpose(X) * s .<= B .* y)  # (3)

# 最適化の実行
optimize!(model)


### 実行結果 ###

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: May 23 2020 

command line - Cbc_C_Interface -solve -quit (default strategy 1)
Continuous objective value is 6 - 0.00 seconds
〜中略〜

Result - Optimal solution found

Objective value:                7.00000000
Enumerated nodes:               0
Total iterations:               83
Time (CPU seconds):             0.37
Time (Wallclock seconds):       0.12

Total time (CPU seconds):       0.37   (Wallclock seconds):       0.12

4.1 実行結果


正常に終了し、各荷物の各ビンへの格納が正しく行われている事を確認できる。
println("TerminationStatus: ", termination_status(model))
# => TerminationStatus: OPTIMAL

println("PrimalStatus: ", primal_status(model))
# => PrimalStatus: FEASIBLE_POINT

println("ObjectiveValue: ", objective_value(model))
# => ObjectiveValue: 7.0

# 荷物の格納されたビンの情報のみを表示
for i = 1:n, j = 1:n
    x = value(X[i,j])
    if x > 0
        println("X[", i, "," , j, "]: ", x)
    end
end
X[1,2]: 1.0
X[2,6]: 1.0
X[3,7]: 1.0
X[4,9]: 1.0
X[5,7]: 1.0
X[6,9]: 1.0
X[7,2]: 1.0
X[8,10]: 1.0
X[9,8]: 1.0
X[10,4]: 1.0

ビンと荷物の対応
  • ビン 2:
    • 荷物1(2kg)
    • 荷物7(3kg)
  • ビン 4:
    • 荷物10(5kg)
  • ビン 6:
    • 荷物2(2kg)
  • ビン 7:
    • 荷物3(2kg)
    • 荷物5(3kg)
  • ビン 8:
    • 荷物9(4kg)
  • ビン 9:
    • 荷物4(2kg)
    • 荷物6(3kg)
  • ビン 10:
    • 荷物8(4kg)


参考