はじめに

ミクロ経済学における基本問題のひとつに「効用最大化」があります。この問題では、消費者は所与の予算制約のもとで効用を最大化すると考えられます。

「支出最小化」は「効用最大化」の双対問題であり、一定の効用水準を達成するために最小限の支出を求めるものです。

この「支出最小化」問題は、例えばマーケティングにおける「広告費を最小化して目標売上や目標ダウンロード数を確保する」問題とアナロジーの関係にあります。本記事ではその数理的対応関係を整理したうえで、RとStanで実装する方法を紹介します。

双対問題

効用最大化

消費者は予算制約のもとで$k$個の財の数量$x=(x_1, \cdots, x_k)$を選び、効用$u(x)$を最大化する。

$$\max_{x\in R^k} u(x)$$

$$s.t. ~\sum_{K=1}^{k}p_{K}\cdot x_{K} = p\cdot x \leq w$$

ただし、$p_K$は財$K$の価格、$w$は予算を表す。
いま、所与の$p,w$のもとで消費者の効用を最大化する需要関数（需要対応）を$x(p,w)\in R^k$で表す。

仮定

ラグランジュ未定乗数法を用いてこの問題を解くために、以下の4つの仮定を置く。

$u(x)$は二回微分可能。
$u(x)$は局所非飽和（locally non-satiated）である。
- いかなる$x\in R^k$と$\epsilon > 0$についても$x^*\in R^k$が存在し、
  $\sqrt{\sum_{K=1}^k(x_K-y_K)^2} < \epsilon$ かつ $u(x^*)>u(x).$
いかなる$p,w, K(=1,\cdots k$)について、$x_K(p,w)>0.$
$x(p,w)$は関数である。
- 「$u(x)$が厳密な準凹関数である。」$\Rightarrow$「$x(p,w)$は関数である。」

ラグランジュ未定乗数法

以下の手順で最適解を求める。

ラグランジアン$\mathfrak{L}=u(x)+\lambda(w-p\cdot x)$を用意する。
それぞれの固定された$\lambda$について、$\mathfrak{L}$を最大化する。
$p\cdot x_{\lambda^{*}}(p,w)=w$を満たすような$\lambda ^ *$を求める。

支出最小化

いま、消費者の行動を別の角度から眺めてみる。

$$\min_{x\in R^k} p \cdot x$$

$$s.t. ~u(x)>u$$

上記の支出最小化問題は、以下のように書き換えてもよい。

$$\max_{x\in R^k} -p \cdot x$$

$$s.t. ~u(x)>u$$

このように、効用を最大化する財の組み合わせは支出を最小化するものであり、逆もまた同様だといえる。

マーケティングへの応用

広告費を投下してダウンロード数を稼ぐことを考え、季節性とトレンドのあるデータを生成する。データは19年1月から24年12月までの月次データとし、コロナウイルスなどの影響は考えない。

ここでは目標のダウンロード数を達成するために必要な最小限の広告費を求める。

set.seed(123) 
data <- tibble::tibble(
  month = seq.Date(from = as.Date("2019-01-01"), 
                   to = as.Date("2024-12-01"), 
                   by = "month"), 
  ad_spend = seq(50, 120, length.out = 72) + rnorm(24, mean = 0, sd = 5), 
  seasonality = rep(c(1.2, 1.1, 1.0, 0.9, 0.8, 1.3, 1.5, 1.4, 1.0, 0.9, 1.1, 1.2), 6), 
  noise = rnorm(72, mean = 0, sd = 300) 
) |>
  dplyr::mutate(
    downloads = round(ad_spend * 100 * seasonality + noise),
    downloads = ifelse(downloads < 0, 0, downloads)
  )

念のため、データを可視化して確認する。

ggplot2::ggplot(data, ggplot2::aes(x = ad_spend, y = downloads)) +
  ggplot2::geom_point(color = "blue") +
  ggplot2::labs(title = "広告費用とアプリダウンロード数",
       x = "広告費用 (万円)",
       y = "アプリダウンロード数")

