混合負の二項分布による広告フリークエンシー分布の推定

こんにちは。TVer の広告事業領域でデータサイエンティストをしている川井です。普段は TVer 広告の配信システムの開発や、広告効果分析、データ基盤構築などを担当しています。

ストリーミング広告の配信において、「誰に何回広告が届いたか」を把握することは非常に重要です。同じユーザーに繰り返し広告が届くよりも、まだ届いていないユーザーに届けたい——そうしたニーズに応えるには、広告接触回数の分布構造を理解する必要があります。

以前の記事「負の二項分布でストリーミング広告のリーチを予測してみた」では、リーチカーブに直接フィットさせるアプローチを紹介しました。今回は別のアプローチとして、フリークエンシー分布（広告接触回数の分布）そのものを混合モデルで分解する手法を紹介します。単一の負の二項分布では「ユーザー全体の異質性」を 1 つのパラメータに集約しますが、混合モデルを使うと「ライト層が○%、ヘビー層が○%」といったユーザー構造を可視化できます。フリークエンシー分布の「なぜこの形になるのか」を解釈しやすくなる点が利点です。

フリークエンシー分布とは

「各ユーザーが何回広告に接触したか」という分布をフリークエンシー (frequency; FQ) 分布といいます。典型的なフリークエンシー分布には以下の特徴があります。

が非常に多い
- 母集団の大半は広告に接触しない
急激に減衰
- 広告に $1$ 回だけ接触したユーザー（ $FQ=1$ ）が最も多く、 $FQ=2$ 以降のユーザー数は急激に減衰する
右に裾が長い
- ごく少数のユーザーが非常に多くの回数接触する場合がある

この分布を統計モデルで表現できれば、「あと○○万回配信できたら何人リーチが増える」といった予測が可能になります。

シミュレーションデータの生成

フリークエンシー分布を推定するため、シミュレーションデータで検証してみましょう。以下のような 3 タイプのユーザーが混在する母集団（100万人）を想定します（ユーザータイプ・構成比・パラメータの数値は全て仮定です）。

フリークエンシー分布は過分散（分散が平均より大きい）を示すことが多く、ポアソン分布では表現できません。負の二項分布はこの過分散を扱えるため、広告接触回数のモデリングに適しています。負の二項分布のパラメータには $\mu$ と $\theta$ があり、 $\mu$ は平均、 $\theta$ は過分散を制御するパラメータです。 $\theta$ が小さいほど分散が大きくなり、 $\theta \longrightarrow \infty$ でポアソン分布に収束します。

ユーザータイプ	構成比	$\mu$	$\theta$	特徴
ライト	60%	0.1	0.5	広告にほとんど接触しない
レギュラー	30%	1.0	0.5	平均 1 回程度接触する
ヘビー	10%	3.0	0.5	頻繁に広告に接触する

負の二項分布の分散は $\mu + \frac{\mu^{2}}{\theta}$ で計算されるため、 $\theta$ が小さいほど分散が大きくなります。今回は簡単のため、全てのユーザータイプで $\theta$ は一律 $0.5$ としました。

なお、私は最近 R にハマっているので、以降は全て R を用いて実験をしています。

Code

set.seed(42)

n_population <- 1000000

# 真のパラメータ
true_params <- tribble(
  ~component, ~prior,  ~mu, ~theta,
  "ライト",     0.60,  0.1,    0.5,
  "レギュラー", 0.30,  1.0,    0.5,
  "ヘビー",     0.10,  3.0,    0.5,
)

sim_data <- true_params |>
  slice_sample(n = n_population, weight_by = prior, replace = TRUE) |>
  mutate(
    user_id = row_number(),
    imp = rnbinom(n(), mu = mu, size = theta)
  )

成分ごとの統計量を確認します。avg_fq は理論上の $\mu$ に、variance は $\mu + \frac{\mu^{2}}{\theta}$ に近い値となるはずです。

Code

sim_data |>
  summarise(
    each_uu = n(),
    each_imp = sum(imp),
    reach_uu = sum(imp > 0),
    avg_fq = each_imp / each_uu,
    variance = var(imp),
    landing_fq = each_imp / reach_uu,
    .by = component
  ) |>
  arrange(each_imp) |>
  gt::gt() |>
  gt::fmt_integer(columns = 2:4) |>
  gt::fmt_number(columns = 5:7, decimals = 2)

component	each_uu	each_imp	reach_uu	avg_fq	variance	landing_fq
ライト	600,076	59,908	52,227	0.10	0.12	1.15
レギュラー	300,033	300,699	127,029	1.00	3.01	2.37
ヘビー	99,891	300,734	62,062	3.01	21.35	4.85

