ほくそ笑む

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

Stan でパラメータに大小関係の制約をつける

EC サイトなどでは、1 ページの読み込み時間の平均値が短いユーザほどコンバージョン率が高くなるという話があります。
この原因については、サクサク読み込みが終わる速いサイトの方がユーザは気前よく購入する、というわけではなく、単に速い回線を持つほどお金持ちなので購入しやすいという身も蓋もない話もあります。
今回はこの原因についてではなく、制約条件が付いた時のパラメータ推定の話をします。
とりあえず擬似データを作成しましょう。

true_cvr <- c(0.209, 0.126, 0.096, 0.093, 0.086, 0.077, 0.067, 0.057)

load_time <- c("0-1", "1-3", "3-7", "7-13", "13-21", "21-35", "35-60", "60+")
session <- c(1000, 6000, 4000, 1500, 700, 500, 200, 150)
set.seed(71)
cv <- unlist(Map(function(n, p) rbinom(1, n, p), session, true_cvr))

data.frame(load_time, cv, session)
  load_time  cv session
1       0-1 208    1000
2       1-3 769    6000
3       3-7 366    4000
4      7-13 142    1500
5     13-21  54     700
6     21-35  35     500
7     35-60  17     200
8       60+  11     150

ロード時間ごとに真のコンバージョン率(true_cvr)を設定して、コンバージョン数(cv)を二項分布からの乱数(rbinom)として生成しています。

このときのコンバージョン率を Stan で推定してみましょう。

// load_time.stan
data {
  int<lower=0> N;
  int<lower=0> cv[N];
  int<lower=0> session[N];
}
parameters {
  real<lower=0, upper=1> cvr[N];
}
model {
  cv ~ binomial(session, cvr);
}
library(rstan)
model <- stan_model("load_time.stan")
standata <- list(N = length(cv), cv = cv, session = session)
result <- sampling(model, data = standata, chains = 3)

この結果を見やすく可視化してみます。

library(tidyr)
df <- rstan::extract(result, pars="cvr")$cvr %>%
  data.frame %>%
  setNames(load_time) %>%
  gather(load_time, cvr)

cvr_hat <- colMeans(rstan::extract(result)$cvr)
point_df <- data.frame(load_time, true_cvr, cvr_hat) %>%
  gather(load_time, cvr) %>%
  setNames(c("load_time", "cvr_type", "cvr"))

library(ggplot2)
ggplot(df, aes(x = load_time, y = cvr)) + 
  geom_violin() + ylim(0, NA) +
  geom_point(data = point_df, aes(x=load_time, y=cvr, color=cvr_type))

パラメータの推定値(cvr_hat)と真値(true_cvr)には開きがあります。

制約条件の追加

ここからが本題です。
今、「読み込み時間が短いユーザほどコンバージョン率が高くなる」という条件が分かっています。
この条件を考慮してパラメータを推定することはできないでしょうか?
Stan では、パラメータを ordered 型で定義することによって、パラメータに順序関係の制約を加えることができます。
このような制約条件を付けた場合の Stan コードは次のようになります。

// load_time_ordered.stan
data {
  int<lower=0> N;
  int<lower=0> cv[N];
  int<lower=0> session[N];
}
parameters {
  ordered[N] cvr_rev;
}
transformed parameters {
  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- inv_logit(cvr_rev[N - i + 1]);
  }
}
model {
  cv ~ binomial(session, cvr);
}
library(rstan)
model <- stan_model("load_time_orderd.stan")
standata <- list(N = length(cv), cv = cv, session = session)
result <- sampling(model, data = standata, chains = 3)

非常にきれいにパラメータを推定することができました。

注意点

まず、parameters ブロックにおいて、ordered 型のパラメータを作成しています。

  ordered[N] cvr_rev;

この際、ordered 型は、「小さい順」という制約であり、求めたいパラメータは大きい順なので、これを逆順にする必要があります。
この逆順にするという操作を行っているのが transformed parameters ブロックです。
ここで、新しいパラメータ cvr に cvr_rev を逆からひとつひとつ代入する、つまり、

  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- cvr_rev[N - i + 1];
  }

のようにすれば良いように思いますが、それではうまくいきません。
ordered 型には lower や upper のような最小値・最大値の制約がつけられないため、−∞〜+∞までの値を取ります。
しかし、cvr には 0〜1 という制約が付いているで、cvr_rev を適当にサンプリングすると、ほとんどがはみ出してしまいます。
このような場合に便利なのが inv_logit() です。
これは、logit 関数の逆関数(logistic 関数!)で、−∞〜+∞ の値を 0〜1 にマッピングしてくれます。
したがって、

  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- inv_logit(cvr_rev[N - i + 1]);
  }

と書けば、cvr_rev がどのような値を取ろうが、cvr は 0〜1 の値となり、cvr の制約条件に違反しなくなります。

まとめ

Stan でパラメータに大小関係の制約を付ける方法を紹介しました。
このように、パラメータ間の制約を柔軟に加えることができるのも、ベイズ統計モデリング(Stan)の魅力かなと思います。
Enjoy!