09-CCDF.Rmd

# emp pokemon_name vs sim pokemon_name, df_game_gen$game_release_date

## I/O

- input: `sim_random_df`, `df_gene`
- output: `temporal_change_per_gen_df`
- for: 10: p_pandoras_aggregate
## summary

1. 第一世代までのリリースのCCDF、第二世代までのリリースのCCDF、…と累積的に、シミュレーションとデータ描いていく。
1. 第一世代のリリースのCCDF、第二世代の時期のリリースのCCDF、…と経時的に、シミュレーションとデータ描いていく。
どういう差がでますかね？
まずは累積的に


## cumulative

filter dfs until specific year (~2000, ~2014, etc.), merge them


```{r}
sim_random_df # 1000 runs
df_gene # `card_name` = group by distinct names, `pokemon_name` = group by distinct pokemons. 
# initialize
temporal_change_cumulative_df <- tibble(
  n = integer(),
  ecdf = double(),
  run = integer(),
  until_date = Date()
)
# run
# i starts from 2 to skip "until 1999-01-01" which doesn't exist
# i ends with 2022-11-18, which will hold all time data.
for (i in 2:length(c_series_gen_dates)) { 
  sim_random_df_until_threshold_date <-
    sim_random_df |> 
    filter(release_date <= c_series_gen_dates[i])
  emp_df_until_threshold_date <- df_gene |> 
    filter(is_pokemon) |> 
    filter(release_date <= c_series_gen_dates[i])
  
  sim_df <-
    sim_random_df_until_threshold_date |> 
    group_by(run, pokedex_id) |> 
    mutate(n = n()) |> 
    ungroup() |> 
    group_by(run) |> 
    arrange(desc(n)) |>
    mutate(ecdf = 1 - ecdf(n)(n - .01)) |> 
    select(run, n, ecdf) |> 
    distinct() |> 
    ungroup() %>%
    mutate(
      condition = rep("simulations", nrow(.))
    )
  emp_card_type_df <- 
    emp_df_until_threshold_date |> 
    mutate(name_type = paste(card_name, card_type)) |> 
    group_by(name_type) |> 
    summarise(n = n()) |> 
    arrange(desc(n)) |> 
    mutate(ecdf = 1 - ecdf(n)(n - .01)) |>
    select(-name_type) |> 
    distinct() %>%
    mutate(condition = rep("card name + type", nrow(.)))
  emp_card_df <- 
    emp_df_until_threshold_date |> 
    group_by(card_name) |> 
    summarise(n = n()) |> 
    arrange(desc(n)) |> 
    mutate(ecdf = 1 - ecdf(n)(n - .01)) |>
    select(-card_name) |> 
    distinct() %>%
    mutate(condition = rep("card name", nrow(.)))
  emp_poke_df <- 
    emp_df_until_threshold_date |> 
    group_by(pokedex_id) |> 
    summarise(n = n()) |> 
    arrange(desc(n)) |> 
    mutate(ecdf = 1 - ecdf(n)(n - .01)) |>
    select(-pokedex_id) |> 
    distinct() %>%
    mutate(condition = rep("pokemon name", nrow(.)))
  merge_df <- bind_rows(
    sim_df, 
    emp_card_type_df,
    emp_card_df,
    emp_poke_df
  ) |> 
    mutate(
      until_date = c_series_gen_dates[i]
    )
  temporal_change_cumulative_df <- bind_rows(temporal_change_cumulative_df, merge_df)
} # 30 sec
temporal_change_cumulative_df |> dim() # 159039 ~~187990 x 5~~

```


#### plot ccdf dataVSsim, game gen, cumulative

論文掲載予定だったけどとりあえずヤメ
compare card name+type name~card name~pokemon name~simulation, ~ cumulative, ~game gen



