大沢統計第1章(2)

前回に続いて大沢統計をプログラミングしてみました。
内容は前回の記事とほぼ一緒です。最初にnum_person人にmoney=30をランダムで配ります。その後にランダムで選ばれた人が、ランダムに選ばれた人にmoney=1を渡すという過程を10000回繰り返します。結果として、どんな人も一度は大きなmoneyを持つことが出来ます。しかし、その確率はexp(-所有するmoney)に比例し、指数分布に従う状態で現れるので、大きなmoneyを所得している時間は、悲しいことかな、短いのです。この結果は、長い時間をかければ、ある状態が等確率で起こることを示唆しています。(簡単にいうと、長い時間をかければ、お金持ちになっている時間やお金を持っていない時間はみんな大体等しくなる。)

#! /usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import copy

### initial set
np.random.seed(0)
money=30
num_person=5
ave=money/num_person
num_exchange=10000
persons=np.zeros(num_person)
persons_trajectory=[]

### initial fig set
fig, axPlot = plt.subplots(figsize=(5.5,5.5))
divider = make_axes_locatable(axPlot)
axHisty = divider.append_axes("right", 1.2, pad=0.1, sharey=axPlot)
plt.setp( axHisty.get_yticklabels(),visible=False)

for i in range(money):
    persons[np.random.randint(0,num_person)]+=1
persons_trajectory.append(persons)

### exhange money
for i in range(num_exchange):
    out=np.random.randint(0,num_person)
    if persons[out]==0:
        tmp=copy.copy(persons)
        persons_trajectory.append(tmp)
        continue
    persons[out]-=1
    persons[np.random.randint(0,num_person)]+=1
    tmp=copy.copy(persons)
    persons_trajectory.append(tmp)

### for fig
persons_trajectory=np.array(persons_trajectory).T

for i in range(num_person):
    axPlot.plot(persons_trajectory[i])
    axHisty.hist(persons_trajectory[i], bins=25, orientation='horizontal',alpha=0.6)

plt.savefig("ch1-2.png")
plt.draw()
plt.show()

下記の図は、左の図の縦軸がお金、横軸はステップを示しています。それぞれの色は5人に対応しています。ヒストグラムで示した分布が大体似ていることから、５人全員が大体同じ確率で高所得、低所得を経験することがわかります。