資源の最適な配分を、ディリクレ分布を用いたシミュレーションから求める
やりたいこと
なにかしらの資源をどう分配したら効用を最大化できるか、ということを解析ではなくシミュレーションによって求める。
例題
今回、例として「総予算上限10万円分で各時間帯にCMを打つ。そのとき、3回以上CMに接触した人数が多いCM投下パターンを求めたい。」という例題で考える
条件
CMを打つ時間によって広告費が変わる
設定
テキトーに以下のような値段設定にする(なんとなく人が多そうな時間ほど高くした)。
時間 | 広告費 |
---|---|
0時 | 1000 |
1時 | 800 |
2時 | 600 |
3時 | 400 |
4時 | 200 |
5時 | 600 |
6時 | 900 |
7時 | 1500 |
8時 | 2000 |
9時 | 1500 |
10時 | 1500 |
11時 | 1500 |
12時 | 2000 |
13時 | 1500 |
14時 | 1500 |
15時 | 1500 |
16時 | 1500 |
17時 | 1500 |
18時 | 3000 |
19時 | 5000 |
20時 | 4500 |
21時 | 4500 |
22時 | 3000 |
23時 | 2000 |
実装
CM投下パターン
各時間帯にどれくらいの金額を使うか、ランダムに配分をおこなう。
その際、ディリクレ分布を用いる。ディリクレ分布とは簡単にいうと「各面の出やすさをパラメータで設定した多面サイコロを振ったときの各面の出る確率」を出すことができる。「各面の出やすさが同じなら確率は等しいのでは?」と思うし、実際無限回行えばそうなるが乱数生成する場合は回数に制限があるため各面の確率はゆらぎがある。
以下に5種類の3面サイコロ結果を示す。
rdirichlet(5, rep(1, 3)) [,1] [,2] [,3] [1,] 0.44650941 0.3393241 0.2141665 [2,] 0.23953497 0.5266908 0.2337743 [3,] 0.05629639 0.4216316 0.5220720 [4,] 0.03596511 0.1777634 0.7862715 [5,] 0.40526176 0.3653088 0.2294295
確率なので各面を足し合わせると1.0になる。
ディリクレ分布のパラメータをサイコロに例える | コード7区
今回、(サイコロで例えると)どの面も出る確率が等しい24面サイコロを振ることで各面、つまり各時間帯が選ばれる確率を生成できる。
前述のように、確率なので各面を足し合わせると1.0になる。考え方を変えると、これは1.0という数字を24つに対してランダムに配分していることになる。そのため、この確率に今回使用する予定の10万円をかけると各時間に10万円をランダムに配分した結果が得られる。
このランダム配分を1000回シミュレーションをする。
library(gtools) all_money = 10^5 # 使用する金額 n_simulation = 10^3 # シミュレーション回数 # 各シミュレーションでの、時間毎使用料金 use_simulation = tibble(simulation_no = rep(seq_len(n_simulation), each = 24), # id=セット数 hour = rep(0:23, n_simulation), money_rate = rdirichlet(n_simulation, rep(1, 24)) %>% t() %>% matrix(ncol = 1, byrow = T) %>% c()) %>% mutate(use_money = money_rate * all_money) # # A tibble: 240 x 4 # simulation_no hour money_rate use_money # <int> <int> <dbl> <dbl> # 1 1 0 0.0432 4322. # 2 1 1 0.103 10302. # 3 1 2 0.0165 1646. # 4 1 3 0.0251 2508. # 5 1 4 0.0472 4724. # 6 1 5 0.00851 851. # 7 1 6 0.0316 3157. # 8 1 7 0.156 15562. # 9 1 8 0.0373 3728. # 10 1 9 0.171 17112. # 11 1 10 0.0531 5306. # 12 1 11 0.00363 363. # 13 1 12 0.0112 1121. # 14 1 13 0.0884 8845. # 15 1 14 0.0184 1844. # 16 1 15 0.00720 720. # 17 1 16 0.00863 863. # 18 1 17 0.00893 893. # 19 1 18 0.000429 42.9 # 20 1 19 0.0668 6685. # 21 1 20 0.0152 1515. # 22 1 21 0.00644 644. # 23 1 22 0.0546 5456. # 24 1 23 0.0179 1790. # 25 2 0 0.0563 5627. # # … with 215 more rows
各時間のCM放送回数
時間毎CMの料金を先程求めた、時間毎使用料金にJOINする。各時間に使用できる料金(予算上限)は先程求めたので1回あたりのCM料金で割ることで、各時間で何回CMを放送するか決まる。
# 時間毎CM料金 hourly_price <- read_csv('./cm_price.csv') # # A tibble: 24 x 2 # hour price # <dbl> <dbl> # 1 0 1000 # 2 1 800 # 3 2 600 # 4 3 400 # 5 4 200 # 6 5 600 # 7 6 900 # 8 7 1500 # 9 8 2000 # 10 9 1500 # # … with 14 more rows simulation_hourly_cm_cnt = use_simulation %>% inner_join(hourly_price) %>% # 各時間の1回あたりCM料金 rename(hourly_price = price) %>% mutate(cm_cnt = floor(use_money/hourly_price)) # 放送可能回数 # # A tibble: 240 x 6 # simulation_no hour money_rate use_money hourly_price cm_cnt # <int> <dbl> <dbl> <dbl> <dbl> <dbl> # 1 1 0 0.0442 4423. 1000 4 # 2 1 1 0.111 11096. 800 13 # 3 1 2 0.0881 8806. 600 14 # 4 1 3 0.0497 4968. 400 12 # 5 1 4 0.0252 2524. 200 12 # 6 1 5 0.0932 9319. 600 15 # 7 1 6 0.0253 2525. 900 2 # 8 1 7 0.0113 1127. 1500 0 # 9 1 8 0.0463 4627. 2000 2 # 10 1 9 0.0238 2383. 1500 1 # # … with 230 more rows
時間毎の人の有無
今回はテキトーに乱数で作成するが、実務においては実際のユーザー行動からユーザー毎に各時間帯の滞在確率を求める。
観察できるユーザー人数は1000として、各ユーザーの各時間帯の滞在率を乱数で作成。
library(tidyverse) n_user = 1000 # 対象ユーザー数 # ユーザー毎の、各時間の滞在確率 user_prob = tibble(user_id = rep(seq_len(n_user), each = 24), # ユーザーIDを24時間分 hour = rep(0:23, n_user), #時間 prob = rep(sample(0:100, 24 * n_user, replace=TRUE)/100) # 滞在確率 ) # # A tibble: 24,000 x 3 # user_id hour prob # <int> <int> <dbl> # 1 1 0 0.84 # 2 1 1 0.1 # 3 1 2 0.82 # 4 1 3 0 # 5 1 4 0.14 # 6 1 5 0.51 # 7 1 6 0.69 # 8 1 7 0.55 # 9 1 8 0 # 10 1 9 0.75 # # … with 23,990 more rows
接触回数
各時間毎のCM放送回数と、各時間毎のユーザー毎滞在確率が求まった。そのため、これらをかけ合わせると各時間毎のユーザー毎滞在期待値が求まる。
このとき、期待値が1以上の人は1回以上接触するユーザーと考えることができるし、期待値が3以上の人は3回以上接触するユーザーと考えることができる。
simulation_resulut = simulation_hourly_cm_cnt %>% inner_join(user_prob) %>% mutate(exp_reach = cm_cnt * prob) %>% mutate(reach_1over = if_else(exp_reach >= 1.0,1,0), reach_3over = if_else(exp_reach >= 3.0,1,0)) %>% group_by(simulation_no) %>% summarise(reach_1over_cnt = sum(reach_1over), reach_3over_cnt = sum(reach_3over)) %>% arrange(-reach_3over_cnt) # # A tibble: 1,000 x 3 # simulation_no reach_1over_cnt reach_3over_cnt # <int> <dbl> <dbl> # 1 226 9106 5781 # 2 721 10618 5549 # 3 882 9026 5535 # 4 290 9649 5502 # 5 491 9694 5472 # 6 229 10332 5438 # 7 470 9940 5376 # 8 640 8784 5356 # 9 180 8813 5332 # 10 78 8286 5322 # # … with 990 more rows
3回以上接触者が多い投下パターン
上記結果より、simulation_no 226において3回以上接触者が5781を獲得していて、1000回シミュレーションをした結果では一番多い。
このときの10万円の使用配分(CM投下パターン)は以下のようになる。
simulation_hourly_cm_cnt %>% filter(simulation_no == 226) %>% print(n=24) # # A tibble: 24 x 6 # simulation_no hour money_rate use_money hourly_price cm_cnt # <int> <dbl> <dbl> <dbl> <dbl> <dbl> # 1 226 0 0.107 10677. 1000 10 # 2 226 1 0.0695 6947. 800 8 # 3 226 2 0.0597 5966. 600 9 # 4 226 3 0.0270 2698. 400 6 # 5 226 4 0.0410 4096. 200 20 # 6 226 5 0.152 15150. 600 25 # 7 226 6 0.0595 5951. 900 6 # 8 226 7 0.00883 883. 1500 0 # 9 226 8 0.00989 989. 2000 0 # 10 226 9 0.0182 1816. 1500 1 # 11 226 10 0.0163 1630. 1500 1 # 12 226 11 0.0534 5337. 1500 3 # 13 226 12 0.0289 2895. 2000 1 # 14 226 13 0.00924 924. 1500 0 # 15 226 14 0.00603 603. 1500 0 # 16 226 15 0.130 12969. 1500 8 # 17 226 16 0.0312 3125. 1500 2 # 18 226 17 0.0810 8097. 1500 5 # 19 226 18 0.0257 2566. 3000 0 # 20 226 19 0.00753 753. 5000 0 # 21 226 20 0.0184 1842. 4500 0 # 22 226 21 0.0295 2946. 4500 0 # 23 226 22 0.00714 714. 3000 0 # 24 226 23 0.00427 427. 2000 0
おまけ
実際のCM投下はここまで単純ではないです。
また、時間帯ごとCM回数は切り捨てにしているので各シミュレーションで全使用料金が変わる問題があります。
さらに言えば、そもそも投下回数が1回の時間帯では、接触期待値部分は1 * 確率(1以下) ≦ 1 となるので100%じゃない限り誰も接触していないことになるのでそれは直感とは反していることになるので少し不適切な気もします。