```{r fig.width = 3, fig.height = 3}
p_temporal_change_cumulative <-
  temporal_change_cumulative_df |> 
  filter(condition == "simulations") |> 
  # filter(until_date %in% five_years_seq) |> 
  ggplot(aes(x = n,y = ecdf, group = paste(until_date, run))) +
  # geom_point(
  #   # aes(alpha = condition, shape = condition),
  #   size = .1
  # ) +
  geom_line(size = .3, alpha = .4, aes(colour = run)) + 
  geom_point(
    data = temporal_change_cumulative_df |> 
  filter(condition == "pokemon name"),
    colour = "#3561B9", size = 1
  ) +
  # scale_alpha_manual(values = c(1,1, 1, .2)) +
  # scale_shape_manual(values = c(16, 16,16, 21)) +
  scale_x_continuous(
     expand = c(0.1,0),
     limits = c(1, NA),
    trans = "log10",
    breaks = 10^(0:10),
    name = "Frequency of Pokémon"
  ) +
  scale_y_continuous(
    expand = c(0.1,0),
    name = expression(italic("Pr") ( X>= x) ),
    # trans = "log10",
    breaks = 10^c(0, -3),
    labels = trans_format("log10", math_format(10^.x))
    ) +
  # scale_colour_manual(values = c("black", "grey50")) +
  # scale_size_manual(values = c(.5, .3/.751)) +
   scale_colour_viridis_c() +
  # facet_wrap(
  #   vars(until_date),
  #   scales = "free_x"
  # ) +
  theme_pokemon +
  theme(aspect.ratio = 1, 
        legend.position = "none")
p_temporal_change_cumulative
ggsave("./output/p_emp_vs_sim_cumulative_ccdf_semilog.png", width = 86, height = 35, unit = "mm", dpi = 600)
```


## per gen