可視化してみましょう。接触回数が $0$ 回のユーザーがほとんどなため、可視化からは除外し、縦軸を対数表記にしています。

Code

sim_data |>
  filter(imp > 0) |>
  mutate(component = factor(component, levels = c("ライト", "レギュラー", "ヘビー"))) |>
  bind_rows(
    sim_data |> filter(imp > 0) |> mutate(component = "全体")
  ) |>
  mutate(component = factor(component, levels = c("全体", "ライト", "レギュラー", "ヘビー"))) |>
  ggplot(aes(x = imp, fill = component)) +
  geom_histogram(binwidth = 1, alpha = 0.5) +
  scale_x_continuous(breaks = scales::breaks_width(10)) +
  scale_y_log10(labels = scales::comma) +
  facet_wrap(~ component, ncol = 1) +
  labs(x = "接触回数", y = "人数")

全体としては右に裾の長いフリークエンシー分布になっています。ライト層は数回の広告接触に留まるのに対し、ヘビー層には広告にたくさん接触しているユーザーが存在することがわかります。なんとなくそれっぽいデータが生成されましたね。

全体の統計量についても確認します。今回はシミュレーションデータなので上記のように内訳がわかりますが、実際の広告配信後に得られるデータは以下のようなものになります。

Code

sim_data |>
  summarise(
    total_uu = n(),
    total_imp = sum(imp),
    reach_uu = sum(imp > 0),
    avg_fq = total_imp / total_uu,
    variance = var(imp),
    landing_fq = total_imp / reach_uu,
    max_fq = max(imp)
  ) |>
  gt::gt() |>
  gt::fmt_integer(columns = 1:3) |>
  gt::fmt_number(columns = 4:6, decimals = 2)

total_uu	total_imp	reach_uu	avg_fq	variance	landing_fq	max_fq
1,000,000	661,341	241,318	0.66	3.88	2.74	79

つまり、母集団サイズ、総インプレッション、リーチ人数、そしてフリークエンシー分布（接触回数ごとの人数）です。この情報だけから、背後にあるユーザー構造（ライト・レギュラー・ヘビー）とその構成比を推定できるでしょうか。

単一の負の二項分布によるフィッティング

まず、単一の負の二項分布でフリークエンシー分布をモデリングしてみます。

ここでは glmmTMB パッケージを使用します。 glmmTMB は一般化線形混合モデル（GLMM）をフィッティングするためのパッケージで、負の二項分布を含む様々な分布族をサポートしています。今回は混合効果を使わないシンプルなケースですが、weights 引数で集計データを扱える点が便利です。weights = cnt とすることで、「接触回数 imp のユーザーが cnt 人いる」という集計データをそのまま渡せます。実務では、ベイズモデリングを使うと推定の不確実性を信用区間で表現できて有用なことが多いですが、今回は真のパラメータがわかっているシミュレーションなので、計算の速い頻度論的アプローチで進めます。

Code

freq_dist <- sim_data |> count(imp, name = "cnt")

fit_single <- glmmTMB::glmmTMB(
  data = freq_dist,
  formula = imp ~ 1,
  family = glmmTMB::nbinom2,
  weights = cnt
)

summary(fit_single)

 Family: nbinom2  ( log )
Formula:          imp ~ 1
Data: freq_dist
Weights: cnt

      AIC       BIC    logLik -2*log(L)  df.resid
1941961.4 1941965.5 -970978.7 1941957.4        57


Dispersion parameter for nbinom2 family (): 0.174

Conditional model:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.413487   0.002695  -153.4   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

モデルは正常に収束しました。推定されたパラメータを確認します。

Code

single_params <- get_parameters(fit_single) |>
  mutate(value = exp(Estimate)) |>
  select(Component, value) |>
  pivot_wider(names_from = Component, values_from = value) |>
  rename(mu = conditional, theta = dispersion)

single_params |>
  gt::gt() |>
  gt::fmt_number(decimals = 2)

mu	theta
0.66	0.17

$\mu$ の推定値
- 母集団平均は約 $0.66$ なので、正しく捉えられています。
$\theta$ の推定値
- 真の $\theta$ は成分ごとに $0.5, 0.5, 0.5$ ですが、推定値は $0.17$ でした。
- 単一分布で混合分布の過分散を表現しようとした結果、 $\theta$ が不自然に小さくなってしまいました。

このパラメータを使って、推定されたフリークエンシー分布と実際の分布を比較してみましょう。

Code

