ほくそ笑む

R言語と統計解析について

非心ベータ分布を Stan で推定する (1)

1. はじめに

ベータ分布は範囲が固定されている連続データの分布を柔軟に表現できるため便利である。 非心ベータ分布 (Noncentral Beta Distribution) はベータ分布を一般化した確率分布であり柔軟性はさらに上がる。

R ではベータ分布に従う乱数を生成する rbeta() 関数に ncp 引数1を指定することで非心ベータ分布に従う乱数を生成できる。

x <- rbeta(2000, 10, 5, ncp = 3)

hist(x, probability = TRUE)
curve(dbeta(x, 10, 5), add = TRUE, col = "red")
curve(dbeta(x, 10, 5, ncp = 3), add = TRUE, col = "blue")

f:id:hoxo_m:20180823223727p:plain

赤線がベータ分布、青線が非心ベータ分布であり、ヒストグラムへの当てはまりは青線の方が良い。

非心ベータ分布はベータ分布の中心をずらした分布になる。 右にずらす場合は Type I、左にずらす場合は Type II と呼ばれる。 rbeta() は Type I にしか対応していないが、Type I と Type II は単に左右反転しただけなので 1 - x を取ることで Type II を表せる。

本記事では Stan を使って非心ベータ分布のパラメータを推定する方法について述べる。

2. 対数尤度関数

Stan には非心ベータ分布を推定するための組み込み関数はない。 したがって、target += 記法を用いてモデルを記述する。 そのために、非心ベータ分布の対数尤度関数を求める必要がある。

定義より、非心ベータ分布の密度関数は次である。

\begin{align} \mathrm{NCBeta}(x \mid \alpha, \beta, \lambda) = \sum_{j = 0}^\infty \mathrm{Pois}(j \mid \lambda) \mathrm{Beta}(x \mid \alpha + j, \beta) \end{align}

ただし、Pois はポアソン分布の密度関数、Beta はベータ分布の密度関数である。

ここから尤度関数を求めたいが、 j についての無限和の存在が問題になる。 ここで、ベータ分布の密度関数を  j に対するオーダーで書き直すと

\begin{align} \mathrm{Beta}(x \mid \alpha + j, \beta) &= \frac{x^{\alpha+j-1}(1-x)^{\beta-1}}{\mathrm{B}(\alpha + j, \beta)} \\ &= \frac{O(x^j)}{O\left(\frac{(\alpha+j-1)!(\beta - 1)!}{(\alpha + j + \beta - 1)!}\right)} \\ &= \frac{O(x^j)}{O\left(\frac{1}{j^\beta}\right)} \\ &= O(x^j j^\beta) \end{align}

となる。ただし B はベータ関数であり、そのオーダーとして  \alpha \beta が整数の場合を考えた。

一方、ポアソン分布の密度関数の  j に対するオーダーは次のようになる。

\begin{align} \mathrm{Pois}(j \mid \lambda) &= \frac{\left(\lambda/2\right)^j}{j!} e^{-\lambda} \\ &= \frac{O(\lambda^j)}{O(j!)} \end{align}

すなわち、ベータ分布は  x = 1 のときでも高々  O(j^\beta) でしか増加しないのに対し、ポアソン分布は  O(j!) の速度で減少する。

したがって、非心ベータ分布の密度関数の  j に関して次のことが言える。

  •  j が十分大きいとき項は 0 とみなせる。
  • ポアソン分布の重みが支配的である。

すなわち、0 とみなせる  j についてはポアソン分布の重みのみを考えれば良い。

ポアソン分布の相対的な重みを比較するには次のようにする。

lambda <- 7
j <- rpois(10000, lambda / 2)
round(table(j) / max(table(j)), 2)
j
   0    1    2    3    4    5    6    7    8    9   10   11   12   13 
0.14 0.49 0.81 1.00 0.82 0.58 0.37 0.18 0.08 0.03 0.01 0.00 0.00 0.00 

これより、 \lambda \leq 7 に限れば、 j は 0 から 10 までを考えれば十分である。  \lambda の上限を上げると  j の範囲も広がり、そのぶん計算にかかる時間は増加することに注意。

以上によりベータ分布の対数尤度関数は次のように表される。

\begin{align} \log\mathrm{NCBeta}(x \mid \alpha, \beta, \lambda) &= \log\left( \sum_{j = 0}^{10} \mathrm{Pois}(j \mid \lambda) \mathrm{Beta}(x \mid \alpha + j, \beta) \right) \\ &= \log \sum_{j=0}^{10} \exp\left( \log \mathrm{Pois}(j \mid \lambda) + \log \mathrm{Beta}(x \mid \alpha + j, \beta) \right) \end{align}

これは log_sum_exp の形になっていることに注意せよ。

3. Stan による実装

非心ベータ分布の対数尤度関数がわかったので、これを Stan で実装すると次のようになる。

