まずは蝋の翼から。

学んだことを書きながら確認・整理するためのメモブログ。こういうことなのかな?といったことをふわっと書いたりしていますが、理解が浅いゆえに的はずれなことも多々あると思うのでツッコミ歓迎

資源の最適な配分を、ディリクレ分布を用いたシミュレーションから求める

やりたいこと

なにかしらの資源をどう分配したら効用を最大化できるか、ということを解析ではなくシミュレーションによって求める。

例題

今回、例として「総予算上限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%じゃない限り誰も接触していないことになるのでそれは直感とは反していることになるので少し不適切な気もします。