freq_dist |>
  mutate(
    total_cnt = sum(cnt),
    expected = total_cnt * dnbinom(imp, mu = single_params$mu, size = single_params$theta),
    residual = cnt - expected
  ) |>
  filter(imp <= 20) |>
  ggplot(aes(x = imp, y = residual)) +
  geom_col(fill = "steelblue") +
  geom_hline(yintercept = 0, color = "coral") +
  scale_x_continuous(breaks = scales::breaks_width(1)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "接触回数", y = "残差（実測 - 推定）")

接触回数が 1 回の部分で大きな残差が生じており、単一の負の二項分布ではフリークエンシー分布の構造を十分に捉えられていないことがわかります。

混合負の二項分布によるフィッティング

単一の負の二項分布では残差が大きく、フリークエンシー分布の構造を十分に捉えられていませんでした。ここでは、複数の負の二項分布を混合したモデルを試してみます。混合負の二項分布のフィッティングには、flexmixパッケージを使用します。flexmix は有限混合モデルをフィッティングするためのパッケージで、countreg パッケージの FLXMRnegbin と組み合わせることで混合負の二項分布を推定できます。

成分数を 2〜4 で試し、BIC で比較します。BIC（ベイズ情報量規準）はモデルの当てはまりと複雑さのバランスを評価する指標で、値が小さいほど良いモデルとされます。

stepFlexmix の主な引数は以下の通りです。

k: 試す成分数。ベクトルで複数指定可能。
model: 各成分の分布。FLXMRnegbin() で負の二項分布を指定。
weights: 集計データの頻度ウェイト
nrep: 初期値を変えて繰り返す回数。EMアルゴリズムは初期値依存なので、複数回試行して最良の結果を採用
control: minprior は成分の最小構成比。これを下回る成分は除外される

Code

tidy_flexmix <- function(fit) {
   flexmix::parameters(fit) |>
    t() |>
    as_tibble() |>
    rename(log_mu = `coef.(Intercept)`) |>
    mutate(
      component = row_number(),
      prior = flexmix::prior(fit),
      mu = exp(log_mu),
      .before = 1
    )
}

set.seed(42)

fit_mix <- flexmix::stepFlexmix(
  data = freq_dist,
  formula = imp ~ 1,
  k = 2:4,
  model = countreg::FLXMRnegbin(),
  weights = ~ cnt,
  control = list(iter.max = 500, minprior = 0.01),
  nrep = 20,
  verbose = FALSE
)

tibble(
  k = 2:4,
  BIC = sapply(fit_mix@models, BIC)
) |>
  gt::gt() |>
  gt::fmt_integer()

k	BIC
2	1,936,738
3	1,936,448
4	1,936,502

$k=3$ で BIC が最小となっています。 $k=4$ では改善せずむしろ悪化しており、3 成分が最適であることを示しています。今回のシミュレーションでは真の成分数が 3 だったので、正しく特定できています。 BIC が最小だった、 $k=3$ のモデルのパラメータを確認します。

Code

fit_best <- flexmix::getModel(fit_mix, which = "BIC")
mix_params <- tidy_flexmix(fit_best)
mix_params |>
  gt::gt() |>
  gt::fmt_number(columns = 2:5, decimals = 2)

component	prior	mu	log_mu	theta
1	0.12	2.67	0.98	0.46
2	0.20	1.22	0.20	0.66
3	0.67	0.12	−2.12	0.40

真のパラメータと比較したのが以下です。

真の成分	構成比	$\mu$	$\theta$	→	推定	構成比	$\mu$	$\theta$
ライト	0.60	0.1	0.50	→	Comp3	0.67	0.12	0.40
レギュラー	0.30	1.0	0.50	→	Comp2	0.20	1.22	0.66
ヘビー	0.10	3.0	0.50	→	Comp1	0.12	2.67	0.46

完全には一致しませんでしたが、構造の傾向は捉えられています。 $\mu$ の大小関係（ライト < レギュラー < ヘビー）は正しく識別されており、構成比もライト層が多数を占めるという特徴を反映しています。ただし、レギュラー層の構成比が過小推定され、その分ライト層に吸収されています。中間的な成分は境界が曖昧になりやすく、隣接成分と混同されやすい傾向があります。

混合負の二項分布の残差確認

推定されたパラメータを使って、フリークエンシー分布の残差を確認します。

Code

