pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

トップ > Python > Pythonで学ぶ確率・統計（ランダムウォーク）

Pythonで学ぶ確率・統計（ランダムウォーク）

Python Python-応用 Python-応用-確率・統計

数学&物理学の具体的な計算にPythonを使って、数学&物理学もPythonも同時に学んでしまいましょう。今回はPythonを使った確率・統計の問題として、ランダムウォークの問題をみてみたいと思います。本記事のPythonの学習事項としては、離散データの作成・処理の方法、散布図の作成を学びます。ただし、動作の確認はJupyterNotebookを用いた場合で行っています。インターフェイスなどの動作環境の違いによって適宜変更点があるかもしれません。




コインの裏表（乱数）

Pythonには（擬似）乱数が用意されています：

import random as rd
rd.randint(0, 1)  # randint(a,b) で a 以上 b 以下の整数をランダムに出力

そこで、１を表に、０を裏に対応させることで、ランダムな出力をコイントスだと思うことができます。たとえば、コイントスを10回だけ試行する場合には

import random as rd
def r(k): return rd.randint(0, 1)   # ゼロか１をランダムに出力する関数
R = []
for k in range(10):
      R.append(r(k))
print(R)

とします。たとえば

\ \displaystyle{[1, 1, 1, 1, 0, 1, 0, 1, 0, 1]}

のような結果が出力されます。もちろんいま乱数を使っているので、出力される結果は入力のたびに変わります。実際に何度か同じコマンドを実行してみると、毎回出力が変わることが確認できると思います。

本記事は、以下の過去記事の発展的な内容になっています。あわせて読んでいただければさいわいです：

pianofisica.hatenablog.com

pianofisica.hatenablog.com


ランダムウォーク

いま考えているランダムウォークでは、コインの表（１）が出たときには粒子を１歩右側に、裏（０）が出たときには粒子を１歩左側に移動させるようにして、これを n ステップだけ繰り返したとき、粒子がどの位置にいるか？という問題を考えます。まずは10ステップ実行してみます

import random as rd
def r(k): return rd.randint(0, 1)
P = [0]   # 粒子の初期位置
def p(k): 
   if r(k) == 1: return P[k]+1   # １が出たとき１進む
   else: return P[k]-1   # １以外（ゼロ）が出たとき１戻る
for k in range(10):
      P.append(p(k))

このときの粒子の移動の様子を図示してみます：

import matplotlib.pyplot as plt
plt.plot(P)
plt.show()

最終的にどの位置に粒子がいるかは

P[10]

として確認することができます。問題のランダム性から、（初期条件としての０ステップ目以外の）ある時点で粒子がどの位置にいるかというのは再現性はなく、入力の都度その答えが変わってきます。

たとえば（リストの名前を変えて）もう一度同じコードを実行して、粒子の移動の様子を１回目の試行と合わせてプロットすると

P2 = [0]   # 粒子の初期位置
def p(k): 
   if r(k) == 1: return P2[k]+1
   else: return P2[k]-1
for k in range(10):
      P2.append(p(k))
plt.plot(P)
plt.plot(P2)
plt.show()

より、次のような図を得ます：

青線は１回目の試行時、オレンジ線は２回目の試行時です。


平均移動距離

ブラウン運動の理論

以上のような粒子のランダムな振る舞いに"法則性"というのがあるのでしょうか？じつは、このようなランダムな中にも法則があって、１次元的なブラウン運動の理論がこの場合に当てはまります。ブラウン運動の理論は、1905年のアインシュタインの『奇跡の年』に書かれた論文のうちの一報にさかのぼることができます（同年にはほかに、特殊相対性理論を基礎付ける論文、光電効果に関する論文なども出版されていて、まさに物理学史上における『奇跡』といわれるゆえんです）。ブラウン運動の理論からは、ある時刻 t に粒子をある位置 x に見出す確率分布関数 \rho(t,x) が次の拡散方程式

\qquad\displaystyle{\frac{\partial \rho}{\partial t}=D\frac{\partial^2 \rho}{\partial x^2}}

にしたがい、平均移動距離が

\qquad\displaystyle{\lambda=\left\langle\left(x-x_{0}\right)^{2}\right\rangle^{1 / 2}=\sqrt{2Dt}}

で与えられると予言されます。ここで x_0 は時刻ゼロでの粒子の位置（初期条件）です。いま離散化した問題を扱っているので、前者を直接的に議論するのは難しいですが、後者については統計的な平均をとって議論することができるでしょう。

統計的平均

例として100ステップをとるという試行を10,000回実行してみます：

import random as rd
t = 10000   # 試行数
N = 100   # ステップ数
def r(k): return rd.randint(0, 1)
P = [0]   # 粒子の初期位置
T = []
def p(k): 
      if r(k) == 1: return P[k-1]+1
      else: return P[k-1]-1
for l in range(t):
      for k in range(N):
         P.append(p(k))
      T.append(P[N])
      P = [0]   # 粒子を初期位置に戻す

最終的に粒子がどの位置にいるかをヒストグラムで表してみます：

