MMMモデルを推定した。で、結局どうやって予算最適化すればいいの？

# はじめに

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

https://dippeople.dip-net.jp/19888/

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

私自身もそうですが、ミクロ経済学のバックグラウンドを持つ人は、往々にして生産・効用最大化の問題を紙とペンで解くことには慣れています。しかし、実務で扱うような複雑な生産関数や効用関数を前にすると、紙とペンだけで最適な資源配分を導出するのはほぼ不可能です。

それにもかかわらず、残念ながらミクロ経済学の教育では、コンピュータを用いた最適化手法を体系的に学ぶ機会はほとんどありません。

本記事では、そうした実務と理論研究のギャップを感じている方に向けて、R言語とStanを用いてベイズ統計モデルを推定し、その結果を活用して資源配分の最適化を行う方法を紹介します。

ミクロ経済学とMMMの関連性について、こちらの記事をご参照ください：

https://qiita.com/Gotoubun_taiwan/items/b92ea7bed5f1a1cbf66b

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

ここでは、会社の売上が次のような構造で生成されていると仮定します：

$$
売上 = OOH広告^{\beta_{OOH}} \times TVCM^{\beta_{TVCM}} \times \epsilon
$$

ここで、$\epsilon$ は対数正規分布に従うノイズ項です。このような形式の関数は、コブ＝ダグラス型生産関数（Cobb-Douglas production function）として知られています。

この記事では、経済学以外のバックグラウンドを持つ読者の方にも理解していただけるように、コブ＝ダグラス型生産関数とは何か、そしてそこから**予算制約の下での最適な広告出稿配分**をどのように導くかを、丁寧に説明していきます。

経済学に馴染みのある方にとっては、ちょうど良い復習の機会として読んでいただければと思います。

なお、広告にかけられる予算の制約は、次のように表現できます：


$$
予算制約 >= OOH広告単価 \times OOH広告 + TVCM単価 \times TVCM
$$

なお、ここでは制約条件を「$≥$」と表記していますが、通常は予算をすべて使い切ることが多いため、「$=$」と置いても問題ありません。

## 最適化問題を数式で解く

さて、マーケティング部門の責任者が直面する制約条件付きの最大化問題は、次のように表現できます：

$$
\max \quad OOH広告^{\beta_{OOH}} \times TVCM^{\beta_{TVCM}}
$$

$$
\text{subject to} \quad 予算制約 \geq OOH広告単価 \times OOH広告 + TVCM単価 \times TVCM
$$

ここで、$\epsilon$ はノイズ項であり、マーケティング部門の責任者はその値を事前に知ることができないため、この最適化問題からは除外しています。

このような制約条件付き最大化問題は、ラグランジュの未定乗数法（Lagrange multiplier method）を用いることで、通常の最大化問題へと変換できます。

ラグランジュ関数は以下のように定義されます：

$$
L = OOH広告^{\beta_{OOH}} \times TVCM^{\beta_{TVCM}} + 
$$


$$
\lambda \left(予算制約 - OOH広告単価 \times OOH広告 - TVCM単価 \times TVCM \right)
$$

それでは、まずこのラグランジュ関数に対して一階条件を求めてみましょう。

OOH広告についての偏微分は次の通りです：

$$
\frac{\partial L}{\partial OOH広告} = \beta_{OOH} \cdot OOH広告^{\beta_{OOH} - 1} \cdot TVCM^{\beta_{TVCM}} - \lambda \cdot OOH広告単価 = 0
$$

整理すると：

$$
\lambda = \frac{\beta_{OOH} \cdot OOH広告^{\beta_{OOH} - 1} \cdot TVCM^{\beta_{TVCM}}}{OOH広告単価}
$$

同様に、TVCMについては：

$$
\frac{\partial L}{\partial TVCM} = \beta_{TVCM} \cdot OOH広告^{\beta_{OOH}} \cdot TVCM^{\beta_{TVCM} - 1} - \lambda \cdot TVCM単価 = 0
$$

整理すると：

$$
\lambda = \frac{\beta_{TVCM} \cdot OOH広告^{\beta_{OOH}} \cdot TVCM^{\beta_{TVCM} - 1}}{TVCM単価}
$$

最後に、$\lambda$について微分すると：

$$
\frac{\partial L}{\partial \lambda} = 予算制約 - OOH広告単価 \times OOH広告 - TVCM単価 \times TVCM = 0
$$