freq_dist |>
  crossing(mix_params) |>
  mutate(prob = prior * dnbinom(imp, mu = mu, size = theta)) |>
  summarise(cnt = max(cnt), prob = sum(prob), .by = imp) |>
  mutate(
    total_cnt = sum(cnt),
    expected = total_cnt * prob,
    residual = cnt - expected
  ) |>
  filter(imp <= 20) |>
  ggplot(aes(x = imp, y = residual)) +
  geom_col(fill = "steelblue") +
  geom_hline(yintercept = 0, color = "coral") +
  scale_x_continuous(breaks = scales::breaks_width(1)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "接触回数", y = "残差（実測 - 推定）")

単一の負の二項分布と比較して、残差が大幅に改善されています。 $FQ=1$ 付近の大きなズレが解消され、全体的に残差が小さくなっています。 $FQ=3, FQ=4$ で若干の過大推定がありますが、母集団 100 万人に対して誤差 200 ~300 人程度なので実用上は問題ないレベルです。

リーチ予測への応用

推定したパラメータを使って、任意のインプレッション数に対するリーチを予測できます。混合負の二項分布では、リーチ率（1人以上に接触する確率）は以下のように計算されます。

$\displaystyle \text{リーチ率} = 1 - \sum_{i} \pi_i \cdot P(X=0 | \mu_i, \theta_i)$

ここで $\pi_i$ は各成分の構成比、 $P(X=0)$ は負の二項分布で接触回数が $0$ となる確率です。

インプレッション数を変化させたときのリーチ予測を行うには $\mu$ をスケーリングします。配信量が 2 倍になれば、各ユーザーへの接触機会も2倍になると仮定し、 $\mu$ を比例させてリーチカーブを描いてみました。

Code

calc_reach <- function(target_imp, params_df, total_uu, actual_imp) {
  scale <- target_imp / actual_imp

  params_df |>
    mutate(
      scaled_mu = mu * scale,
      p_zero = prior * dnbinom(0, mu = scaled_mu, size = theta)
    ) |>
    summarise(p_zero = sum(p_zero)) |>
    pull(p_zero) |>
    {\(p) total_uu * (1 - p)}()
}

freq_dist |>
  summarise(
    actual_imp = sum(imp * cnt),
    actual_reach = sum(cnt[imp > 0]),
    total_uu = sum(cnt)
  ) |>
  crossing(target_imp = seq(10000, 1000000, by = 10000)) |>
  mutate(
    predicted_reach = pmap_dbl(
      list(target_imp, total_uu, actual_imp),
      \(target_imp, total_uu, actual_imp) calc_reach(target_imp, mix_params, total_uu, actual_imp)
    ),
    reach_rate = predicted_reach / total_uu
  ) |>
  ggplot(aes(x = target_imp, y = predicted_reach)) +
  geom_line(color = "steelblue") +
  geom_point(aes(x = actual_imp, y = actual_reach), color = "coral", size = 3) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "インプレッション数", y = "リーチ人数")

点は実測値（約66万 imp → 約24万リーチ）で、予測カーブが実測値を通過しており、モデルの妥当性が確認できます。また、カーブの形状からリーチの逓減効果も見て取れます。インプレッション数を増やしても、リーチの伸びは徐々に緩やかになります。これは一部のヘビー層が繰り返しインプレッションを吸収するためで、混合モデルがこの構造を捉えていることを示しています。

このように、パラメータさえ推定できれば「○○万 imp でリーチ何人？」という問いに即答できるようになります。

まとめ

今回は、フリークエンシー分布を混合負の二項分布でモデリングする手法を紹介しました。シミュレーションを通じて、以下のことが確認できました。

単一の負の二項分布では、異なるユーザー層が混在するフリークエンシー分布を十分に表現できない
混合負の二項分布を用いることで、ライト・レギュラー・ヘビーといったユーザー構造を識別できる
推定したパラメータから、任意のインプレッション数に対するリーチ予測が可能になる

この手法を実データに応用することで、リーチ予測やユーザー構造の定量的な比較が可能になるほか、ターゲティング戦略やフリークエンシーキャップ設定の検討に活用できそうです。

なお、今回はシミュレーションデータを使用しましたが、実際のキャンペーンデータに適用する際には、配信期間やターゲティング条件によるパラメータの変動など、より突っ込んだ検討が必要になります。

We’re Hiring!

TVer では、データサイエンスの力で広告事業の成長を支えるメンバーを募集しています。今回紹介したようなリーチ予測モデルの構築・改善や、広告配信の最適化など、チャレンジングな課題に取り組んでいます。また、広告データ基盤も現在絶賛整備中のため、データエンジニアの募集も行っています。

ご興味のある方は、ぜひカジュアル面談からでも気軽にお話しましょう！