幸运世界模型的蒙卡模拟（附python代码）
笔者对2022年搞笑诺贝尔经济学奖论文TALENT VERSUS LUCK: THE ROLE OF RANDOMNESS IN SUCCESS AND FAILURE有点兴趣，遂在下班前花半小时写了个简单的python代码。
笔者先做了一个简单的一维模拟，对论文的情况简化很多：有1000人，智商正态分布IQ=Gauss（100,20）, 启动资金10000，每年发年薪1000。幸运降临的时候收入1000*智商/100, 不幸降临的时候收入减去1000元。假设每年幸运和不幸都会随机降临在700人头上。追踪60年，重复实验1000次。
结果显示和赢家的智商高、输家的智商低，与智商还是相关的。对这个简单的小世界模型，因为正态分布的关系，实际上通过简单的算术就可以估计出结果。要赢的话智商要高，但是智商过高的人数比较少，所以就稳定在130左右。同理，输家的正态分布均值在60左右。比较有意思的是60年后的资产分布，是不对称的（可以计算skew量化），看着像朗道分布，不知道修改模型后会不会出现原文提到的Pareto分布（二八法则）？
等有空详细看看那篇文章，对模型进行修改，估计会有不少有趣的结果。原文作者是意大利天体物理的，这种通过简单的模型来设定世界的思路还比较有意思，特别是社科、经济、复杂现象等科研人员应该可以做出很多文章来。
输家智商分布
赢家智商分布
输家资产
赢家资产
import os, sys
from stat import *
import numpy as np
from math import sqrt
from scipy.stats import norm
from pylab import plot, show, grid, axis, xlabel, ylabel, title
import matplotlib.pyplot as plt
def trial():
## label 1000 persons, they change from each year
N = 1000
asset0 = 10000
salary = 1000
luck = 1000
unluck = -1000
Nt = 60
label = [i for i in range(N)]
mu, sigma = 100, 20 # mean and standard deviation
s = np.random.normal(mu, sigma, N)
ability = [i/100 for i in s]
start_asset = [asset0 for i in range(N)]
individual_track = []
track_year = []
## print start_asset
for yy in range(Nt):
temp = [0 for i in range(N)]
## print "year", yy
for i in range(int(N*0.7)):
randLuck = np.random.uniform(0,N)
randUnlucky = np.random.uniform(0,N)
whoisLucky = int(randLuck)
whoisUnlucky = int(randUnlucky)
start_asset[whoisLucky] = start_asset[whoisLucky] + ability[whoisLucky]*luck
start_asset[whoisUnlucky] = start_asset[whoisUnlucky] + -1*np.random.uniform(1000,5000)
#print "year=", year
#print whoisLucky, whoisUnlucky, ability[whoisLucky]
for i in range(N):
start_asset[i] += salary
## print "check", i, start_asset[i]
for i in range(N):
temp[i] = start_asset[i]
## print "end of year", start_asset
track_year.append(temp)
#print track_year
final = track_year[Nt-1]
sum1 = 0
#print "who is most unlucky", final.index(min(final)), "with money: ", min(final)
#print "who is most luck", final.index(max(final)), "with money: ", max(final)
badMan = final.index(min(final))
goodMan = final.index(max(final))
badManIQ = ability[badMan]
goodManIQ = ability[goodMan]
badManAsset = min(final)
goodManAsset = max(final)
return badManIQ, badManAsset, goodManIQ, goodManAsset
Nexp = 1000
result_badManIQ = []
result_badManAsset = []
result_goodManIQ = []
result_goodManAsset = []
for i in range(Nexp):
badManIQ, badManAsset, goodManIQ, goodManAsset = trial()
result_badManIQ.append(badManIQ)
result_badManAsset.append(badManAsset)
result_goodManIQ.append(goodManIQ)
result_goodManAsset.append(goodManAsset)
plt.hist(result_badManIQ)
plt.show()
plt.hist(result_goodManIQ)
plt.show()
plt.hist(result_badManAsset)
plt.show()
plt.hist(result_goodManAsset)
plt.show()