ここからわかるように：

$$
\frac{\beta_{OOH} \cdot OOH広告^{\beta_{OOH} - 1} \cdot TVCM^{\beta_{TVCM}}}{OOH広告単価} = \frac{\beta_{TVCM} \cdot OOH広告^{\beta_{OOH}} \cdot TVCM^{\beta_{TVCM} - 1}}{TVCM単価}
$$

$$
\frac{TVCM単価}{OOH広告単価} = \frac{\beta_{TVCM} \cdot OOH広告^{\beta_{OOH}} \cdot TVCM^{\beta_{TVCM} - 1}}{\beta_{OOH} \cdot OOH広告^{\beta_{OOH} - 1} \cdot TVCM^{\beta_{TVCM}}}
$$

$$
\frac{TVCM単価}{OOH広告単価} = \frac{\beta_{TVCM} \cdot OOH広告}{\beta_{OOH} \cdot TVCM}
$$

この式を整理すると、最適なOOH広告とTVCMの投入量の比率がわかります：

$$
\frac{TVCM}{OOH広告} = \frac{\beta_{TVCM} \cdot OOH広告単価}{\beta_{OOH} \cdot TVCM単価}
$$

すなわち、最適な広告配分は各広告の価格と限界貢献度（限界生産力、$\beta$）のバランスで決まることが分かります。

さらに整理すると、支出金額（単価×利用量）の比に変換できます：

$$
\frac{TVCM \cdot TVCM単価}{OOH広告 \cdot OOH広告単価} = \frac{\beta_{TVCM}}{\beta_{OOH}}
$$

つまり、支出ベースで表すと：

$$
\frac{TVCM出稿金額}{OOH広告出稿金額} = \frac{\beta_{TVCM}}{\beta_{OOH}}
$$

要するに、TVCMの売上貢献（$\beta_{TVCM}$）が大きければ、最適なTVCM広告への支出金額も相対的に大きくなるということです。これが経済学における**限界効率に基づく資源配分**の考え方です。

この関係式と予算制約式：

$$
OOH広告単価 \cdot OOH広告 + TVCM単価 \cdot TVCM = 予算制約
$$

を連立することで、需要関数（最適なOOH広告とTVCMの水準）を明示的に導出できます。

### 需要関数

まず、上記の最適投入比：

$$
TVCM = \frac{\beta_{TVCM} \cdot OOH広告単価}{\beta_{OOH} \cdot TVCM単価} \cdot OOH広告
$$

を予算制約に代入します：

$$
OOH広告単価 \cdot OOH広告 + TVCM単価 \cdot \left( \frac{\beta_{TVCM} \cdot OOH広告単価}{\beta_{OOH} \cdot TVCM単価} \cdot OOH広告 \right) = 予算制約
$$

$$
OOH広告単価 \cdot OOH広告 \left(1 + \frac{\beta_{TVCM}}{\beta_{OOH}} \right) = 予算制約
$$

よって、OOH広告の最適水準は：

$$
OOH広告^* = \frac{\beta_{OOH}}{\beta_{OOH} + \beta_{TVCM}} \cdot \frac{予算制約}{OOH広告単価}
$$

同様に、TVCMの最適水準は：

$$
TVCM^* = \frac{\beta_{TVCM}}{\beta_{OOH} + \beta_{TVCM}} \cdot \frac{予算制約}{TVCM単価}
$$

このように、コブ=ダグラス型の関数における**需要関数**は、各広告に割り当てるべき「予算の割合」で決まるのが特徴です。

具体的には、まず広告商材の強さの割合（OOH広告の場合、$\frac{\beta_{OOH}}{\beta_{OOH} + \beta_{TVCM}}$）で予算を決めて（$\frac{\beta_{OOH}}{\beta_{OOH} + \beta_{TVCM}} \cdot 予算制約$）、最後に単価で割って出稿量を決めます。

### 間接効用関数

上記の最適なOOH広告とTVCMの値を元の目的関数（売上関数）に代入することで、最大化された売上、すなわち間接効用関数が得られます：

$$
売上^* = 
\left( \frac{\beta_{OOH}}{\beta_{OOH} + \beta_{TVCM}} \cdot \frac{予算制約}{OOH広告単価} \right)^{\beta_{OOH}} \cdot
\left( \frac{\beta_{TVCM}}{\beta_{OOH} + \beta_{TVCM}} \cdot \frac{予算制約}{TVCM単価} \right)^{\beta_{TVCM}}
$$

