08-simrun-random-draw.Rmd

# random draw

## I/O

- in: `df_gene` etc.
- out: `sim_random_df` `df_series_gen`: それぞれのシリーズのgenつきのdf_series, etc.

## summary

ここでの目的は、株式会社ポケモンのデザイナーの思考パターンをランダムでどこまで模せるかを調べること。そのためにはいくつかのpresumptionが必要になる。まず、あるタイムステップ（具体的には`release_date`のどこか）`t = t_n`で、すでにリリースされているポケモンのゲーム・すでに公開されているポケモンのリストがあるとする。`t_n`においてリリースする枚数は興味がないので、すでに決定しているとする（`release_date`における`total_cards`に相当する）。そのとき、デザイナーは次のように適当にポケモンを選ぶことが考えられる：

1. すでに発表されているポケモンのプールの中からランダムに枚数ぶん選ぶ。重複を許す。
1. すでに発表されているポケモンのプールの中からランダムに枚数ぶん選ぶ。重複をなるべく避ける（その回にすでにひいたポケモンにネガティブなバイアスがかかる）。毎回のリリースで偏りが生じる。
1. すでに発表されているポケモンのプールの中からランダムに枚数ぶん選ぶ。重複を好む（すでにひいたポケモンにポジティブなバイアスがかかる）
1. すでに発表されているポケモンのプールの中からランダムに枚数ぶん選ぶのだが、そのとき`t_n`が最も最近（頭痛が痛い…）の`game_release_date`に近い場合、新しく追加されたポケモンを好んで選ぶ（新しいポケモンの発表直後はどうも多くリリースされているように見えるので）。pro-noveltyな変異 bias.
1. すでに発表されているポケモンのプールの中からランダムに枚数ぶん選ぶのだが、Generation Iの151匹から好んで選ぶ。第1世代へポジティブなバイアス。やっぱり人気だから…イーブイとかピカチュウとか…。
1. すでに発表されているポケモンのプールの中からランダムに枚数ぶん選ぶするのだが、今までのリリースで多く選ばれたものを好んで選ぶ。正の頻度依存変異バイアス。
1. すでに発表されているポケモンのプールの中からランダムに枚数ぶん選ぶするのだが、前回のリリースで選ばれたポケモンを好んで選ぶ。伝達バイアス。

これだけあればとりあえず十分かと・・・。まず１だけ実装してみよう。
nested dfでやるかmatrixでやるかはたまたdfでやるか。毎回population_sizeが違うのでnested dfがいいのかな



```{r}
series_gen_df3 <- series_gen_df2 |> 
  left_join(df_game_gen |> select(game_gen, pokemon_cumulative_n), by = "game_gen") |> 
  filter(!is.na(pokemon_data_count))
# series_gen_df3
```

## actual gen of releases

実際には新しいgameがリリースされてすぐにそのGenのポケカがでてくるわけではない（たとえば第二世代10/15の翌日のポケカリリースGym Challenge 2000-10-16は第一世代のものばかり）。それにのっとったシミュレーションをまわさないとtotal debutantsのデータがおかしくなる。

まずは実データから、最大どのPokedex_idのカードがでているのかをみる：

```{r}
df_gene |> 
  filter(is_pokemon) |> 
  # filter(n > 40) |> # 一度に４０枚以上リリースされているもののみ
  # filter(pokedex_id == max(pokedex_id)) |> 
  ggplot(aes(release_date, pokedex_id)) +
  geom_point(size = .3) +
   geom_text(
     data = df_gene |> 
       filter(is_pokemon) |> 
       group_by(release_date) |> 
       filter(pokedex_id == max(pokedex_id)),
    aes(label = series, y = pokedex_id),
    size = 2 * pt_convert,
    colour = pokemon_blue,
    hjust = 0,
    # angle = 90
  ) +
  geom_point(
    data = df_game_gen,
    aes(game_release_date, pokemon_cumulative_n),
    colour = "red",
    size = .2
  ) +
 
  geom_hline(
    data = df_game_gen,
    aes(yintercept = pokemon_cumulative_n),
    colour = "red",
    size = .2
  ) +
  # scale_x_date(
  #   limits = c(as.Date("2016-01-01"), NA)
  # ) +
  # scale_y_continuous(
  #   limits = c(600, NA)
  # ) +
  coord_flip() +
  theme_pokemon 
ggsave("./output/max_pokedex_id.svg", width = 200, height = 100, units = "mm")
```
やっぱりだいぶありますね。２０１７年？くらいのEvolutionsなどは第１世代しかいない。cumulative nを頼りに、たとえばmax pokedex idが251以下151以上であればGen II（までのポケモンからはとれたはずだよね）とみなす、というのがよさそう。

