ブラウン運動についてのアインシュタインの関係式のシミュレーション
「ブラウン運動における微粒子の平均移動距離（出発点からの距離）が時間の平方根に比例する」というアインシュタインの関係式（Newton別冊「無とは何か、改訂第２版」、p. 60）のシミュレーションを行ってみた。なお、numpyのndarray型オブジェクトの計算例をこのサイトの後半に上げておいた。
ブラウン運動を次式のランダムウォーク（Random Walk）で表した。
ここで、出発点を原点とし、座標値の変化量
および
がランダムに独立に標準正規分に従うとする。すなわち、
とする。このとき、
となる。
式（２）に基づいて時刻0から時刻200までの50個の微粒子の動きを描画するスクリプトをリスト１のように作成した。
リスト１ ランダムウォーク（酔歩）のシミュレーション
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
T = 200
N = 50
for i in range(N):
Xs =
np.cumsum(ss.norm.rvs(size = T))
Ys =
np.cumsum(ss.norm.rvs(size = T))
plt.plot(Xs, Ys)
plt.title(f'N = {N}, T = {T}')
plt.xlabel('X', fontsize = 14)
plt.ylabel('Y', fontsize = 14)
plt.show()
リスト１のスクリプトを実行すると、図１のグラフが表示される。
図１ ランダムウォーク
リスト１のランダムウォークは、式（２）で与えられているので、時刻ｔにおける座標
は平均が０の独立な確率変数の和である。したがって、それらの平均
は0収束する（大数の強法則；P. Billingsley, 1995. Probability and Measure, third edition, John Wiley & Sons. p. 85）。
微粒子20個について
および
のグラフを
までの範囲で描くスクリプトをリスト２のように用意した。
リスト２ 平均座標値の収束
import scipy.stats as ss
import numpy as np
import matplotlib.pyplot as plt
T = 10000
N = 20
traces = []
for i in range(N):
Xs =
np.cumsum(ss.norm.rvs(size = T))
Ys =
np.cumsum(ss.norm.rvs(size = T))
traces.append([Xs, Ys])
traces = np.array(traces)
Ts = np.array([v+1 for v in range(T)])
mean_traces = traces / Ts
for i in range(N):
plt.plot(range(T),
mean_traces[i][0])
plt.plot(range(T),
mean_traces[i][1])
plt.plot([0, T-1], [0, 0], c='k')
plt.xlabel('t', fontsize = 14)
plt.ylabel(r'$\bar{X}_t, \bar{Y}_t$',
fontsize = 14)
plt.title(f'N={N}')
plt.tight_layout()
plt.show()
リスト２を実行すると、図２のようなグラフが得られる。
図２ 平均座標の収束（大数の強法則）
平均座標およびが0に収束していることが示されている。
アインシュタインの関係式では、出発点からの距離の平均値の時間変化が扱われる。時刻１から2000にわたる700個の微粒子の出発点からの距離の時刻ごとの平均値を描くスクリプトをリスト３のように用意した。
リスト３ アインシュタインの関係式
import scipy.stats as ss
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
from scipy.optimize import minimize_scalar as
minimize
T = 2000
N = 700
dist_traces = []
for i in range(N):
Xs =
np.cumsum(ss.norm.rvs(size = T))
Ys =
np.cumsum(ss.norm.rvs(size = T))
dist = (Xs**2 + Ys**2)**0.5
dist_traces.append(dist)
Ts = np.arange(1, T+1)
mean_dist_traces =
np.mean(dist_traces, axis=0)
plt.plot(range(T), mean_dist_traces,
label = 'Simulation')
def diff(a):
return
np.sum((mean_dist_traces - a * Ts**0.5)**2)
res = minimize(diff, bracket = [0,
100], method = 'Brent')
a = res.x
preds = a*(Ts**0.5)
plt.plot(range(T), preds, label =
'Model')
plt.xlabel('t', fontsize = 14)
plt.ylabel('Mean Dsitance', fontsize =
14)
plt.legend(fontsize = 14)
plt.show()
リスト３のスクリプトを実行すると図３のようなグラフが表示される。
図３ アインシュタインの関係式
シミュレーションによる距離の時刻ごとの平均値のグラフが青色の折れ線で描かれている。この折れ線に対して最小２乗基準で求めた時刻の平方根で表される関数
gが橙色の曲線で描かれている。
配列を扱う型としてnumpyのndarray型がよく用いられている。簡単な計算例を以下に挙げる。
Ndarray型と数値との演算は、ndarray型の各要素と数値との計算になるが、ndarray型同士のときは対応する要素ごとの演算になる。
例えば、次のスクリプト
b = numpy.array([1, 2, 3])
print('b =', b)
b5 = b * 5
print('b*5 =', b5)
a = numpy.array([1,10,100])
print('a =', a)
ba = b * a
print('b*a =', ba)
実行すると、次の出力を得る。
b = [1 2 3]
b*5 = [ 5 10 15]
a = [ 1
10 100]
b*a = [ 1
20 300]
内積の値は、関数numpy.dot()あるいは行列の掛け算演算子@により得ることができる。
次のスクリプト
dotba = numpy.dot(b, a)
print('dotba =', dotba)
b_at_a = b @ a.T
print('b_at_a =', b_at_a)
を実行すると、次の出力を得る。
dotba = 321
b_at_a = 321
numpyの関数はndarray型に適用すると、各要素に対して適用した結果が得られる。
例えば、次のスクリプト
pis = numpy.array([0.0, numpy.pi*0.5,
numpy.pi, numpy.pi*1.5, numpy.pi*2.0])
print('pis =', pis)
vsin = numpy.sin(pis)
print('vsin =', vsin)
を実行すると、次の出力を得る。
pis = [0.
1.57079633 3.14159265 4.71238898 6.28318531]
vsin = [ 0.0000000e+00 1.0000000e+00 1.2246468e-16 -1.0000000e+00
-2.4492936e-16]
上のランダムウォークの例では、配列の要素の累積和を求めるcumsumが使われている。
次のスクリプト例
a = numpy.array([1,10,100])
cumsum = numpy.cumsum(a)
print('a =', a)
print('cumsum =', cumsum)
を実行すると、次の出力を得る。
a = [ 1
10 100]
cumsum = [ 1
11 111]
ndarray型については、このウェブサイトでも説明している。
参考文献
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版