この式は、経営会議でよく出てくる、**今の広告予算で売上をどこまで伸ばせるか**を解答する非常に重要な存在です。

また、最後に、強化学習界隈の方にとって、間接効用関数を**価値関数**、需要関数を**方策関数**と理解していただいた方がわかりやすいかもしれません。

では、次にコブ=ダグラス型の関数の売上予測モデルを推定し、予算配分まで実施する方法を紹介します。

# シミュレーション

そもそも、すでに解析的な解が得られているのであれば、コブ＝ダグラス型の売上予測モデルを推定し、その需要関数を使って広告出稿量を決定すればよいのではないか、と思う方もいるかもしれません。

しかし実務では、予測精度の向上や欠落変数によるバイアスの回避を目的として、単純なコブ＝ダグラス型だけではなく、何らかの追加項目を加えたより複雑な構造が用いられるのが一般的です。

たとえば、以下のような拡張がよく見られます：

- 長期トレンドや季節性の影響を除去するためのガウス過程
- 広告の効果が時間とともに減衰する構造を捉えるadstock構造

このような拡張が加わると、モデルは解析的に解くことが難しくなり、最適化や推定には数値的手法が必要となります。

したがって本記事では、まず解析解が既知であるシンプルなコブ＝ダグラス型の売上予測モデルを用いて、数値的に最適化を行い、その結果が理論解と一致するかを確認します。この検証を通じて、実務で扱うより複雑なモデルにも同様の手法が適用可能かどうかを検討します。

まず、学習データを作ります：

```R
set.seed(12345)
ad_df <- 50000 |>
  rnorm() |>
  exp() |>
  tibble::tibble(
    e = _
  ) |>
  dplyr::mutate(
    ooh = abs(rnorm(dplyr::n())),
    tvcm = abs(rnorm(dplyr::n())),
    sales_amount = (ooh ^ 0.3) * (tvcm ^ 0.7) * e
  )
```

$\beta_{OOH}$が0.3、$\beta_{TCVM}$が0.7になります。

続いて、Stanでモデルを記述します：

```cpp:cobb_douglas.stan
data {
  int N;
  
  vector[N] a;
  vector[N] b;
  array[N] real y;
}
parameters {
  simplex[2] beta; // ミクロ経済学では、コブ＝ダグラス型関数のbetaが足して1になると仮定することが多い
  real<lower=0> sigma;
}
model {
  beta ~ dirichlet(rep_vector(0.5, 2));
  sigma ~ inv_gamma(0.1, 0.1);
  log(y) ~ normal(beta[1] * log(a) + beta[2] * log(b), sigma);
}

```

次に、モデルをコンパイルして、パラメーターを変分推論で推定します：

```R
m_cd_init <- cmdstanr::cmdstan_model("cobb_douglas.stan")

m_cd_estimate <- m_cd_init$variational(
  data = list(
    N = nrow(ad_df),
    a = ad_df$ooh,
    b = ad_df$tvcm,
    y = ad_df$sales_amount
  )
)
```

さて、ここではおそらくあまり知られていない、Stanで事後分布のみを抽出するコードを紹介します：

```cpp:cobb_douglas_predict.stan
data {
  real a;
  real b;
}
parameters {
  simplex[2] beta;
  real<lower=0> sigma;
}
generated quantities {
  real predict;
  
  predict = exp(normal_rng(beta[1] * log(a) + beta[2] * log(b), sigma));
}

```

書き方の考えはシンプルです。`cobb_douglas_predict.stan`は、`cobb_douglas.stan`の事後分布のパラメーターを利用するため、`cobb_douglas.stan`と同様の`parameters`を持つ必要があります。また、推定はせずにただ事後分布を出すだけなので、`model`ブロックは不要です。

「`cobb_douglas_predict.stan`をわざわざ分けずに、`cobb_douglas.stan`に組み込めばいいのでは？」と思う方もいらっしゃるかもしれません。しかし、両者を統合してしまうと、予測を行うたびにパラメータの推定からやり直す必要があり、非常に多くの時間がかかってしまいます。

本記事では単純なトイデータを用いていますが、実務では数十万〜数百万レコード規模の学習データを扱うケースも珍しくありません。そのような大規模データで毎回推定をやり直していたら、きっと予測結果が出た頃には、もうTVCMを打つタイミングを逃してしまっていることでしょう、、、、、、