```{r}
# c_game_release_dates <- df_game_gen$game_release_date # inaccurate
c_series_gen_dates <- c(df_series_first_release_date_corrected$release_date, "2022-11-18")
```
シリーズのわけめで回す。第1世代、第2世代、でわけるのが一番いいわということで。`c_series_gen_dates`に関しては、第９世代より前ですよという日付が必要なので追加している。
こちらを主体にしている。上のもやってもいいけど…。
ここでそもそもn=0なやつらをカウントしてしまってもいい気がする(`possible_pokemons_dfをつかってcompleteなりjoinなりする）。0からはじまるCCDFは今回に限っては意味がある気がしないでもない。ガクっとなるかもしれませんが…。→しました
`complete()`はこうすればMissingな行もつくってくれる：
> When used with continuous variables, you may need to fill in values that do not appear in the data: to do so use expressions like year = 2010:2020 or year = full_seq(year,1).

```{r}
# run
for (i in 1:(length(c_series_gen_dates) - 1)) { 
  # i from 1 to 8
  
  sim_random_df_in_threshold_duration <-
    sim_random_df |> 
    filter(
      release_date < c_series_gen_dates[i + 1] & # 第２，３，４，。。。世代のリリース日よりも前で（含まない）
      release_date >= c_series_gen_dates[i] # 第１，２，３，。。。世代のリリース日よりもあと(含む)
        )
  emp_df_in_threshold_duration <- df_gene |> 
    filter(is_pokemon) |> 
    filter(
      # 第２，３，４，。。。世代のリリース日よりも前で
      release_date < c_series_gen_dates[i + 1] & 
      # 第１，２，３，。。。世代のリリース日よりもあと
      release_date >= c_series_gen_dates[i] 
    )
  # 第１，２，３，。。。世代のポケモンの総数
  pokemon_max_pokedex_id <- df_game_gen$pokemon_cumulative_n[i] 
  sim_df <-
    sim_random_df_in_threshold_duration |> 
    count(run, pokedex_id) |> 
    complete(
      pokedex_id = 1:pokemon_max_pokedex_id, 
      run = 1:1000,
      fill = list(n = 0)
      ) |> 
    mutate(
      condition ="simulations"
    )
  # unlike emp_poke_df below, these dfs cannot count absent pokemons
  # because... too many possible combinations...
  # e.g., Galarian delta flying pichu VMAX huh?
  emp_card_type_df <- 
    emp_df_in_threshold_duration |> 
    count(card_name, card_type) |> # not card_type2 
    mutate(
      condition = "card name + type"
    )
  emp_card_df <- 
    emp_df_in_threshold_duration |> 
    count(card_name) |> 
    mutate( condition = "card name")
  emp_poke_df <- 
    emp_df_in_threshold_duration |> 
    count(pokedex_id) |> 
    complete(
      pokedex_id = 1:pokemon_max_pokedex_id,
      fill = list(n = 0)
    ) |>
    mutate(condition = "pokemon name")
  
  # merge and calculate ECDF here
  
  merge_df <- bind_rows(sim_df, emp_card_type_df, emp_card_df, emp_poke_df) |> 
    group_by(condition, run) |> 
    arrange(desc(n)) |>
    mutate(ecdf = 1 - ecdf(n)(n - .01)) |> 
    left_join(df_names, by = "pokedex_id") |> 
    mutate(series_gen_era = i) 
  if (i == 1) {
    # initialize
    temporal_change_per_gen_df <- tibble(
      n = integer(),
      ecdf = double(),
      run = integer(),
      series_gen_era = integer()
    )
  } 
    temporal_change_per_gen_df <- bind_rows(
      temporal_change_per_gen_df,
      merge_df
    )
} # 1min?
temporal_change_per_gen_df |> dim() 
# 4381065 x 14 when all possible pokemons 
# (fill n = 0 for pokemons that had never appeared in that era), 
# 4048655 x 14 originally.
temporal_change_per_gen_df |> 
  filter(condition == "simulations", 
    run == 1,
    series_gen_era == 1
  ) 
```
### temporal change per gen, CCDF for each pokemon_gen

新世代がどれくらいバイアスかかってるかみてみたい

```{r}
# run
for (i in 1:(length(c_series_gen_dates) - 1)) { 
  # i from 1 to 8
  
  sim_random_df_in_threshold_duration <-
    sim_random_df |> 
    filter(
      release_date < c_series_gen_dates[i + 1] & 
      release_date >= c_series_gen_dates[i] 
    ) 
  
  emp_df_in_threshold_duration <- df_gene |> 
    filter(is_pokemon) |> 
    filter(
      release_date < c_series_gen_dates[i + 1] & 
      release_date >= c_series_gen_dates[i] 
    )
  
  # 第１，２，３，。。。世代のポケモンの総数
  pokemon_max_pokedex_id <- df_game_gen$pokemon_cumulative_n[i] 
  sim_df <-
  # df_view <- 
    sim_random_df_in_threshold_duration |> 
    count(run, pokedex_id) |> 
    complete(
      pokedex_id = 1:pokemon_max_pokedex_id, 
      run = 1:1000,
      fill = list(n = 0)
      ) |> 
    mutate(
      condition ="simulations"
    ) |> 
    left_join(
      df_names |> select(pokedex_id, pokemon_gen), 
      by = "pokedex_id"
    ) |> 
    group_by(pokemon_gen, run) |> 
    arrange(desc(n)) |>
    mutate(ecdf = 1 - ecdf(n)(n - .01)) 
  
  emp_poke_df <-
  # df_view <-
    emp_df_in_threshold_duration |> 
    count(pokedex_id) |> 
    complete(
      pokedex_id = 1:pokemon_max_pokedex_id,
      fill = list(n = 0)
    ) |> 
    mutate(condition = "pokemon name") |> 
    left_join(
      df_names |> select(pokedex_id, pokemon_gen), 
      by = "pokedex_id"
    ) |> 
    ungroup() |> 
    group_by(pokemon_gen) |> 
    arrange(desc(n)) |>
    mutate(
      # rank = rank(-n, ties.method = "max"),
      ecdf = 1 - ecdf(n)(n - .01)
    ) 
  # merge and calculate ECDF here
  
  merge_df <-
  # df_view <- 
    bind_rows(sim_df, emp_poke_df) |> 
    mutate(series_gen_era = i) 
  if (i == 1) {
    # initialize
    temporal_change_per_gen_per_pokemon_gen_df <- tibble(
      n = integer(),
      ecdf = double(),
      run = integer(),
      series_gen_era = integer(),
      pokemon_gen = integer()
    )
  } 
    temporal_change_per_gen_per_pokemon_gen_df <- bind_rows(
      temporal_change_per_gen_per_pokemon_gen_df,
      merge_df
    )
} # 1min?
```

```{r}
temporal_change_per_gen_per_pokemon_gen_df |> 
  filter(condition == "pokemon name", pokemon_gen == 1, series_gen_era == 2)
temporal_change_per_gen_per_pokemon_gen_df |> 
  filter(condition == "pokemon name", pokemon_gen == 2, series_gen_era == 2)
```



### temporal change, temporal

これも悪くはないと思うんだけど別に大したこと言えないのでとりあえず放置。・・・と思っていたけど、重複がどれくらい起きているかというのはわりと重要な気がする。重複を許していれば毎回ある程度は重複するけど、それがどれくらい偏っているか。リリースによってはかなり偏っているのでは？それと帰無モデルとの差をみられないか。

```{r}
# sim_random_df # 100 run
# df_gene # `card_name` = group by distinct names, `pokemon_name` = group by distinct pokemons.
release_dates <- df_series2 |> 
  filter(!is.na(pokemon_data_count)) |> 
  pull(release_date)
# number_of_releases <- length(release_dates) # 148
release_dates |> min() # 1999 Jan
release_dates |> max()# 2022 Sep
```


```{r}
# initialize
temporal_change_each_release_df <- tibble(
  n = integer(),
  ecdf = double(),
  run = integer(),
  release_date = Date()
)
# run length(months_seq)
for (i in 1:length(release_dates)) {
  sim_random_df_each_release <-
    sim_random_df |> 
    filter(release_date == release_dates[i])
  emp_df_each_release <- df_gene |> 
    filter(is_pokemon) |> 
    filter(release_date == release_dates[i])
  
  sim_df <- 
    sim_random_df_each_release |> 
    group_by(run, pokedex_id) |> 
    mutate(n = n()) |> 
    ungroup() |> 
    group_by(run) |> 
    arrange(desc(n)) |>
    mutate(ecdf = 1 - ecdf(n)(n - .01)) |> 
    select(run, n, ecdf) |> 
    distinct() |> 
    ungroup() %>%
    mutate(
      condition = rep("simulations", nrow(.))
    )
  emp_long_df <- 
    emp_df_each_release |> 
    group_by(card_name) |> 
    summarise(n = n()) |> 
    arrange(desc(n)) |> 
    mutate(ecdf = 1 - ecdf(n)(n - .01)) |>
    select(-card_name) |> 
    distinct() %>%
    mutate(condition = rep("card name", nrow(.)))
  emp_poke_df <- 
    emp_df_each_release |> 
    group_by(pokedex_id) |> 
    summarise(n = n()) |> 
    arrange(desc(n)) |> 
    mutate(ecdf = 1 - ecdf(n)(n - .01)) |>
    select(-pokedex_id) |> 
    distinct() %>%
    mutate(condition = rep("pokemon name", nrow(.)))
  merge_df <- bind_rows(sim_df, emp_long_df, emp_poke_df) %>%
    mutate(
      release_date = rep(release_dates[i], nrow(.))
    )
  cat(i, ", ")
  temporal_change_each_release_df <- bind_rows(temporal_change_each_release_df, merge_df)
} # 1min?
temporal_change_each_release_df |> dim() # 30915 x 5

```

plot

```{r}
temporal_change_each_release_df_with_release_gen <- temporal_change_each_release_df |> 
  left_join(
    df_series_releaes_gen_world_release_only |> select(release_date, release_gen),
    by = "release_date"
  ) |> 
  filter(condition != "card name")
```

すべてのリリースDateに対して、そのシミュレーション結果とポケモン名CCDFを重ねてみる。スナップショット。なんかだんだん最近になるにつれテールが長い気がする。
```{r}
temporal_change_each_release_df_with_release_gen |> 
  filter(condition == "simulations") |> 
  ggplot(
    aes(
      x = n,
      y = ecdf, 
      group = run
      )) +
  geom_line(alpha = .3, size = .3, colour = "grey20") +
  geom_count(alpha = .3, size = .3, colour = "grey20") +
  geom_point(
    data = temporal_change_each_release_df_with_release_gen |> 
      filter(condition == "pokemon name"),
    aes(colour = release_date),
    alpha = 1,
    size = .5
  ) +
  geom_line(
    data = temporal_change_each_release_df_with_release_gen |> 
      filter(condition == "pokemon name"),
    aes(colour = release_date),
    alpha = 1,
    size = .5
  ) +
  scale_x_continuous(
    # limits = c(1, 10000),
    trans = "log10", 
    breaks = 10^(0:10)
  ) +
  scale_y_continuous(
    # limits = c(0.0001, 1),
    trans = "log10", 
    breaks = 10^(0:-10)
    ) +
  scale_colour_viridis_c() + 
  facet_wrap(vars(release_date),
             nrow = 8
  ) +
  theme_pokemon +
  theme(
    aspect.ratio = 1
  )
ggsave("./output/temporal-change-each_release2.png", width = 5120, height = 2880, unit = "px")

```

2009-12-01とか何が起きてるんだろう。なんで1のやつが一人もいないんだろう
```{r}
df_gene |> filter(
  release_date == as.Date("2009-12-01")
)
```



### omake: frequency distribution of 1000 runs

```{r}
df_count_pokemon_random_draw <- sim_random_df |> 
  group_by(run, pokedex_id) |> 
  mutate(n = n()) |> 
  ungroup() |> 
  group_by(run) |> 
  arrange(desc(n)) |>
  mutate(ecdf = 1 - ecdf(n)(n - .01)) |> 
  select(run, n, ecdf) |> 
  distinct() |> 
  ungroup()
df_count_pokemon_random_draw_mean <- df_count_pokemon_random_draw |> 
  group_by(n) |> 
  summarise(mean_ecdf = mean(ecdf))
df_count_pokemon_random_draw |> 
  # filter(run == 3) 
  ggplot(
    aes(
      x = n, 
      y = ecdf, 
      colour = run
    )
  ) +
  geom_line(aes(group = run), alpha = .01, size = .1) +
  geom_point(aes(group = run), alpha = .01, size = .1) +
  geom_line(
    data = df_count_pokemon_random_draw_mean,
    aes(y = mean_ecdf),
    colour = "red"
  ) +
  scale_x_continuous(
    limits = c(1, 130),
    trans = "log10", 
    breaks = 10^(0:10)
  ) +
  scale_y_continuous(
    limits = c(0.001, 1),
    trans = "log10", 
    breaks = 10^(0:-10)
    ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(size = .3),
    axis.ticks = element_line(size = .1),
    aspect.ratio = 1
  )
ggsave("./output/random-draw-1000runs.png", width = 3000, height = 2000, unit = "px")

```






### 多様性はどうか

1年くらいのウィンドウをずらしていって、なんとかの多様性を測る。個体数とVariant数が一致しないので難しいのだが…。個体数がでればなぁ　試合なら出るか。


## remix?

CharizardのほうがPokemonとしては多いが、Eeveeという名前のカードとCharizardという名前のカードだとEeveeのほうが多い、というのを比べてみる。

```{r fig.height = 15, fig.width = 4}
df_empirical_count_remix <- 
  df_gene |> 
  filter(!is.na(id)) |> # filter only pokemons. is_
  group_by(id) |> 
  mutate(count_pokemon = n()) |> 
  ungroup() |> 
  group_by(card_name) |> 
  mutate(count_name = n()) |> 
  ungroup() |> 
  mutate(name_type = paste(card_name, card_type)) |> 
  group_by(name_type) |> 
  mutate(count_name_type = n()) |> 
  ungroup() |> 
  select(starts_with("count_"), id, name_type, card_name, pokemon_name) |> 
  distinct() |> 
  mutate(pokemon_name = fct_reorder(pokemon_name, desc(count_pokemon))) |> 
  
```

### original vs remix


```{r}
df_empirical_count_remix2 <-
  df_empirical_count_remix |> 
  select(-count_name_type) |> 
  distinct() |> 
  filter(card_name == pokemon_name) |> 
  mutate(count_remix = count_pokemon - count_name) 
  # count_remix: ピカチュウのなかでまさにピカチュウという名前のカードを除いたカードの数
```

```{r}
# count_pokemon > count_name > count_name_type
df_empirical_count_remix2 |> 
  ggplot(aes(x = count_pokemon, y = count_remix)) +
  geom_count() +
  geom_abline(intercept = 1) +
  xlim(c(0, 140)) +
  ylim(c(0, 140)) +
  theme_pokemon +
  theme(aspect.ratio = 1)
```
格差を拡大させているのはひとつにはリミックスの多さであろうから（ぴかちゅうとかは人気が高い→リミックスたくさん、3枚とかしかまだ出てないようなポケモンはリミックスほとんどない場合が多い）、それを除いたらちょっと格差が是正されるわな→そんなにでした