そのデータを作ろう：
このやり方だとEvolutionsがactual_release_genが１になってしまい、まあそうなんだけど、それこそがバイアスのある作り方・drawの仕方なので、、、、、
```{r}
series_gen_actual_df1_too_accurate <-
# df_view <- 
  df_gene |> 
  filter(is_pokemon) |> 
  group_by(release_date) |> 
  filter(pokedex_id == max(pokedex_id)) |> 
  mutate(
    # mapを使えばもう少しスマートに書けそうだが、８個ならいいよね？
    # 最大のpokedex_idがある世代の合計のポケモン数以下であればその世代を割り振る
    actual_release_gen = case_when( 
      pokedex_id <= df_game_gen$pokemon_cumulative_n[1] ~ 1,
      pokedex_id <= df_game_gen$pokemon_cumulative_n[2] ~ 2,
      pokedex_id <= df_game_gen$pokemon_cumulative_n[3] ~ 3,
      pokedex_id <= df_game_gen$pokemon_cumulative_n[4] ~ 4,
      pokedex_id <= df_game_gen$pokemon_cumulative_n[5] ~ 5,
      pokedex_id <= df_game_gen$pokemon_cumulative_n[6] ~ 6,
      release_date <= df_game_gen$game_release_date[7] ~ 6,
      pokedex_id <= df_game_gen$pokemon_cumulative_n[7] ~ 7, 
      # マギアナのいるSteam Siegeは~6を満たさず~7を満たしてしまう。~6を満たすようにする
      pokedex_id <= df_game_gen$pokemon_cumulative_n[8] ~ 8,
    )
  ) |> 
  select(series, release_date, actual_release_gen) |> 
  distinct()

  series_gen_df2 |> 
  left_join(df_game_gen |> select(game_gen, pokemon_cumulative_n), by = "game_gen") |> 
  filter(!is.na(pokemon_data_count)) 
```
→ゲームがリリースされ、Eraが更新されたとしても、ポケカは常に第一世代からその最新世代までのポケモンをもとに作られるわけではないことがわかる。しかし、たとえばEvolutionsとかもすでにたくさんの新しい世代があるのに、それを無視して古い世代をやろうとしたわけで、まさに変異のバイアスといえる。そのバイアスまでランダムモデルに反映させてしまうのはよくない。あくまで反映させるべきは「ゲームがでて翌日に新しいポケカがだせるわけではない」点なので、新しい世代のポケモンが登場したポケカシリーズ以降はその世代からランダムどろーできるものとする、とするのがよいだろう。つまり、
```{r}
df_series_first_release_date <-
  series_gen_actual_df1 |> 
  ungroup() |> 
  group_by(actual_release_gen) |> 
  filter(release_date == min(release_date)) |> 
  rename(series_gen = actual_release_gen, first_series = series) 
  
```
この日付以降。ただし、このdfには２つの問題がある。

1. まず、あまりに実態にあわせている。Evolutionsはgen 1ということになっていて、まあ正しいといえば正しいのだがRandom drawのNull modelにはあわない。
1. これが`df_game_gen$release_date`よりも早ければ、そちらを優先する。なぜなら、Steam Siegeに関してはマギアナが例外的に第７世代登場以前に出ており、ここを基準にすると大変なことになるので、変える必要がある。