data {
  int<lower=1> N;
  real<lower=0, upper=1> values[N];
}
parameters {
  real<lower=1> alpha;
  real<lower=1> beta;
  real<lower=0, upper=7> lambda;
}
model {
  real logliks[11];
  for (i in 1:N) {
    for (j in 0:10) {
      logliks[j + 1] = poisson_lpmf(j | lambda / 2);
      logliks[j + 1] += beta_lpdf(values[i] | alpha + j, beta);
    }
    target += log_sum_exp(logliks);
  }
}

これを ncbeta01.stan という名前でファイルに保存する。

推定を実行する R コードは以下である。

set.seed(314)

alpha <- 10
beta <- 5
ncp <- 3
N <- 200

x <- rbeta(N, alpha, beta, ncp)

standata <- list(N = N, values = x) # Type I 
# standata <- list(N = N, values = 1 - x) # Type II

library(rstan)
# options(mc.cores = 3)
model <- stan_model("ncbeta01.stan")
fit <- sampling(model, data = standata, chains = 3)

print(fit)
         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha    8.93    0.03 1.15   6.70   8.18   8.91   9.70  11.29  1160    1
beta     4.57    0.01 0.42   3.80   4.28   4.56   4.86   5.43  1491    1
lambda   3.17    0.06 1.97   0.17   1.45   3.04   4.81   6.73  1255    1
lp__   150.20    0.04 1.25 147.12 149.60 150.52 151.13 151.65   928    1

真の値  \alpha = 10,  \beta = 5,  \lambda = 3 がうまく推定できたようである。

まとめ

非心ベータ分布を Stan で推定する方法を述べた。 上記のコードだと私の PC では 1 chain 80秒程度かかる。

次回はこれを高速化して 1 chain 10秒程度にする方法を書こうと思う。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)


  1. npc は Non-Centrality Parameter の略

メモ: Local Explain における Fidelity の評価

機械学習モデルの解釈において、Local Explain が流行っている。 Global Explain の場合、解釈手法の評価は Fidelity によってなされる。 Fidelity (忠実度) は、解釈したいモデル f(x) を説明 g(x) がどれくらい忠実に模倣しているかを測る指標である。 通常は f(x) を正解として g(x) の精度 (AUC や F-measure などの普通の指標) として定義される。 しかし、Local Explain の場合、Fidelity の評価は難しそうだ。

LIME の場合、5.2 で次のように評価している。

  1. 特徴量の重要度のわかるモデル (スパースロジスティック回帰と決定木) に対して、訓練データを学習する。
  2. 重要度の高い上位 10個をゴールドフィーチャとする。
  3. LIME で説明を作り、説明で使われる特徴量がゴールドフィーチャに含まれる割合を出す。

LIME では平均して 90% 以上がゴールドフィーチャに含まれたようだが、なんというか、コレジャナイ感が強い。 Local Explain (というか Outcome Explain) における Fidelity の評価方法については一考の余地がありそうだ。

useR! 2018 メモ (3) 後半

SAR: a practical, rating-free hybrid recommender for large data

Microsoft のレコメンド手法 SAR (Smart Adaptive Recommendation) について。 Azure の API の紹介だが、R 実装もあるとのこと。

Jumpstart Machine Learning with Pre-Trained Models

こちらも Microsoft の人の発表で、訓練済み機械学習モデルを共有する MLHub の紹介。 https://mlhub.ai

FastR: an alternative R language implementation

Oracle の人で高速な R 実装 FastR の紹介。 Renjin が有名だが、それの対抗馬。GraalVM 上で動く。

Estimating individual Customer Lifetime Values with R: The CLVTools Package

顧客生涯価値 (Customer Life Value; CLV) を確率的な顧客離脱モデルを使って推定するパッケージ CLVTools の紹介。このパッケージまだ公開されていないらしい。

Modelling Field Operation Capacity using Generalised Additive Model and Random Forest

イギリスの電力会社かなんかで顧客の修理依頼の数を予測してスタッフを十分確保しておくという話。

iml: A new Package for Model-Agnostic Interpretable Machine Learning

機械学習の解釈可能性についての論文をいくつか実装したパッケージ。 LIME の独自実装を含む。

Teaching R to New Users: From tapply to Tidyverse

しゃべりのうまいおじさん

Disciplined Convex Optimization with CVXR

凸最適化用のパッケージ CVXR の紹介。 凸最適化問題はプログラムで記述するのが面倒だけどこれだと自然な数式で表現できる。

Robustness criterion for the derivative kriging-based optimization

ロバストなクリギング手法。より安定な点に最適化するようにロバスト基準を加える。

Augmented Lagrangian for constrained optimizations in Empirical Likelihood estimations

Empirical Likelihood 法というので制約条件つける方法。 拡張ラグランジュで最適化するらしい。

