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

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
  )

data {
  int N;
  
  vector[N] a;
  vector[N] b;
  array[N] real y;
}
parameters {
  simplex[2] beta; // ミクロ経済学では、コブ＝ダグラス型関数のbetaが足して1になると仮定することが多い
  reallower=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);
}

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
  )
)

data {
  real a;
  real b;
}
parameters {
  simplex[2] beta;
  reallower=0> sigma;
}
generated quantities {
  real predict;
  
  predict = exp(normal_rng(beta[1] * log(a) + beta[2] * log(b), sigma));
}

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]
}

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
  )

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))

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))