イグノーベル賞Talent vs Luck論文にPythonとRで反論を試みる

10月3日のノーベル賞発表直前、駆け込み記事です。今回は2022年イグ・ノーベル経済学賞受賞論文を扱ってみます。

※使用環境
Python:3.10.6 (numpy:1.23.2)
R:4.2.0 (tidyverse:1.3.2, sqldf:0.4.11)
## 元論文の概要

さて、この論文、「なぜ成功は、最も才能のある人ではなく、最も幸運な人にもたらされることが多いのかを数学的に説明した」(Improbable Research, 2022)ことが受賞理由と書かれています。この世は才能なのか運なのか、いわゆる「Talent vs Luck」論争に対する数値シミュレーションを行い、やっぱり運の方が重要だと言っているわけです。

プルチーノ氏は、以下のような条件でシミュレーションしてます(Pluchino et al., 2018)。
・1000人の人間が500のイベントに遭遇する世界のシミュレーションを実施
・人間の才能は0～1間の値が正規分布で与えられ、一生変わらない
・成功を示すための資産は平等に10からスタート
・20歳から60歳までの年齢で半年に一回、つまり80回のイテレーションを行う
・遭遇するイベントは不幸か幸運かに50%の確率で分かれる
・幸運イベントは才能に応じて勝ち取ることができ、資産が2倍
・不幸イベントは才能に関わらず無条件で資産が半減

シミュレーションの結果は以下の通り
・成功は冪乗則とりわけパレートの法則を再現した(上位20%が富の80%を持っている)
・成功したのは幸運イベントの回数が多かった人であり才能の高かった人ではない

論文の後半は、幸運の格差の是正のためには才能に応じた補助金よりも一律の補助金の方が効率が良いことを示唆しています。

## 反論の要旨
私は素人ですが、モデルに恣意的な単純化が行われている気がしてなりません。すぐに思いつく疑問は二つです。まず、イベントに遭遇できる契機は、その人間の外交性、つまり才能に依存しているのではないかということ。社交の才能がある人間の方がより多くイベントに遭遇するようなモデルで計算してみたらどうなるでしょうか。次に、不幸の時に資産の半減を回避することもまた、才能があれば可能なのでは、ということです。神経症傾向が強く事故を避けやすい、事故に遭っても適切な保険に入っているかどうかなどですね。これも、そのようなモデルで計算してみたくなります。（なお、外交性、神経症傾向という言葉は心理学のビッグファイブを意識しています。これらの性格特性は、ある程度まで遺伝で決まると言われています。そもそもの話として、性格の方向性だけでも5つあるのに、才能が1つの軸上に並べられることなどないわけですよ。）

いやー、どう考えても、世の中運だけなんてことはないと思うんですけどね。受賞論文以外は調べてないのですが、それが2018年だし、既にいろいろ議論されてるものなのかもですが。

とにかく、以下の手順で進めていきましょう。
・まずPythonで上記モデルに近いものを実装してみる
・出力したデータをRでグラフにして、論文の結果が再現できるか確認
・上記疑問点を加味した別モデルに書き換えて、比べてみる

なお、最初、Javaで書き始めたんですが、確率を生成するためにApache Mavenを入れたり大変そうだったので、Pythonにしました。Pythonの中でPandas等グラフが書けるツールも試してみようかと思ったのですが、Rの方が素早くできそうなので、今回はやめておきました。(急がないとノーベル賞が発表されてしまって、イグノーベル賞について書く期を逃します！)

## Python実装

まず、人間クラスを作ります。論文と同じ値を変数に持たせて、そのインスタンスに対するイベント判定も同クラス内にメソッドで書いてしまいます。

```python sumulator.py
import numpy as np

class Human():
    def __init__(self):
        self.talent = np.random.normal(0.6, 0.1, None)
        self.wealth = 10
        self.lucky_cnt = 0
        self.unlucky_cnt = 0
        
    def face_event(self):
        freq = np.random.poisson(lam=0.45)
        for i in range(1,freq):
            lucky = bool(np.random.randint(0,2))
            if lucky:
                if(np.random.random() < self.talent):
                    self.wealth *= 2
                self.lucky_cnt +=1
            else:
                self.wealth =  self.wealth/2
                self.unlucky_cnt +=1

    def header(): return "wealth,talent,lucky_cnt,unlucky_cnt"
    def to_str(self): return f"{self.wealth},{self.talent},{self.lucky_cnt},{self.unlucky_cnt}" 

```

才能の正規分布には、numpy.random.normal()を用います。meanが0.6、標準偏差が0.1というのは論文と同じ値になります。コンストラクタで、才能と、初期値10の資産、あとは後でグラフに書くイベント遭遇回数を持たせます。

