今日書いたコレを元に、シミュレーションっぽいものを書いて遊んでいました。
社会シミュレーション的なことをやるとき、オブジェクト指向の考え方で、人や粒子をオブジェクトとしてプログラムを組む人が多いようです。経営なんちゃら学科で社会シミュレーションやってるって人はそのためにJava使ってる言ってたのを聞いたことあります。ただ、自分はもともとC言語でシミュレーションをやっていたクセで、行列に値をもたせたほうが慣れてるのでそうします。
ブラジルナッツ効果を質量（というか移動しやすい粒子としにくい粒子の違い）によって説明できたらスマートなんじゃないかなという感じのシミュレーションです。
二次元平面上のランダムな場所を選び、その上にある粒子が上下左右ランダムな位置に移動しようとする→移動先が空白なら移動が完了、移動先に粒子がいた場合移動できず、その粒子が自分より軽ければ弾き飛ばす、というふうに設定しています。
import numpy as np
from pylab import *
def initial(N,A,B): #初期配置
grid = np.zeros(N-(A+B)) #空の領域
grid = np.r_[grid,(np.ones(A))]
grid = np.r_[grid,(np.ones(B))*2]
np.random.shuffle(grid)
grid = grid.reshape(L,L)
return grid
def move(x,y,coin): #移動
if coin < 1:
x_trans = (x+1) % L
y_trans = y
elif coin < 2:
x_trans = x
y_trans = (y+1) % L
elif coin < 3:
x_trans = (x-1)
y_trans = y
else:
x_trans = x
y_trans = (y-1)
if space[x_trans][y_trans] == 0:
space[x][y], space[x_trans][y_trans] = 0, space[x][y]
elif space[x][y] > space[x_trans][y_trans]: #移動先の人の質量が小さい時
move(x_trans, y_trans, coin) #移動先の人を同じ方向にはじき出す
else: #移動先の人の質量が大きいまたは同等のとき
pass #はじき出せない
L = 50
N = L*L
A = N/3 #質量1の人数
B = N/3 #質量2の人数
space = initial(N,A,B)
for x in range(N*N): #N×N回の試行
x = np.random.random()*L
y = np.random.random()*L
coin = np.random.random()*4
move(x,y,coin)また、同じ質量どうしで固まっている具合を表すパラメータを解析する関数も用意しておきます。
def search_all_directions(x,y):
searcher = space[x][y]
count = 0
xi = (x+1) % L
yi = (y+1) % L
xd = (x-1)
yd = (y-1)
if space[xi][y] == searcher:
count += 1
elif space[x][yi] == searcher:
count += 1
elif space[xd][y] == searcher:
count += 1
elif space[x][yd] == searcher:
count += 1
return count
def homophily(): #このパラメータが大きいほど、同じ質量の人が固まっているはず
homo = np.zeros(3)
for x in range(L):
for y in range(L):
if space[x][y] == 0:
homo[0] += search_all_directions(x,y)
elif space[x][y] == 1:
homo[1] += search_all_directions(x,y)
elif space[x][y] == 2:
homo[2] += search_all_directions(x,y)
else:
print "error"
return homo初期状態をプロットした結果。
hot()
pcolor(space)
colorbar()
show()
シミュレーション後。心なしか、同じ種類で固まっているようにも見えます。特に黒色（空白のスペース）。
先ほど定義した homophily() 関数を見ると、シミュレーション前に同じ種類のセルが固まっていることを示すパラメータが、0（空白のセル）、1（軽い粒子）、2（重い粒子）がそれぞれ、
[ 0.81534772 0.81032413 0.79231693]
シミュレーション後が
[ 0.85131894 0.83913565 0.83913565]
となっており、どの種類もシミュレーション前より固まっており、ブラジルナッツ効果っぽいものがあらわれてるんじゃないかなーって感じの結果が出ている気がします。ただし、ちゃんと説明すると言えるようになるには、もっと厳密に検証しないといけないでしょう。計算時間やシミュレーション回数が足りているかどうかも怪しいです（そもそもどっかにバグがあってもおかしくない）。
そもそもこの前の『データの見えざる手』に書いていたシミュレーションのように、粒子をランダムに動かす試行をするだけである程度は偏った分布（指数分布）になるはずなので。
【Python】『データの見えざる手』第一章にあったU分布（指数分布）のシミュレーション - 歩いたら休め
思いつきでやったのですが、もしちゃんと研究するなら、比較対象として質量の大小による相互作用の違いを消してシミュレーションをしたり、粒子の割合を変えてみたり、いろいろやらなきゃいけないのでダルいでしょう。誰かやっといてください。
追記
バグがありました。やけに左下に粒子が溜まると思ったら、粒子の位置を指定するxとyの値を整数にするのを忘れており、[x]と[x-1]で指定するセルが同じものになって左端からの移動ができていませんでした。直したところ、そこまで目立つ値は出てこなくなりました。
for x in range(N*N): #N×N回の試行
x = np.floor(np.random.random()*L)
y = np.floor(np.random.random()*L)
coin = np.random.random()*4
move(x,y,coin)