ggplot2::ggplot(data, ggplot2::aes(x = month, y = downloads)) +
  ggplot2::geom_line(color = "green") +
  ggplot2::labs(title = "月ごとのアプリダウンロード数",
       x = "月",
       y = "アプリダウンロード数")

上が広告費とダウンロード数の相関を示すプロット、下が月次のダウンロード数の推移を示すプロットになっている。

25年1月から6月までの6か月間について、目標とするダウンロード数と投下予定の広告費をデータフレームとして追加する。

prediction <- tibble::tibble(month = seq(from = as.Date("2025-01-01"),
                                                     to = as.Date("2025-06-01"),
                                                     by = "month")) |>
  dplyr::mutate(ad_spend  = c(120, 100, 80, 95, 100, 110),
                seasonality = NA,
                noise = NA,
                downloads = rep(14000, 6)) 

merged <- data |>
  dplyr::bind_rows(prediction)

Stanファイルの実装

季節性とトレンドはガウス過程でモデリングする。

data {
  int time_type;
  array[time_type] real time_seq; 
  
  int N;
  
  vector[time_type] ad;
  vector[time_type] outcome;
}

parameters {
  real<lower=0> rho_gp_exp;                                  // ガウス過程のパラメータ
  real<lower=0> alpha_gp_exp;                                // ガウス過程のパラメータ
  real<lower=0> sigma_gp_exp;                                // ガウス過程のパラメータ
  vector[time_type] eta_gp_exp;
  
  real<lower=0> rho_gp_cycle;                                  // ガウス過程のパラメータ
  real<lower=0> alpha_gp_cycle;                                // ガウス過程のパラメータ
  real<lower=0> sigma_gp_cycle;                                // ガウス過程のパラメータ
  vector[time_type] eta_gp_cycle;

  
  real intercept;
  real<lower=0> sigma;
  
  real<lower=0> rho_global;
  real<lower=0> beta_ad;
}

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_matrix = gp_exp_quad_cov(time_seq, alpha_gp_exp, rho_gp_exp);

    // diagonal elements
    for (i in 1:time_type) {
      K_matrix[i, i] = K_matrix[i, i] + 0.001;
    }

    L_K = cholesky_decompose(K_matrix);
    f_trend = (L_K * eta_gp_exp) - mean(L_K * eta_gp_exp);
  }
  
  {
    matrix[time_type, time_type] L_K;
    matrix[time_type, time_type] K_matrix = gp_periodic_cov(time_seq, alpha_gp_cycle, rho_gp_cycle, 12.0/time_type);

    // diagonal elements
    for (i in 1:time_type) {
      K_matrix[i, i] = K_matrix[i, i] + 0.001;
    }

    L_K = cholesky_decompose(K_matrix);
    f_cycle = (L_K * eta_gp_cycle) - mean(L_K * eta_gp_cycle);
  }
}

model {
  rho_gp_exp ~ inv_gamma(5, 5);
  alpha_gp_exp ~ normal(0, 1);
  sigma_gp_exp ~ normal(0, 1);
  eta_gp_exp ~ normal(0, 1);
    
  rho_gp_cycle ~ inv_gamma(5, 5);
  alpha_gp_cycle ~ normal(0, 1);
  sigma_gp_cycle ~ normal(0, 1);
  eta_gp_cycle ~ normal(0, 1);
  
  rho_global ~ cauchy(0, 1);
  beta_ad ~ double_exponential(0, rho_global);
  
  intercept ~ normal(0, 10);
  sigma ~ inv_gamma(0.01, 0.01);

  log(outcome) ~ normal(
      intercept + 
      f_trend[1:N] + 
      f_cycle[1:N] + 
      beta_ad * log(ad),
      sigma
    );
}

変分推論でパラメータの事後分布を推定する。

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

m_estimate <- m_init$variational(
  data <- list(
    time_type = nrow(merged),
    time_seq = (1:nrow(merged))/nrow(merged) - 0.5,
    N = nrow(merged),
    ad = merged$ad_spend,
    outcome = merged$downloads
  )
)

上のstanファイルに規定されたパラメータを使って25年1月～6月のダウンロード数の予測値が分かればよく、この期間のパラメータは必要ない。

そこで、genrated quantitiesブロックをstand aloneさせる。

