はじめに

こんにちは、事業会社で働いているデータサイエンティストです。

最近の業務では、経営の意思決定における資源配分の最適化という、古典的なミクロ経済学の課題に取り組むことが中心になっています。

前回の記事では、Stanで推定したベイズMMMモデルを利用して、広告の最適予算配分を推定する方法を説明しました：

今回の記事では、ミクロ経済学由来のコブ＝ダグラス型（Cobb-Douglas）売上関数に、さらにガウス過程というベイズ機械学習の要素を入れて、長期間の広告予算を最適配分する方法を説明します。

シミュレーションシナリオ

本記事では、売上が以下のように生成されると仮定します：

$$
売上 = OOH広告^{\beta_{OOH}} \cdot TVCM^{\beta_{TVCM}} \cdot トレンド要素 \cdot 周期要素 \cdot \epsilon
$$

ここで、$\epsilon$ は独立かつ同一に分布する対数正規分布に従うノイズ項です。

広告出稿量に関する情報は機密性が高いため、本物のデータをそのまま用いて記事を書くことはできません。しかし、ガウス過程を必要とするほどのトレンド性・周期性が含まれないデータでは、現実味に欠けるシナリオとなってしまいます。

そこで、Googleトレンドからあるキーワードの検索ボリューム（以下 search）を取得し、それを用いて以下のように合成データを作成しました：

$$
売上 = OOH広告^{0.3} \cdot TVCM^{0.7} \cdot search \cdot \epsilon
$$

この設定のもとで、searchの推移をガウス過程でトレンドと周期的な変動に要因分解します。

このガウス過程を用いたモデルが$\beta$やsearchのトレンドと周期性を適切に推定できるか、また最適化を通じてどのような結果が得られるのかを検証してみましょう！

データ準備

本記事で使用するデータは、以下のようにして生成しました。

まず、Googleトレンドから取得した検索ボリュームのデータを、そのまま保存したものが timeseries.csv です。どのキーワードを使用したかについては、申し訳ありませんが非公開とさせていただきます。ご了承ください。

set.seed(12345)
sales_data <- readr::read_csv("timeseries.csv", skip = 2) |>
  `colnames<-`(c("date", "latent")) |>
  dplyr::mutate(
    # latentが記事内のsearchにあたる。潜在変数だからlatentにした
    latent = exp((latent-mean(latent))/sd(latent)),
    ooh = abs(rnorm(dplyr::n())),
    tvcm = abs(rnorm(dplyr::n())),
    e = exp(rnorm(dplyr::n())),
    sales = (ooh ^ 0.3) * (tvcm ^ 0.7) * latent * e
  )

モデル推定

ここでは、下記のStan言語のコードで、モデルを推定します；

cobb_douglas_with_gp.stan

data {
  int time_type;
  array[time_type] real time_seq;
  
  int N;
  
  vector[N] a;
  vector[N] b;
  array[N] real y;
}
parameters {
  real<lower=0> rho_quad;
  real<lower=0> alpha_quad;
  real<lower=0> sigma_quad;
  vector[time_type] eta_quad;
  
  real<lower=0> rho_cycle;
  real<lower=0> alpha_cycle;
  real<lower=0> sigma_cycle;
  vector[time_type] eta_cycle;
  
  real intercept;
  simplex[2] beta;
  real<lower=0> sigma;
}
transformed parameters {
  vector[time_type] f_trend;
  vector[time_type] f_cycle;
  {
    matrix[time_type, time_type] L_K;
    matrix[time_type, time_type] K = gp_exp_quad_cov(time_seq, alpha_quad, rho_quad);

    for (n in 1:time_type) {
      K[n, n] = K[n, n] + sigma_quad;
    }

    L_K = cholesky_decompose(K);
    f_trend = L_K * eta_quad - mean(L_K * eta_quad);
  }
  {
    matrix[time_type, time_type] L_K;
    matrix[time_type, time_type] K = gp_periodic_cov(time_seq, alpha_cycle, rho_cycle, 52.0/time_type);

    for (n in 1:time_type) {
      K[n, n] = K[n, n] + sigma_cycle;
    }

    L_K = cholesky_decompose(K);
    f_cycle = L_K * eta_cycle - mean(L_K * eta_cycle);
  }
}
model {
  rho_quad ~ inv_gamma(5, 5);
  alpha_quad ~ normal(0, 1);
  sigma_quad ~ normal(0, 1);
  eta_quad ~ normal(0, 1);
  
  rho_cycle ~ inv_gamma(5, 5);
  alpha_cycle ~ normal(0, 1);
  sigma_cycle ~ normal(0, 1);
  eta_cycle ~ normal(0, 1);
  
  intercept ~ normal(0, 10);
  beta ~ dirichlet(rep_vector(0.5, 2));
  sigma ~ inv_gamma(0.1, 0.1);
  log(y) ~ normal(intercept + f_trend[1:N] + f_cycle[1:N] + beta[1] * log(a) + beta[2] * log(b), sigma);
}
generated quantities {
  vector[time_type] f_total;
  f_total = f_trend + f_cycle;
}