続いて、コンパイルして、予測用の関数とコスト計算用の関数を作ります：

```R
gq_cd_init <- cmdstanr::cmdstan_model("cobb_douglas_predict.stan")

predict_cobb_douglas <- function(S){
  gq_cd_init$generate_quantities(
    fitted_params = m_cd_estimate,
    data = list(
      a = S[1],
      b = S[2]
    )
  )$
    summary() |>
    dplyr::pull(median)
}

compute_cost <- function(S){
  30 * S[1] + 15 * S[2]
}
```

要するに、OOH広告の単価が30円（？）、TVCMの単価が15円の設定になります。

ここで、R言語の`RSlolnp::solnp`関数でラグランジュ未定条数法を実際に数値的に解きましょう！また、広告予算を10円から100円まで変動させて、達成できる最大の売上とその時の広告出稿配分も併せて可視化用に保存します。

```R
future::plan(future::multisession, workers = 16)
opt_results <- seq(10, 100, length.out = 300) |>
  furrr::future_map(
    \(this_cost){
      this_result <- Rsolnp::solnp(
        # 初期値は適当に入れました、、、
        pars = c(10, 10),
        # 最小化を想定しているため、マイナスをつけましょう
        fun = \(S) -predict_cobb_douglas(S),
        eqfun = \(S) compute_cost(S),
        eqB = this_cost,
        # 出稿量がマイナスになってはいけない
        ineqfun = \(S){S},
        ineqLB = c(0, 0),
        ineqUB = c(Inf, Inf),
        control = list(
          outer.iter = 5,
          inner.iter = 5
        )
      )
      
      return(
        list(
          value = this_result$values[length(this_result$values)],
          pars = this_result$pars
        )
      )
    },
    .progress = TRUE
  )
```

かなり時間がかかります。より複雑なモデルを利用する際は、Pythonなど他言語の数値最適化パッケージを利用するか、モデルの近似版を最適化してみてください。モデルの近似版の考え方はまた別の記事で説明する予定です。

では、早速数値的に解いた間接効用関数（与えられた予算のもとで達成できる最大の売上）は、解析的に解いた正解と合致しているかを確認しましょう：

```R
opt_results |> 
  purrr::map_dbl(\(x){-x$value}) |>
  tibble::tibble(
    estimated_value = _
  ) |>
  dplyr::mutate(
    cost = seq(10, 100, length.out = 300),
    theoretical_value = (((0.3 * cost)/30)^0.3) * (((0.7 * cost)/15)^0.7)
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = cost, y = estimated_value), color = "blue", alpha = 0.3) + 
  ggplot2::geom_line(ggplot2::aes(x = cost, y = theoretical_value))

```

![value_function.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/f79249cb-1e46-4d6c-a80e-6e8ae3ba617d.png)

黒い線が解析的に解いた正解で、青い点がR言語の`RSlolnp::solnp`関数で求めた値です。少しずれている箇所もありますが、おおむね正解の直線にフィットしているといえるのではないかと思います。

最後に、TVCMの需要関数（与えられた予算のもとの最適な出稿量）を確認しましょう：

```R
opt_results |> 
  purrr::map_dbl(\(x){x$pars[2]}) |>
  tibble::tibble(
    estimated_value = _
  ) |>
  dplyr::mutate(
    cost = seq(10, 100, length.out = 300),
    theoretical_value = (0.7 * cost)/15
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = cost, y = estimated_value), color = "blue", alpha = 0.3) + 
  ggplot2::geom_line(ggplot2::aes(x = cost, y = theoretical_value))
```

![tvcm_opt.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/3049b16b-c606-4807-a517-d176be6188ac.png)

こちらも、少しずれている箇所がありますが、概ね数学で証明された解析解と同じ形になっています。

# 結論

いかがでしたか？

このように、数値解析の手法を活用することで、推定した売上予測モデルに基づいて、ある程度の妥当性が担保された最適な広告配分を導出することが可能になります。

さらに、間接効用関数を経営会議に提供すれば、マーケティング部門内でのTVCMとOOH広告の配分だけでなく、「そもそも全社としてどの程度の予算をマーケティングに割くべきか」といった、より上位の意思決定を行う際の重要な材料となります。特に上場企業であれば、数十億円規模のマーケティング予算を、より合理的かつ効率的に配分することが可能になるでしょう。

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

https://recruit-dip.com/newgrad/page-1004/