幸か不幸かのイベントの分け方は、random.randint(0,2)をbooleanにキャストしました。これで50%です。幸運なとき、才能値を確率として資産が2倍になるという条件は、0～1のランダム値が才能値よりも小さいかどうかで判定できます。困ったのは、イベントへの遭遇です。

元論文によると、正方形の2次元空間に1000人の人間が配置され、その中を500個のイベントが動き回っているというシミュレーションのようです。動くんすか！？ さすがに難易度高くないですか・・・。でも、特定の期待値で遭遇するイベントの確率を求めるんだから、今回はポアソン分布でやってみましょう(厳密には異なるのでしょうが・・・)。 ポアソン分布のラムダ値は、イベントに遭遇する期待値です。1000人の人間の中を500個のイベントが動き回っているのだから、期待値は0.5・・・。ということで、λ=0.5のポアソン分布に従って遭遇回数を決め、その分だけ幸・不幸判定をループする、と。

※これは試した後の調整ですが、0.5だと、パレートの法則(上位20%が80%の富を持つ)がなかなか再現できません。80%になることもありますが、平均としては85%近辺に落ち着いてしまいます。やはり、500個のイベントは全部ぶつからないということでしょうか。450個だけぶつかる・・・λ=0.45で、同値を80%くらいにできたので、そちらに書き換えて採用します。

同じファイルの続きです。

```python sumulator.py
def run_simulation(class_arg):
    population = []
    for n in range(1000):        
        population.append(class_arg())

    for i in range(80):
        for man in population:
            man.face_event()
    try:
        with open("data.csv", "w") as file:
            file.write(Human.header()+"\n")
            for man in population:
                file.write(man.to_str()+"\n")
    except:
        print("write file error.")

run_simulation(Human)
```
人間インスタンスを1000個作って、80回イベントを回します。その結果を、CSVファイルに書き込みます。今後、Humanクラスを継承して別のタイプの人間を作るつもりなので、クラスを引数にして走らせることにしました。pythonの関数ってクラス名を引数にして、その引数に「()」付けたらインスタンスをnewできるんすね・・・。とにかく急ぎましょう。次はRです。
## R実装
```R script.R
library(tidyverse)
library(sqldf)

dat <- read.csv("data.csv")
# Check talent distribution
ggplot(data = dat, aes(x = talent)) +
geom_histogram(binwidth = 0.01, fill = "white", colour = "black")
```

まずロードして、才能が元論文と近い正規分布になっていることをチェックします。
<img width="320" alt="plot1.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/26b73bb5-da82-e834-4567-893fab1d1530.png">


```R script.R
# Check wealth distribution
# log-lin
ggplot(data = dat, aes(x = wealth)) +
geom_histogram(binwidth = 50, fill = "white", colour = "black") +
scale_y_log10()
```

そして、資産の分布。（左のy軸は対数表記に変えています。）
<img width="320" alt="plot2.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/998d5deb-c47f-2554-a6bd-03c47f1ef673.png">

```R script.R
# Check Pareto's "80-20" rule
sum(sqldf("SELECT wealth FROM dat 
    ORDER BY wealth DESC
    LIMIT (SELECT 0.20 * COUNT(*) FROM dat)")) /
(sqldf("SELECT SUM(wealth) FROM dat"))

sum(sqldf("SELECT wealth FROM dat
       ORDER BY wealth
       LIMIT (SELECT 0.80 * COUNT(*) FROM dat)")) /
(sqldf("SELECT SUM(wealth) FROM dat"))
```
上位20%の富が何割か、計算します。(ややこしかったのでSQL使ってしまいましたごめんなさい。)
```
0.7992045
0.2007955
```
79%。これが近ければ、モデルとして合格にします！ 次が論文の大事なポイント、富と才能の相関です。