import matplotlib.pyplot as plt
plt.hist(T,bins=100, range=(-30, 30))
plt.show()

から

のようなヒストグラムが作成されます。原点を中心とするきれいなガウス分布であることが一見してわかると思います。

平均（mean value）標準偏差（standard deviation）を求めましょう：

import statistics as st
m = st.mean(T)
s = st.stdev(T)
print('平均：　{}'.format(m))
print('偏差：　{}'.format(s))

から平均値（ \mu）と標準偏差（ \sigma）がそれぞれ

\qquad\displaystyle{\mu=0.018 \qquad\sigma=7.042}

と求まります（この値自体は再現性のあるものではありません）。この標準偏差が、法則

\qquad\displaystyle{\lambda=\left\langle\left(x-x_{0}\right)^{2}\right\rangle^{1 / 2}=\sqrt{2Dt}}

にしたがう平均移動距離を与えると考えられます。ここで係数 D については未知ですが、時間 t はステップ数が対応します。そこで、いま N_1=100 としているステップ数を N_2 に変化させたときの平均移動距離の比 {\lambda_{N_2}}/{\lambda_{N_1}} を考えて、これが

\qquad\displaystyle{\frac{\lambda_{N_2}}{\lambda_{N_1}}= \sqrt{\frac{N_2}{N_1}} }

と振る舞うかどうかをみることによって、法則の正否を確認することができそうです。

理論の実証

ステップ数 N_2 として N_2=n^2\times 100 \ (n=2,\,3,\,4,\,5) という４個の値を考えます。したがって理論の予言が正しければ、平均移動距離の比は

\qquad\displaystyle{\frac{\lambda_{N_2}}{\lambda_{N_1}}= \sqrt{\frac{N_2}{N_1}}=n }

と n に一致するはずです。これをみてみましょう（計算負荷の軽減のために、統計をとる試行数は各ステップ数で1,000回にしています）：

import random as rd
import statistics as st
S = 5   # n のラベルの最大値
d = 1   # n の刻み幅
t = 1000   # 統計のための試行数
N = 100   # N_1
def r(k): return rd.randint(0, 1)
def n(s): return int((1+s*d)**2*N)
def p(k): 
      if r(k) == 1: return P[k-1]+1
      else: return P[k-1]-1
def q(s): return st.stdev(T)
P = [0]
T = []
U = [0]
for s in range(S):
      for l in range(t):
            for k in range(n(s)):
               P.append(p(k))
            T.append(P[n(s)])
            P = [0]   # リスト P の初期化
      U.append(q(s))
      T = []   # リスト T の初期化

各ステップ数の場合で得られた標準偏差（平均移動距離）のリストは

U

で与えられます。ここで、上で作成したリストからグラフ（縦軸：平均移動距離、横軸：n）を描画してみます。

import matplotlib.pyplot as plt
plt.plot(U, marker="o", linestyle='None')
plt.show()

によって出力されるのが

のような図です。ほぼ直線状に並んでいることが観測できます。


これを１次関数でフィッティングすると

import numpy as np
X = [0,1,2,3,4,5]
f = np.polyfit(X,U,1)
y = np.poly1d(f)(X)
plt.plot(U, marker="o", linestyle='None')
plt.plot(y)
plt.show()

より

というグラフを得ます。ここで

f

から、線形近似の関数 y=ax+b が

\qquad\displaystyle{y=7.06\times n+0.05}

で与えられることがわかります。

また、 \lambda_{N_1=100} との比をプロットしてみましょう：

import numpy as np
V = []
for s in range(S+1):
   V.append(U[s]/U[1])
fr = np.polyfit(X,V,1)
yr = np.poly1d(fr)(X)
plt.plot(V, marker="o", linestyle='None')
plt.plot(yr)
plt.show()

から

のような図が得られます。ほぼ傾きが１の直線が得られ、関係式

\qquad\displaystyle{\frac{\lambda_{n^2\times N_1}}{\lambda_{N_1}}=n }

が成り立っていることがわかります。


数値計算の検証

10,000ステップ（ n=10 の場合）をとる試行を実行してその平均移動距離をみてみます：

import random as rd
t = 1000
N = 10000
def r(k): return rd.randint(0, 1)
P = [0]
T = []
def p(k): 
      if r(k) == 1: return P[k-1]+1
      else: return P[k-1]-1
for l in range(t):
      for k in range(N):
         P.append(p(k))
      T.append(P[N])
      P = [0]

より平均移動距離が

import statistics as st
st.stdev(T)

から

\qquad\displaystyle{\lambda_{N_2=10000}=72.5}

などと求まります。上の線形近似の予想からは

\qquad 7.06\times 10+0.05\sim 70.6

なので、上の関係が確かに成り立っていることがわかります。


今回は、乱数を使ってコンピュータ上でランダムウォークを実現してみました。ある現象について理論から予想される結論を、コンピュータを使って検証できるというのはなかなか面白いですね。





キーワード：Python、確率・統計、乱数、ランダムウォーク、ブラウン運動、データ作成、プロット

プライバシーポリシー