Data Preprocessing using Recipes

その他 useR! 2018 参加報告

engineering.linecorp.com gingi99.hatenablog.com

施策効果測定におけるメタアナリシスの応用

はじめに

マーケティング施策を行うときに、その施策効果を測定するために、コントロールグループ(施策を適用しないユーザ)を作る場合がある。 例えば、販促メールを送るという施策を行うときに、一部のユーザには送らないようにする。 仮にメールを送らなかったユーザの平均売上が1000円であり、メールを送ったユーザは平均1100円だとすると、その差である100円が1人あたりの施策効果となる。

しかし、施策を適用しないユーザの数を増やすと、全体の売上効果はそれだけ減少してしまう。 つまり、全体のユーザ数を1万人とすると、全員にメールを送れば100万円の売上効果があるが、半分の5千人をコントロールとすると売上効果も半分の50万円となる。 したがって、コントロールグループに割り当てる人数はなるべく小さくしたいという要求がある。

ただし、コントロールグループの人数を少なくすると、効果測定の精度が落ちるという問題がある。 今、1万人のユーザがいるECサイトで、メール施策を行うことを考えよう。 このとき、コントロールグループは全体の 1% である 100人しか取れないとする。 ユーザの作り出す売上がポアソン分布に従うとして、

  • メールを送らないユーザは平均 1000 円の売上
  • メールを送ると平均 1100 の売上

すなわちメール施策の真の効果は 100円としてシミュレーションしてみる。

set.seed(17744)
# 施策未実施ユーザ 100人の個々人の売上
control_group   <- rpois( 100, 1.0) * 1000
# 施策実施ユーザの 9900人の個々人の売上
treatment_group <- rpois(9900, 1.1) * 1000

# それぞれのグループに対して10人の売上を表示
head(control_group, 10)
#> [1] 1000 2000 2000 2000 1000 1000    0    0    0 4000
head(treatment_group, 10)
#> [1] 4000 1000 3000 1000    0 1000 2000 1000 1000 1000

こうして得られた各ユーザの売上を使って、施策効果を推定する。 各グループの平均売上の差が施策効果である。

effect <- mean(treatment_group) - mean(control_group)
effect
#> [1] 47.3

施策効果は 47 円と算出され、真の効果である 100円から大きくずれてしまった。

また、施策に本当に効果があったかは統計的仮説検定を用いて判断することが多い。 ここでは Welch t-test を使ってみよう。

t.test(treatment_group, control_group)
 Welch Two Sample t-test

data:  treatment_group and control_group
t = 0.41425, df = 100.68, p-value = 0.6796
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -179.1100  273.6555
sample estimates:
mean of x mean of y 
 1087.273  1040.000 

検定の結果より、p-value = 0.68 なのでこの施策には有意な効果はないと判断されてしまう。

つまり、コントロールグループの人数が少ないことが原因で

  • 施策効果の推定値が真の値から大きくずれてしまう
  • 検定を行なっても有意な効果がないと判断されてしまう

という問題が起こる。

ただし、このようなマーケティング施策は繰り返し行われることが多い。 例えば、販促メールを毎週水曜日に送るなどである。 この場合、メタアナリシスを用いて上記の問題を解決できる。

今、上記のメール施策を毎週水曜日に5回繰り返したとする。 このときのメタアナリシスの結果を下図に示す。

f:id:hoxo_m:20180809192359p:plain

各週のメール施策において算出された効果は Effect 列に記載されている。 これはばらつきが大きく、マイナスになることもあることがわかる(2018-08-08)。 右側のグラフは、四角が Effect であり、そこから左右に伸びる線が信頼区間を表している。縦線は 0円であり、2018-08-22を除いて各施策は 0 をまたいでいるため、検定結果は有意でない。

メタアナリシスによってこれらの結果を統合することで、推定精度を向上させることができる。 メタアナリシスによる推定結果は上図の Summary に記されている。 これによると、効果の推定値は 114 円であり、真の効果 100円に近い。 また、信頼区間は 26円から 202円となり、0円をまたがないので検定は有意になる。

本記事ではRでこれを行う方法を説明する。

メタアナリシスの実践

次のように、5回の施策に対して、その効果と信頼区間を求める。

# 8/1(水) から 8/29(水) までの水曜日の日付を取得
begin <- as.Date("2018-08-01")
end   <- as.Date("2018-08-29")
dates <- seq(begin, end, by = "weeks")