これらを勘案して、次のように調整する。それぞれのシリーズにGenを付与したものを`df_series_gen`と、最初のリリースを`df_series_first_release_date_corrected`としよう。
```{r}
df_series_gen <-
  series_gen_df2 |> 
  mutate(
    series_gen = case_when(
      release_date < df_release_gen$release_date[2] ~ 1,
      release_date < df_release_gen$release_date[3] ~ 2,
      release_date < df_release_gen$release_date[4] ~ 3,
      release_date < df_release_gen$release_date[5] ~ 4,
      release_date < df_release_gen$release_date[6] ~ 5,
      release_date <= df_release_gen$release_date[7] ~ 6, # magiana taisaku de <=
      release_date < df_release_gen$release_date[8] ~ 7,
      TRUE ~ 8,
    )
  ) |> 
  left_join(
    df_game_gen |> select(game_gen, pokemon_cumulative_n), 
    by = c("series_gen" = "game_gen")) |> 
  filter(!is.na(pokemon_data_count)) |> 
  select(series, series_gen, release_date, everything()) |> 
  arrange(release_date)

df_series_first_release_date_corrected <- df_series_gen |> 
  group_by(series_gen) |> 
  rename(first_series = series) |> 
  filter(release_date == min(release_date))
```
Steam Siegeのseries_genが6になっていることをかくにん！よかった。その直後があのEvolutionsなので困りますがまあこれで正しい。エッジケースの対処としては非常にマニュアルで、あまり良くないが時間がないので仕方ない。


## run