```R script.R
# Check wealth distribution with talent
ggplot(data = dat, aes(x = wealth, y = talent)) +
geom_point() +
scale_x_log10()

ggplot(data = dat, aes(x = talent, y = wealth)) +
geom_point()
```
ここで、富と才能が必ずしも対応していないこと、最も富を持っているのは、最も才能を持ったグループではないこと、というのが論文著者の言いたいことのようです。再現できました。
<img width="320" alt="plot3.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/0fe621f6-c1cf-ca39-7c36-d270e32cf3d7.png"><img width="320" alt="plot4.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/4deb917d-b7cf-d32f-5a81-b1a318bb8339.png">
ただ、才能があるほうが富を持ちやすい傾向はあり、何回かやると、才能のある人物だけ富が突出しているケースにも遭遇します。(ゆえに論文はチェリーピッキングだ、と言っているわけではないです・・・。)　最後は蛇足ですが、運の寄与を示すグラフも再現しておきます。
```R script.R
# Check wealth distribution with luck
ggplot(data = dat, aes(x = wealth, y = lucky_cnt)) +
geom_point() +
scale_x_log10()

# Check wealth distribution with unluck
ggplot(data = dat, aes(x = wealth, y = unlucky_cnt)) +
geom_point() +
scale_x_log10()
```
富と幸運イベントの正の相関、富と不幸イベントの負の相関が確認できます。
<img width="320" alt="plot5.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/027fe9e4-7ace-7d0f-78cb-1e7137605ea6.png"><img width="320" alt="plot6.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/74b650a5-040a-b448-2973-1dfc3270a4ea.png">
以上、論文の主張と近いグラフが描かれるモデルを作ることができました！

## 別モデル1

では、社交性を才能としたケースを追加しましょう。Humanクラスを継承してSociableHumanクラスを作ります。

```python sumulator.py
class SociableHuman(Human):
    def __init__(self):
        super().__init__()

    def face_event(self): # overrides
        freq = np.random.poisson(lam=self.talent-0.15) 
        for i in range(1,freq):
            lucky = bool(np.random.randint(0,2))
            if lucky:
                if(np.random.random() < self.talent):
                    self.wealth *= 2
                self.lucky_cnt +=1
            else:
                self.wealth =  self.wealth/2
                self.unlucky_cnt +=1

run_simulation(SociableHuman)
```
イベント遭遇のポアソン分布の期待値を、才能に従って分布させます。才能の中央値が0.6、元のラムダが0.45なので、才能から0.15を引いたところを中心に、期待値をバラつかせます。つまり、良く動く才能のある人は、それだけイベントに遭遇しやすい。で、このモデルを実行すると、富と才能の分布は以下の形になります。
![plot8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/b155c783-2801-0d39-478e-69b44a71ca4d.png)
才能が高いほど、富が減りやすく増えやすいと。でも、これだけだと幸運の方が主要因であることに変わりはなく、反証になっていませんね。さらに悪いことに、才能があり幸運なケースが滅茶苦茶稼ぐので、パレートの法則が崩れてしまいます。なので、モデルとしては使えないことにしておきます。元論文では、富の分布を現実に近くできたことが説得材料になっているからです。
```
0.9114258
0.08857422
```

## 別モデル2

今度は、不幸に備える才能を加味した場合です。才能があれば、資産が半額になる不幸イベントを防げるコードにしました。

```python sumulator.py
class PreventiveHuman(Human):
    def __init__(self):
        super().__init__()

    def face_event(self): # overrides
        freq = np.random.poisson(lam=0.5) 
        for i in range(1,freq):
            lucky = bool(np.random.randint(0,2))
            if lucky:
                if(np.random.random() < self.talent):
                    self.wealth *= 2
                self.lucky_cnt +=1
            else:
                if(np.random.random() > self.talent):
                    self.wealth =  self.wealth/2
                self.unlucky_cnt +=1

run_simulation(PreventiveHuman)
```

このモデルはパレートの法則をほぼ満たします！
```
0.8172335
0.1827665
```
そして、才能と資産の分布をみてください。
<img width="320" alt="plot10.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/397a658f-6735-81ea-5d5b-dc672fbf0416.png"><img width="320" alt="plot11.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2866483/5f74c4a8-553a-f2f5-f9cb-c2e132a05d9e.png">

この場合には、才能と成功には相関関係があるし、最も稼いでいるのは、かなり才能のある方の人物です。個人的には、こちらの方が現実のモデルに近いのではないかなと思ったりします。

## 結論

とにかく、今日はここまでにします。論文のように厳密な方法ではないですが、シミュレーションを再現したうえで別の意見を述べるという目的は果たせたかと思います。結局のところ、才能が関係あるかどうかは、適用するモデルにによるということですね。ここまで読んでくださった方、ありがとうございましたー。

追記：GitHubにソースを上げておきました。
https://github.com/KentAnak/TvLSimulator

<b>Reference:</b>
Improbable Research. (2022). <i>Ignobel Prize Winners|2022.</i> https://improbable.com/ig/winners/#ig2022
Pluchino, A., Biondo, A. E., & Rapisarda, A.(2018). Talent vs. Luck: The Role of Randomness in Success and Failure,<i>Advances in Complex Systems</i>, 21[3,4]. https://doi.org/10.1142/S0219525918500145