result <- data.frame()
for (i in seq_along(dates)) {
  date <- dates[i]
  
  # 売上のシミュレーション
  set.seed(as.integer(date))
  control_group   <- rpois( 100, 1.0) * 1000
  treatment_group <- rpois(9900, 1.1) * 1000
  
  # 売上効果の算出
  effect <- mean(treatment_group) - mean(control_group)
  
  # Welch t-test による信頼区間の推定
  t <- t.test(treatment_group, control_group, conf.level = 0.95)
  lower <- t$conf.int[1]
  upper <- t$conf.int[2]

  # # ブートストラップ法による信頼区間の推定
  # library(simpleboot)
  # boot <- two.boot(treatment_group, control_group, mean, R = 1000)
  # ci <- boot.ci(boot, conf = 0.95, type = "perc")
  # lower <- ci$perc[4]
  # upper <- ci$perc[5]

  res <- data.frame(date, effect, lower, upper)
  
  result <- result %>% bind_rows(res)
}

print(result)
        date     effect     lower    upper
1 2018-08-01  47.272727 -171.8910 268.2418
2 2018-08-08  -4.545455 -224.9838 181.0983
3 2018-08-15  18.181818 -171.9945 194.4213
4 2018-08-22 306.969697  150.0565 465.5196
5 2018-08-29  92.323232 -129.1142 297.0604

それぞれの施策実施日について、効果と信頼区間が算出された。 ここでは信頼区間の推定にt検定を用いたが、売上が正規分布に従わないことが不安ならブートストラップ法を使うと良い。 ただし以下では信頼区間は平均について対称であることを仮定している。

この結果に対して『Rで楽しむ統計』(p.108) を参考にメタアナリシスを行う。

Rで楽しむ統計 (Wonderful R 1)

Rで楽しむ統計 (Wonderful R 1)

library(dplyr)

summary <- result %>%
  mutate(weight = 1/(upper - lower)^2, 
         effect = effect * weight / sum(weight)) %>%
  summarise(effect = sum(effect), sd = sqrt(1/sum(weight))/2) %>%
  mutate(lower = effect - sd, upper = effect + sd) %>%
  select(-sd) %>%
  mutate(date = "Summary") %>%
  select(date, everything())

print(summary)
     date   effect    lower    upper
1 Summary 114.2419 26.44446 202.0394

それぞれの施策に対して信頼区間の幅を2乗して逆数を取ったものを重みとする。 それぞれの施策効果に対して重み付けを行い、全てを足し合わせたものが統合された施策効果となる。 また、信頼区間は重みの総和の逆数の平方根によって算出される。 メタアナリシスによって算出された施策効果は 114円であり、真の効果100円と近い値が算出された。 また、信頼区間は 26円から202円であり、0円をまたがないため、統計的に有意な結果であることが言える。 ここではメタアナリシスのモデルとして固定効果モデルを用いたが、一般的には施策効果にはぶれが生じるため、それを考慮したい場合はランダム効果モデルを用いる。

メタアナリシスの可視化にはフォレストプロットというものが使われる。 これは次のように作成できる。

library(forestplot)

nrow <- nrow(result) # 施策の回数

# 施策結果とメタアナリシスの結果を統合
df <- result %>% 
  mutate(date = as.character(date)) %>% 
  bind_rows(summary)

# フォレストプロットの左側に表示される表データの作成(日付と効果)
labeltext <- list(c("Date", as.character(df$date)),
                  c("Effect", round(df$effect)))

# フォレストプロットの描画
forestplot(labeltext, mean = c(NA, df$effect), 
           lower = c(NA, df$lower), upper = c(NA, df$upper),
           is.summary = c(TRUE, rep(FALSE, nrow), TRUE),
           hrzl_lines = gpar(col="#444444"))

f:id:hoxo_m:20180809192359p:plain

それぞれの施策における効果と信頼区間、そして統合された結果が可視化された。

まとめ

マーケティング施策効果の測定において、

  • 施策を実施しないコントロールグループが作られる
  • コントロールグループの人数が少なく効果測定の精度が低い
  • 同じ施策が繰り返される

という場合において、メタアナリシスを用いて推定精度を向上させる方法を説明した。

参考文献

Rで楽しむ統計 (Wonderful R 1)

Rで楽しむ統計 (Wonderful R 1)

犬は「さよなら」を言わない

インタラクティブな分析とプログラミングの間には常に葛藤がある

インタラクティブな分析とプログラミングの間には常に葛藤がある。 インタラクティブに作業している場合には、R に期待通りの処理をしてもらい、R が間違った場合は自分で直ちに修正を行いたい。 一方、プログラミングの場合、微妙あるいは不可解な問題が生じたら、すぐにエラーを返す関数が望ましい。 関数を作成する場合も、この葛藤を忘れないようにしよう。 インタラクティブなデータ分析のための関数を書いている場合は、分析者の期待を推測し、多少想定外の処理が行われても自動的に補正するようにしよう。 これに対してプログラミングとしての関数の場合は、より厳格になり、関数を呼び出したユーザが期待する処理を推測して対応するのは間違いである。

Hadley Wickham『R言語徹底解説』(p.148)

R言語徹底解説

R言語徹底解説