[Unnest a list of vectors in a data frame in R](https://stackoverflow.com/questions/66819613/unnest-a-list-of-vectors-in-a-data-frame-in-r)
```{r}
sim_conditions_df <- df_series_gen |> 
  select(
    series_gen, series, release_date, game_gen,
    pokemon_data_count, pokemon_cumulative_n
    )

simulation_runs <- 1000

for (simulation_run in 1:simulation_runs) {
  single_sim_result_df <- sim_conditions_df |> 
    mutate(
      pokedex_id =  map2( # rowwise() does not work...?
        .x = pokemon_data_count,
        .y = pokemon_cumulative_n,
        ~ runif(
          n = .x,
          min = 1,
          max = .y + 1
        ) |> floor()
      ),
      run = simulation_run
    ) |> 
    unnest_longer(pokedex_id)
  
  # initialize before bind_rows
  if (simulation_run == 1) {
    sim_random_df <- tibble(
      series = character(),
      release_date = Date(),
      game_gen = integer(),
      series_gen = integer(),
      pokemon_cumulative_n = integer(),
      pokedex_id = integer(), 
      run = integer()
    )
  }
  sim_random_df <- bind_rows(sim_random_df, single_sim_result_df)
} #  7min? for 1000 runs
sim_random_df |> dim() # 1 2325 000 x 8
```


## compare gen histogram 

### count
```{r}
df_empirical_count_id <- 
  df_gene |> 
  filter(is_pokemon) |> # filter only pokemons. is_
  count(pokedex_id) |> 
  complete(pokedex_id = 1:905, fill = list(n = 0)) |> 
  mutate(
    condition = "empirical"
  )
sim_random_count_pokemon_df <- 
  sim_random_df |> 
  # left_join(df_names, by = "id") |> 
  count(pokedex_id, run) |> 
  complete(
    pokedex_id = 1:905,
    run = 1:simulation_runs, 
    fill = list(n = 0)
  ) |> 
  mutate(
    condition = "simulation"
  ) 


```
### merge

```{r}
emp_sim_count_id_df <- 
  bind_rows(
    df_empirical_count_id, 
    sim_random_count_pokemon_df
  ) |> 
  left_join(df_names, by = "pokedex_id") 
  
# emp_sim_count_id_df |> colnames()
```

### plot the sim~emp difference

SDを準備（重いし、大量に書かれるので）
```{r}
emp_sim_count_id_summary_df <- emp_sim_count_id_df |> 
  group_by(condition, pokemon_gen) |> 
  summarise(
    mean = mean(n),
    sd = sd(n)
  )
```
annotation用のdfを用意。上位３ポケモンの画像を貼る
```{r}
emp_sim_count_id_annotate_df <-
  emp_sim_count_id_df |>
  filter(condition == "empirical") |> 
  group_by(pokemon_gen) |> 
  mutate(rank = rank(-n, ties.method = "min")) |> 
  filter(rank <= 3 | n > 41 | rank == max(rank))
  
```
これ↓を掲載予定
geom_quasirandomをつかっていると非常にプロットが重い（最終的に70MBとかになる）。実データに関してはviolinを使うか、pngで書き出して扱ったほうがよさそう
```{r}
aspect_ratio <- 1/1.5
p_emp_sim_count_id <-
  emp_sim_count_id_df |> 
  ggplot(aes(x = fct_rev(as.factor(pokemon_gen)), y = n, group = pokemon_gen, fill = condition)) +
  geom_violin(
    data = emp_sim_count_id_df |>
      filter(condition == "simulation"),
    colour = "transparent",
    scale = "width",
    bw = .5,
  ) +
  geom_quasirandom(
    data = emp_sim_count_id_df |> filter(condition == "empirical"),
    aes(group = pokemon_gen, fill = condition),
    pch = 21,
    size = .7,
    # dodge.width = .8,
    stroke = 0,
    colour = "transparent"
    ) +
geom_pointrange(
    data = emp_sim_count_id_summary_df,
    aes(
      y = mean,
      ymin = mean-sd,
      ymax = mean+sd,
      group = condition,
      colour = condition
    ),
    position = position_dodge(.8),
    fatten = .1,
    stroke = 1,
    size = .2
  ) +
  geom_image(
    data = emp_sim_count_id_annotate_df,
    aes(image = image_large),
    size = .06,
    height = .05,
    by = "width",
    asp = 1/aspect_ratio,
    nudge_x = .4,
  ) +
  scale_colour_manual(
    values= c("simulation" = pokemon_yellow_stronger, "empirical" = pokemon_blue)
  ) +
  scale_fill_manual(
    values= c(
      "simulation" = pokemon_yellow |> lighten(.5), 
      "empirical" = pokemon_blue |> darken(.4)
    )
  ) +
  scale_x_discrete(
    breaks = 1:8,
    labels = 1:8 |> as.roman(),
    expand = c(.1, .1),
    name = "Pokémon Generation"
  ) +
  scale_y_continuous(
    limits = c(1, 80), # ignore pikachu!
    breaks = c(1, 40, 80),
    expand = c(0, 0),
    name = "Freqeuncy of Pokémon Cards",
    # trans = "log10"
  ) +
  coord_flip() +
  theme_pokemon +
  theme(
    legend.position = c(.8, .7),
    aspect.ratio = aspect_ratio
  )
ggsave("./output/p_emp_sim_count_id.png", width = 1000, height = 1000, unit = "px", dpi = 600)
ggsave("./output/p_emp_sim_count_id.svg", width = 83, height = 45, unit = "mm", dpi = 600) 


```
`geom_quasirandom()` points might not be visible in rstudio, but works fine for svgs.

Caption:
各世代に属するポケモンの、2022年9月までに発売されたカードの数。オレンジが実際のデータ、ダークグレーがシミュレーション結果。白点と線は平均±SD。第一世代に属するPikachuに関しては、実データにおいてn=135で全世代通してトップだが、あまりに大きいためプロットからは除いた（計算にはもちろん含めてある）。

というかEeveeが2位だと思ってたのだがCharizardなのだね。pokemon_nameとcard_nameでだいぶ違うのだろう。Charizardという名前card_nameでは少ないけど、Charizardというポケモンpokemon_nameでは多いということ。つまりEeveeよりもリミックスされがちだということだ！
リミックス度が測れるな　横軸がポケモンの数、縦軸がVariationの数的な


#### variations (ignore)

```{r}
sim_random_count_pokemon_df |> 
  filter(run < 10) |> 
  ggplot(aes(x = pokedex_id, y = n)) +
  geom_jitter(alpha = .1, size = .2, colour = "grey40") +
  geom_point(data = df_empirical_count_id, size = 1, alpha = .5, colour = "orange") +
  theme_pokemon
ggsave("./output/count_emp_vs_sim_total_individual_ID.png", width = 83, height = 35, unit = "mm", dpi = 600)

```