data {
  int time_type;
  array[time_type] real time_seq; 
  
  int new_N;
  array[new_N] int new_time; 
  array[new_N] real ad_new;
}

parameters {
  real<lower=0> rho_gp_exp;
  real<lower=0> alpha_gp_exp;
  real<lower=0> sigma_gp_exp;
  vector[time_type] eta_gp_exp;
  
  real<lower=0> rho_gp_cycle;
  real<lower=0> alpha_gp_cycle;
  real<lower=0> sigma_gp_cycle;
  vector[time_type] eta_gp_cycle;

  real intercept;
  real<lower=0> sigma;
  
  real<lower=0> rho_global;
  real<lower=0> beta_ad;
}

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_matrix = gp_exp_quad_cov(time_seq, alpha_gp_exp, rho_gp_exp);

    // diagonal elements
    for (i in 1:time_type) {
      K_matrix[i, i] = K_matrix[i, i] + 0.001;
    }

    L_K = cholesky_decompose(K_matrix);
    f_trend = (L_K * eta_gp_exp) - mean(L_K * eta_gp_exp);
  }

  {
    matrix[time_type, time_type] L_K;
    matrix[time_type, time_type] K_matrix = gp_periodic_cov(time_seq, alpha_gp_cycle, rho_gp_cycle, 12.0/time_type);

    // diagonal elements
    for (i in 1:time_type) {
      K_matrix[i, i] = K_matrix[i, i] + 0.001;
    }

    L_K = cholesky_decompose(K_matrix);
    f_cycle = (L_K * eta_gp_cycle) - mean(L_K * eta_gp_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_ad * log(ad_new[i]),
        sigma));
        }
}

このstand-alone generated quantitiesなるStanファイルをコンパイルする。

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


gq_estimate <- gq_init$generate_quantities(
  fitted_params = m_estimate,
  data = list(
    time_type = nrow(merged),
    time_seq = (1:nrow(merged))/nrow(merged) - 0.5,
    new_N =6,
    new_time = (nrow(merged)-6+1):nrow(merged), 
    ad_new = c(120, 100, 80, 95, 100, 110)    
  )
)

支出最小化

いま考えているのは、以下の支出最小化問題である。

$$\min_{S} \sum_{t=1}^6S_t$$

$$\mu_t(S)\geq A_t, ~ t = T-5, \cdots, T$$

$$S_t\geq0, \forall t$$

ただし、$S_t$は$t$月の広告費。$A_t$は目標とする$t$月のダウンロード数。

R言語のRSlolnp::solnpでこの問題を解く。

targets <- c(18000, 17000, 14000, 14000, 15000, 16000)  

predicted_means <- function(S){
  out_summary <- gq_init$generate_quantities(
    fitted_params = m_estimate,
    data = list(
      time_type = nrow(merged),
      time_seq = (1:nrow(merged))/nrow(merged) - 0.5,
      new_N = 6,
      new_time = (nrow(merged) - 6 + 1):nrow(merged),
      ad_new = S[1:6]
    )
  )$summary()

predicted <- out_summary |> 
  dplyr::pull(mean)
  return(as.numeric(predicted)) 
}

obj_min_spend <- function(S){
  sum(S)
}

ineqfun_means <- function(S){
  predicted_means(S)
}

start_guess <- c(120, 100, 80, 95, 100, 110)

res <- Rsolnp::solnp(
  pars = start_guess,
  fun = obj_min_spend,          
  ineqfun = ineqfun_means,
  ineqLB = targets,        
  ineqUB = rep(Inf, 6),
  LB = rep(0, 6)                 
)

25年1月～6月までの6か月間の目標ダウンロード数(18000, 17000, 14000, 14000, 15000, 16000) を達成するために必要な最小の広告投資は(109.88237,102.81089,90.77335, 94.90221, 94.63948, 115.08639)
となった。

注意：今回は係数×広告費にガウス過程でモデリングしたトレンドや季節性を加えることで、モデルを拡張している。これによって解析的な解を得るのが難しくなっている可能性がある（つまり、得られた解が最適解とは限らないということ）。その場合は数値的手法が必要となる。