実際にデータを渡して推定しましょう！ここでは、最後の6期分のデータを検証データとして残し、この後の広告予算配分の対象期間とします：

m_init <- cmdstanr::cmdstan_model("cobb_douglas_with_gp.stan")

m_estimate <- m_init$variational(
  data = list(
    time_type = nrow(sales_data),
    time_seq = (1:nrow(sales_data))/nrow(sales_data),
    
    N = nrow(sales_data) - 6,
    a = sales_data$ooh[1:(nrow(sales_data) - 6)],
    b = sales_data$tvcm[1:(nrow(sales_data) - 6)],
    y = sales_data$sales[1:(nrow(sales_data) - 6)]
  )
)


m_summary <- m_estimate$summary()

まず、$\beta$から確認しましょう：

> m_summary |> dplyr::filter(stringr::str_detect(variable, "^beta\\["))
# A tibble: 2 × 7
  variable  mean median     sd    mad    q5   q95
  <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl>
1 beta[1]  0.270  0.265 0.0484 0.0487 0.198 0.352
2 beta[2]  0.730  0.735 0.0484 0.0487 0.648 0.802

0.3と0.7に近い値になっています。

次に、トレンドと周期性のガウス過程がそれぞれどんな形状の関数を学習したのかを確認しましょう：

g_trend <- m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^f_trend\\[")) |>
  dplyr::bind_cols(
    date = sales_data$date
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::labs(
    title = "trend"
  )

g_cycle <- m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^f_cycle\\[")) |>
  dplyr::bind_cols(
    date = sales_data$date
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::labs(
    title = "cycle"
  )

g_total <- m_summary |>
  dplyr::filter(stringr::str_detect(variable, "^f_total\\[")) |>
  dplyr::bind_cols(
    date = sales_data$date
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::labs(
    title = "trend + cycle"
  )

ここで、patchworkパッケージで少し見やすく可視化します：

library(patchwork)
(g_trend + g_cycle)/g_total

トレンド全体としては、明確な下降傾向が見られます。また、周期的なパターンとしては、年末年始にかけて検索ボリュームが落ち込み、3月から4月にかけて年間のピークを迎えるという特徴が確認できます。キーワードを非公開にしているため、このような分析に再現性がない点は承知していますが、それでもドメイン知識的には非常に納得感のある結果となっています。

広告予算の最適化

こちらのStan言語のコードで最適化します：

cobb_douglas_with_gp_predict.stan

data {
  int time_type;
  array[time_type] real time_seq;
  
  int new_N;
  array[new_N] int new_time;
  array[new_N] real new_a;
  array[new_N] real new_b;
}
parameters {
  real<lower=0> rho_quad;
  real<lower=0> alpha_quad;
  real<lower=0> sigma_quad;
  vector[time_type] eta_quad;
  
  real<lower=0> rho_cycle;
  real<lower=0> alpha_cycle;
  real<lower=0> sigma_cycle;
  vector[time_type] eta_cycle;
  
  real intercept;
  simplex[2] beta;
  real<lower=0> sigma;
}
transformed parameters {
  vector[time_type] f_trend;
  vector[time_type] f_cycle;
  {
    matrix[time_type, time_type] L_K;
    matrix[time_type, time_type] K = gp_exp_quad_cov(time_seq, alpha_quad, rho_quad);

    for (n in 1:time_type) {
      K[n, n] = K[n, n] + sigma_quad;
    }

    L_K = cholesky_decompose(K);
    f_trend = L_K * eta_quad - mean(L_K * eta_quad);
  }
  {
    matrix[time_type, time_type] L_K;
    matrix[time_type, time_type] K = gp_periodic_cov(time_seq, alpha_cycle, rho_cycle, 52.0/time_type);

    for (n in 1:time_type) {
      K[n, n] = K[n, n] + sigma_cycle;
    }

    L_K = cholesky_decompose(K);
    f_cycle = L_K * eta_cycle - mean(L_K * eta_cycle);
  }
}
generated quantities {
  vector[new_N] predict;
  for (i in 1:new_N){
    predict[i] = exp(normal_rng(intercept + f_trend[new_time[i]] + f_cycle[new_time[i]] + beta[1] * log(new_a[i]) + beta[2] * log(new_b[i]), sigma));
  }
}

最適化のために、上記のStan言語のコードをコンパイルして、必要な関数群を作ります：

gq_init <- cmdstanr::cmdstan_model("cobb_douglas_with_gp_predict.stan")

predict_cobb_douglas_gp <- function(S){
  gq_init$generate_quantities(
    fitted_params = m_estimate,
    data = list(
      time_type = nrow(sales_data),
      time_seq = (1:nrow(sales_data))/nrow(sales_data),
      
      new_N = 6,
      new_time = (nrow(sales_data) - 6 + 1):nrow(sales_data),
      # S[1:6]にOOHの出稿量が入る
      new_a = S[1:6],
      # S[7:12]にTVCMの出稿量が入る
      new_b = S[7:12]
    )
  )$
    summary() |>
    dplyr::pull(mean) |> 
    sum()
}

compute_cost <- function(S){
  OOH広告の単価が30円（仮の単位）、TVCMは10円
  30 * sum(S[1:6]) + 15 * sum(S[7:12])
}

さて、ここでは、6期間（6週間）分の広告予算が合計で100円（仮の単位）と仮定し、推定された売上予測モデルをもとに、各期間におけるOOHおよびTVCMの出稿量を最適化してみましょう。

result <- Rsolnp::solnp(
  # 初期値は適当に入れました、、、
  pars = rep(10, 12),
  # 最小化を想定しているため、マイナスをつけましょう
  fun = \(S) - predict_cobb_douglas_gp(S),
  eqfun = \(S) compute_cost(S),
  eqB = 100,
  # 出稿量がマイナスになってはいけない
  ineqfun = \(S){S},
  ineqLB = rep(0, 12),
  ineqUB = rep(Inf, 12),
  control = list(
    outer.iter = 5,
    inner.iter = 5
  )
)

この計算は本当にかなり時間がかかります。

最後に、この6期内の最適な広告出稿量を確認しましょう：

tibble::tibble(
  date = sales_data$date[(nrow(sales_data) - 6 + 1):nrow(sales_data)],
  ooh = result$pars[1:6],
  tvcm = result$pars[7:12]
) |>
  tidyr::pivot_longer(!date, names_to = "media", values_to = "amount") |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = amount, color = media))

結局、$\beta$の大きいTVCMの出稿量が全体的に多くなるという結果になっています。これは必ずしも不合理な結果ではありませんが、出稿量の推移を見ると少々ギザギザしており、最適化がうまく機能していない可能性があります。

この問題に対処する一つの方法として、Rsolnp::solnpのcontrol引数に渡すlist内で最大計算回数を増やすことが考えられます。しかし、本質的な解決を目指すのであれば、Rの他パッケージやPythonなど他言語の最適化アルゴリズムやソルバーを活用する方が適切かもしれません。

また、ガウス過程によって推定されたトレンドと周期性を確認すると、予算配分の対象となる最後の6期はほぼ平坦であり、「3ヶ月後にトレンドと周期性のピークが重なる。このタイミングでリソースを集中投下しよう！」といった戦略的な配分がサジェストされるかどうかを検証することはできませんでした。

結論

いかがでしたか？

本記事では、コブ＝ダグラス型の売上関数にガウス過程を組み合わせたベイズモデルを提案し、その有用性を検証しました。

推定結果および最適化アルゴリズムの挙動からも明らかなように、ガウス過程などの機械学習的構造と、コブ＝ダグラス型関数の構造を組み合わせた本モデルは、構造的な解釈可能性と柔軟性を両立する点で非常に有望です。加えて、これは厳密な数学的裏付けがあるわけではなくあくまで推測に過ぎませんが、コブ＝ダグラス型の構造を採用したことにより、完全にブラックボックス化されたモデル（決定木系の手法や深層学習など）に比べて、最適化が比較的容易になっている可能性があります。

また、このモデルを応用することで、将来の広告予算を配分する際に、市況の変動（ガウス過程が学習したトレンドと周期性）を考慮することができるため、将来の広告効果を間接的に予測していると捉えることができます。将来の広告効果は、経営会議でも定番の話題の一つです。データサイエンスの力で経営陣の意思決定を支援し、より質の高い広告戦略形成に貢献できれば、データサイエンスチームは経営から一目置かれる存在となるはずです。

私自身も、まさにそれを目指して、頼もしいメンバーたちとともに日々working hardしています。

最後に、私たちと一緒にデータサイエンスの力を使って社会を改善していきたい方はぜひこちらのページをご確認ください：