センス8の感応者が本当にランダムだったら

少し古いのですが、センス8というNetflixのドラマを見始めました。とても面白いです。世界中でランダムに選ばれた8人が超能力に目覚め、力を合わせて困難に立ち向かう、という内容です。この8人の「感応者」は世界中に散らばっておりナイジェリア人や韓国人を含んでいます。おそらく多様性についての話が展開するのでしょうか、LGBTの人も複数含んでいます。

さて、この8人が本当に「ランダム」なのか気になってきました。どうも白人が多いような気がするのです。世界中からランダムに8人を選ぶとこの組み合わせになる可能性はどのくらいでしょうか。

調べて驚いたのですが、世界の人種構成についてのちゃんとしたデータは存在しません。そもそも「人種」や「民族」というのに明確な定義がないため統計に資さないのでしょう。ですので、人種の代理として、大陸ごとの構成について検討してみました。

すなわち、世界中からランダムに8人を選んだ時に、各大陸の構成がドラマのようになる確率はどの程度でしょうか。[wikipedia](https://en.wikipedia.org/wiki/List_of_continents_by_population) の世界人口のデータと比べてみましょう。（見やすくするために、オセアニアと南極のデータは省いています。）


```python
%config InlineBackend.figure_formats = {'png', 'retina'}

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 12, 6

world_pop = pd.Series({
    'Asia': 4436224000,
    'Africa': 1216130000,
    'Europe': 738849000,
    'North America': 579024000,
    'South America': 422535000,
    'Oceania': 38304000,
    'Antarctica': 1106
                      })

world_pop = world_pop[~world_pop.index.isin(['Oceania', 'Antarctica'])]

sense8_cluster = pd.Series({
    'Capheus': 'Africa',
    'Sun': 'Asia',
    'Nomi': 'North America',
    'Kala': 'Asia',
    'Riley': 'Europe',
    'Wolfgang': 'Europe',
    'Lito': 'North America',
    'Will': 'North America'
                  })

sense8_tabs = sense8_cluster.value_counts()

pop_summary = pd.DataFrame({'real world': world_pop, 'sense8': sense8_tabs}).fillna(0)
pop_summary.plot(kind = 'pie', subplots=True, legend=False)
plt.tight_layout()
```


<img width="855" alt="output_1_0.png" src="https://qiita-image-store.s3.amazonaws.com/0/145500/a1d2c6f2-2e78-2f5f-9580-350878709f4d.png">



世界人口の大陸ごとの割合が左、センス8の8人の大陸ごとの割合が右です。アジア大陸が少なく、北米大陸が多いことが分かります。
世界中から8人を選んだ時に、ドラマの組み合わせになる確率を計算します。これは高校数学ですね。


```python
n_total = sense8_tabs.sum()
pop_summary['real_fraction'] = pop_summary['real world'] / pop_summary['real world'].sum()

denom = 1
for i in sense8_tabs:
    denom *= np.math.factorial(i)

sense8_p = np.prod(pop_summary['real_fraction'] ** pop_summary['sense8']) * (np.math.factorial(n_total) / denom) * 100

print("sense8 cluster would look like this in only {:0.3f} % of the time".format(sense8_p))
```

    sense8 cluster would look like this in only 0.048 % of the time
    

というわけで、なんと0.048%です。非常に低く感じます。ただし、5大陸から8人を選んだときに可能な組み合わせが非常に多いため、0.048%という数字が大きいのか小さいのかよく分かりません。ですので、切り口を変えてみましょう。

ブートストラップ法はどうでしょうか。世界中から8人を選ぶ、というのは典型的な多項分布といえます。```np.random.multinominal()```で容易にシミュレーションできます。ランダムに8人を選ぶゲームを10000回繰り返した時、各大陸が占める割合の分布はどのようになるのでしょうか。その分布から観測値が大きく逸脱していればランダムに選んでいない、と推定できます。感応者がランダムに選ばれたと仮定して、ドラマの割合が出る確率を計算します。


```python
n_sample = 10000
rnd = np.random.RandomState(67)

sim_result = np.zeros([n_sample, pop_summary.shape[0]])

for i in range(n_sample):
    sim_result[i, :] = rnd.multinomial(n_total , pop_summary['real_fraction']) / n_total

sim_mean = np.median(sim_result, axis = 0)
ci_lower = np.quantile(sim_result, axis = 0, q = 0.025)
ci_upper = np.quantile(sim_result, axis = 0, q = 0.975)

p_value_raw = np.sum(sim_result <= sense8_fraction.values, axis = 0) / n_sample
p_value = np.minimum(p_value_raw, (1-p_value_raw)) * 2

sim_summary = pd.DataFrame({'observed': sense8_fraction, 
                           'sim_mean': sim_mean,
                           'ci_lower': ci_lower,
                           'ci_upper': ci_upper, 
                           'p_value': p_value})
sim_summary['ci_lower_centered'] = sim_summary['sim_mean'] - sim_summary['ci_lower']
sim_summary['ci_upper_centered'] = sim_summary['ci_upper'] - sim_summary['sim_mean']

sim_summary = sim_summary.sort_values('sim_mean')
new_index = [i + "**" if sim_summary.loc[i, 'p_value'] < 0.01 else i + '*' if sim_summary.loc[i, 'p_value'] < 0.05 else i for i in sim_summary.index]
sim_summary.index = new_index

sim_summary
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>observed</th>
      <th>sim_mean</th>
      <th>ci_lower</th>
      <th>ci_upper</th>
      <th>p_value</th>
      <th>ci_lower_centered</th>
      <th>ci_upper_centered</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>North America**</th>
      <td>0.375</td>
      <td>0.000</td>
      <td>0.00</td>
      <td>0.250</td>
      <td>0.0038</td>
      <td>0.000</td>
      <td>0.250</td>
    </tr>
    <tr>
      <th>South America</th>
      <td>0.000</td>
      <td>0.000</td>
      <td>0.00</td>
      <td>0.250</td>
      <td>0.7620</td>
      <td>0.000</td>
      <td>0.250</td>
    </tr>
    <tr>
      <th>Africa</th>
      <td>0.125</td>
      <td>0.125</td>
      <td>0.00</td>
      <td>0.500</td>
      <td>0.7726</td>
      <td>0.125</td>
      <td>0.375</td>
    </tr>
    <tr>
      <th>Europe</th>
      <td>0.250</td>
      <td>0.125</td>
      <td>0.00</td>
      <td>0.375</td>
      <td>0.0780</td>
      <td>0.125</td>
      <td>0.250</td>
    </tr>
    <tr>
      <th>Asia</th>
      <td>0.250</td>
      <td>0.625</td>
      <td>0.25</td>
      <td>0.875</td>
      <td>0.1024</td>
      <td>0.375</td>
      <td>0.250</td>
    </tr>
  </tbody>
</table>
</div>



北アメリカ大陸のP valueが0.0038であり、感応者がランダムに選ばれているという帰無仮説を棄却できます。すなわち、北米大陸の感応者が不自然に多い、といえます。これを可視化してみます。


```python
plt.scatter(sim_summary['sim_mean'], sim_summary.index, color = 'blue', s = 200)
plt.errorbar(sim_summary.sim_mean, sim_summary.index, xerr = sim_summary[['ci_lower_centered', 'ci_upper_centered']].values.T, fmt = 'none')
plt.scatter(sim_summary['observed'], sim_summary.index, color = 'red', s = 100)

plt.xlabel("fraction")

plt.show()
```


<img width="776" alt="output_7_0.png" src="https://qiita-image-store.s3.amazonaws.com/0/145500/dd5f95fb-230b-7dc2-cfd6-1ba2c9fae8f7.png">



赤が実際のドラマの割合、青がシミュレーションでの中央値と95%区間です。北米大陸が突出して多いことが分かります。アジア大陸が95%区間の下限あたりにあります。アジア大陸を削って北米大陸を増やしている、とも言えます。北米大陸のキャラクター3人とも白人ですから、アジア人を削って白人を増やした、white-washingだと怒られそうな気もします。では、PC (politically correct)な組み合わせはどのようなものでしょうか。先ほどのシミュレーションで頻度の高い組み合わせを抽出します。


```python
row_count = np.unique(sim_result, axis = 0, return_counts=True)

sim_df = pd.DataFrame(row_count[0], columns=pop_summary.index)
sim_df['freq'] = row_count[1] / n_sample
sim_df = sim_df.sort_values('freq', ascending=False)

sim_df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Africa</th>
      <th>Asia</th>
      <th>Europe</th>
      <th>North America</th>
      <th>South America</th>
      <th>freq</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>124</th>
      <td>0.125</td>
      <td>0.750</td>
      <td>0.125</td>
      <td>0.000</td>
      <td>0.0</td>
      <td>0.0422</td>
    </tr>
    <tr>
      <th>125</th>
      <td>0.125</td>
      <td>0.875</td>
      <td>0.000</td>
      <td>0.000</td>
      <td>0.0</td>
      <td>0.0358</td>
    </tr>
    <tr>
      <th>176</th>
      <td>0.250</td>
      <td>0.750</td>
      <td>0.000</td>
      <td>0.000</td>
      <td>0.0</td>
      <td>0.0347</td>
    </tr>
    <tr>
      <th>175</th>
      <td>0.250</td>
      <td>0.625</td>
      <td>0.125</td>
      <td>0.000</td>
      <td>0.0</td>
      <td>0.0336</td>
    </tr>
    <tr>
      <th>123</th>
      <td>0.125</td>
      <td>0.750</td>
      <td>0.000</td>
      <td>0.125</td>
      <td>0.0</td>
      <td>0.0310</td>
    </tr>
  </tbody>
</table>
</div>



上位5位がこの組み合わせです。すなわちsense8 - PC edition!!を作ると大陸ごとの構成はこうなります。

|大陸|人数|
|:------------|:--------|
|アフリカ大陸|1～2人|
|アジア大陸|5～7人|
|ヨーロッパ人|0～1人|
|北米大陸|0～1人|
|南米大陸|0人|

ランダムに選ぶと、アジアとアフリカの話になりそうです。かえって多様性が減ってしまい、僕らが感じるPCとずれてきますね。というわけで、私はPC edition!!ではなく、ありのままのセンス8を楽しみたいと思うのです。