2017-08-12

体積粘性係数の分子動力学計算

Gromacs

等方的な流体の場合、粘性係数には、shear viscosityとbulk viscosity(体積粘性係数)の2種類がある。水のような液体の場合、非圧縮性近似をするのが一般的なので、Navier-Stokes方程式で体積粘性係数がかかる項は0になる。気体の場合も、流速が音速よりずっと遅い場合には非圧縮近似はよい近似になっているとされている。固体の場合は、弾性が支配的なので、通常、粘性そのものが考慮されない。

Wikipedia情報によれば、体積粘性係数の値は多くの場合、あんまり知られてないらしい

Volume viscosity/Measurement
https://en.wikipedia.org/wiki/Volume_viscosity#Measurement

The volume viscosity of many fluids is inaccurately known, despite its fundamental role for fluid dynamics at high frequencies. The only values for the volume viscosity of simple Newtonian liquids known to us come from the old Litovitz and Davis review, see References. They report the volume viscosity of water at 15 °C is 3.09 centipoise.

そもそも、通常の流体運動の範疇では体積粘性係数項は消えてしまって観測量と関わってこない(というか、NS方程式から消えてるだけだけど)から、測定が難しいのは仕方ない。多分、体積粘性係数を考慮に入れるべき身近な現象は音の減衰だと思う。なので、体積粘性係数を測定する方法として、音を利用するのは自然に思える。

Bulk viscosity and compressibility measurement using acoustic spectroscopy
http://aip.scitation.org/doi/abs/10.1063/1.3095471?journalCode=jcp

という論文には"To date, measurement of ultrasound attenuation is the only known way of experimentally studying bulk viscosity."と書かれている(Section Bの冒頭)。この論文によれば、25度に於ける水の体積粘性係数は、2.4cPであるらしい(TABLE III)

医療分野では、流体ではない生体組織の弾性や粘性(の分布)を測定するエラストグラフィと総称される手法がいくつか存在し、知る限り、いずれの方法でも超音波を何らかの形で利用している。
Elastography
https://en.wikipedia.org/wiki/Elastography
elastographyという言葉は、名前の通り、弾性の計測を主要な目的としている。癌などの腫瘍組織は一般に周囲の組織より硬くなる(乳がんが、しこりとして検出できる理由でもある)ことから、病変検出に有用と考えられている。一方、音の伝播の方程式に粘性項を入れるのは、固体であっても流体と変わるところはないので、同じ方法で粘性係数を決めることも理屈上はできる。実際に、そのような研究は存在するものの、弾性に比べると粘性を決定する医学的意義はよく分からないので、広く行われているとは言いがたいっぽい。しかし、何はともあれ、固体・流体を問わず、実験的に体積粘性係数を知る方法は、超音波の伝播や減衰を調べる以外にないらしい

現在は21世紀なので、簡単な分子の場合は、実験によらず、分子動力学で計算するという道がありえる。原理的には、shear viscosityも体積粘性係数も、Green-Kubo公式で計算できるはずなので。shear viscosityの計算は色々やられていて、手法もいくつか存在する

[1]Determining the Bulk Viscosity of Rigid Water Models
http://pubs.acs.org/doi/pdf/10.1021/jp211952y

は、分子動力学で、shear viscosityとbulk viscosityを決める論文。水のモデルとしては、TIP4P/2005がよいらしい。これは、他の論文でも、精度高いと書いてある

というわけで、Gromacsを使って、上の論文と同じことができないか試した。Gromacsのバージョンは5.1.2(aptで入れたら、これになった)。マニュアルには、gmx viscosityコマンドで、上の論文の手法が実装されてるよって書いてあるけど、そんなコマンドは存在しない。というわけで、手動で頑張ってデータを解析するしかない

また、水のモデルは、デフォルトでは、TIP4P/2005が入ってない

GROMACS files for the TIP4P/2005 model
http://www.sklogwiki.org/SklogWiki/index.php/GROMACS_files_for_the_TIP4P/2005_model

にあるのが、TIP4P/2005モデルらしいけど、面倒なので、普通のTIP4Pを使うことにした。

論文[1]によれば、

At the beginning, a MD simulation was performed at constant temperature and pressure conditions (NPT ensemble) for a period of 2 ns, in order to determine the density of liquid water as it is predicted by each water model.

The average energy ⟨E⟩ of the system at that density was computed from a second constant volume and temperature (NVT ensemble) simulation of 2 ns.

For the calculation of the bulk and shear viscosity at each temperature, we have performed a set of 100 independent constant volume and energy (NVE ensemble) simulations. For each trajectory, initial atomic momenta were sampled from a Maxwell distribution at temperature T.

らしいので、NPTアンサンブルで2ns→NVTアンサンブルで2ns→独立なNVE計算を100回(10psの平衡化の後、300psの計算をやると書いてある)という順番で計算するらしい。NVE計算では、長い時間やると、enerygy driftによって、エネルギー保存が成立しなくなるので短時間の計算を沢山やるということのよう。平衡化する場合、NVT→NPTの順が一般的な気がするけど、今回は逆

論文では、水256分子で計算しているので、適当に250分子を用意した。まあ、大体、以下のような感じのコマンドを叩く。エネルギー最小化は論文には書いてないので、特にやらなくてもいいかもしれない。NVE計算は、面倒なので100回やらずに30回しかやってない。

#-- We will get 250 molecules
gmx solvate -cs tip4p -o water.gro -box 2.1 2.0 2.0 -p conf/topol.top

#-- Energy minimization
gmx grompp -f conf/em1.mdp -o em1 -c water.gro -p conf/topol.top
gmx mdrun -s em1.tpr -deffnm em1 -v 
gmx grompp -f conf/em2.mdp -o em2 -pp em2 -po em2 -c em1 -t em1 -p conf/topol.top
gmx mdrun -s em2.tpr -deffnm em2

##gmx energy -f em1.edr -o min-energy.xvg
##gmx energy -f em2.edr -o min2-energy.xvg 

#-- NPT ensemble (to determine the density of liquid water)
gmx grompp -f conf/npt.mdp -o nptin -c em2.gro -p conf/topol.top -maxwarn 1
gmx mdrun -s nptin.tpr -deffnm nptout

#-- check density ,pressure, energy, temparature
echo -e "7\n8\n10\n14\n15\n" | gmx energy -f nptout.edr -o nptout.xvg

#-- NVT ensemble (to determine the average energy)
gmx grompp -f conf/nvt.mdp -o nvtin -c nptout -t nptout -p conf/topol.top
gmx mdrun -s nvtin.tpr -deffnm nvtout

#-- check energy
echo -e "7\n9\n11\n" | gmx energy -f nvtout.edr -o nvtout.xvg

#-- NVE ensemble
for i in `seq 1 30`;do
  gmx grompp -f conf/nve.mdp -o nvein -c nvtout -t nvtout -p conf/topol.top
  gmx mdrun -s nvein.tpr -deffnm nveout_$i
  echo -e "7\n8\n10\n" | gmx energy -f nveout_$i.edr -o nveout_$i.xvg    #--energy,temperature,pressure
done

使用したmdpファイルとtopol.topは、以下に置いた
https://gist.github.com/vertexoperator/aebb9676509e5eba608dfff6ff74bfc5

NPT計算の結果

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -9772.19        4.7    154.428    -16.766  (kJ/mol)
Temperature                 298.147      0.062    11.6618   -0.23701  (K)
Pressure                   -4.63437          2    868.369   -8.01756  (bar)
Volume                      7.51647     0.0023   0.124593 -0.00247449  (nm^3)
Density                     995.296        0.3    16.3923   0.194343  (kg/m^3)

とかなった。密度が大体正しそうなので多分OK。温度は298.15K=25度を指定しているので、そのへんで動く。論文[1]のTable Iでは、25度でTIP4Pを使用した時、0.9953(g/cm^3)という値を得たと書いてあるので、妥当な数値っぽい

NVT計算の結果は

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -9755.22        3.2    134.505    -16.928  (kJ/mol)
Temperature                 298.143       0.11    11.1177   -0.10815  (K)
Pressure                   -386.763        9.7    777.601    18.4497  (bar)

まぁ、よさげ？

Gromacsには、tcafコマンドというのがあって、私が理解してない謎の手法で、shear viscosityを計算できるらしい。
gmx tcaf
http://manual.gromacs.org/programs/gmx-tcaf.html

それで、NVT計算の結果とかに対して、

gmx tcaf -f nvtout.trr -s nvtin.tpr

とかやると、visc_k.xvgに(一部ヘッダ省略)

@    title "Fits"
@    xaxis  label "k (nm\S-1\N)"
@    yaxis  label "\xh\f{} (10\S-3\N kg m\S-1\N s\S-1\N)"
@TYPE xy
@    s0 symbol 2
@    s0 symbol color 1
@    s0 linestyle 0
 3.113 0.51114
 3.269 0.493288
 3.269 0.48944
 4.514 0.386852
 4.514 0.387317
 4.514 0.398109
 4.514 0.389984
 4.623 0.392635
 4.623 0.385249
 5.573 0.329574
 5.573 0.330219
 5.573 0.326714
 5.573 0.325943
 6.226 0.296062
 6.538 0.28417
 6.538 0.284974

などと出力される(デフォルトでは、visc_k.xvgにkとetaのみの出力がされる)。マニュアルによると、これを適当にfittingすると、shear viscosityが計算できるらしい

以下のようなコードで、実際にfittingしてみた

import numpy as np
from scipy import stats

k,y=[],[]
for line in open("visc_k.xvg"):
    if line[0]=="@":continue
    if line[0]=="#":continue
    ls = line.strip().split()
    k.append(float(ls[0]))
    y.append(float(ls[1]))


x=np.array(k)**2
y=np.array(y)
r=stats.linregress(x,y)
print(r.intercept)

計算すると、例えば0.540536022685という値を得たので、0.54e-3(kg/(m s))という感じであるらしい。水の粘性係数は、25度で0.000894(Pa s)らしいので、2/3くらいの値になっているが、論文[1]の方法では、TIP4Pで計算した場合、25度で、(4.83±0.09)e-4 (Pa s)という値になったらしいし、まぁTIP4Pの粘性に関する精度はこんなものであるらしい。NPT,NVEの結果から、tcafで計算しても、同じような結果を得た。TIP4P/2005なら、もっとよい結果が得られるらしい

Green-Kubo公式で、体積粘性係数を計算するために
$B_{ACF}(t) = \langle \delta P(t) \delta P(0) \rangle$
を計算する。
After an equilibration period of 10 ps, a configuration with total energy equal to the ⟨E⟩ decided before was selected. For that configuration and for the next 300 ps, the instantaneous stress tensor elements were stored in the disk every 20 fs for further analysis.
とか書いてある(はNVT計算の結果を使うらしい)けど、に到達するのに、ものすごく時間がかかるっぽかったので無視した(なんか勘違いしてる?)

よくやるように、
$B_{ACF}{(t) = \langle \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} dt' \delta P(t+t') \delta P(t') \rangle$
で計算する。これが得られたら、論文の式(4)でfittingを行う。次のpythonコードが実装例

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt


def func(t , C , tau_f , tau_s , omega , beta_f , beta_s):
   fast_term = (1-C)*np.exp(-np.power((t/tau_f) , beta_f))*np.cos(omega*t)
   slow_term = C*np.exp(-np.power(t/tau_s , beta_s))
   return (fast_term + slow_term)


BACF = []
for i in range(1,31):
   pressures = []
   for line in open("nveout_{}.xvg".format(i)):
        if line[0]=="@" or line[0]=="#":continue
        ls  = line.strip().split()
        assert(len(ls)==4)
        t,E,_,P = ls
        pressures.append( float(P) )
   Pbar = sum(pressures)/len(pressures)
   deltaP = np.array([P-Pbar for P in pressures])
   ACF = []  #--auto-correlation function
   for t in range(len(deltaP)):
      if t==0:
          sv = np.sum(deltaP*deltaP)/(len(deltaP))
      else:
          sv = np.sum(deltaP[t:]*deltaP[0:-t])/(len(deltaP) - t)
      ACF.append( sv )
   if len(BACF)==0:
      BACF = np.copy(ACF)
   else:
      BACF = BACF+ACF


dt = 0.02  #-- pico seconds
xdata = np.arange(len(BACF))*dt
ydata = BACF/30

popt , pcov = curve_fit(func , xdata , ydata/ydata[0])
plt.plot(xdata[:1000], func(xdata[:1000], *popt))
plt.show()

#--積分値(total relaxation time)
print( np.sum(func(xdata , *popt)) * dt)

あと、単位に注意する(粘性係数の単位をPa sとなるように合わせる)。上のdtはpsなので、10^{-12}(s)で、Gromacsが出力する圧力の単位は、bar(1bar = 10^5 Ps)。体積はNVT計算の平均値、温度も同様にして、計算する。上のコードで得られた諸々の数値をもとに、

(7.51647e-27/(1.38e-23 * 298.15))*(ydata[0]*1.0e10)*np.sum(func(xdata , *popt)) * dt*1.0e-12

で、求めたい体積粘性係数の値が出る。私の場合は、0.0018268112792585417という値を得た。論文のTABLE2によれば、298KでTIP4Pを使った場合、0.001173という値を得ている。オーダーは合っている。ずれは、アンサンブル数の少なさ、方法の若干の違いによるものと思われる。最初のほうで見た2.4cPという実測値と比較しても、悪くない。論文を信じれば、TIP4P/2005の方が精度はよいらしい。おしまい

2017-08-12

Boltzmannの壺とEhrenfestの壺

最近読んだ論文の読書感想文コーナー

[1]Statistical mechanics of money
https://arxiv.org/abs/cond-mat/0001432

N個の壺とE個のボールがあって、ボールはどれかの壺に入っているとする(論文には、ボールではなく、moneyと書いてあるし、壺とは書いてないけど、確率論に於いて、壺とボールはポピュラーなモデルなので壺にする)。(少なくとも一個以上のボールが入っている)壺をランダムに一つ選んで、やはりランダムに選んだ別の壺に移すということを繰り返すと、壺に入っているボールの平均数は当然E/Nだけど、Nがある程度大きい時、壺に入っているボールの数xの従う確率分布は、C exp(-x/T)で近似できる(Tは平均がE/Nになるように決める。Cは規格化定数）。ボール数をエネルギーと思い、壺と壺の間のボールの交換は、粒子と粒子の間のエネルギー交換の単純化と思うと、C exp(-x/T)はBoltzmann-Gibbs分布と見なせるというようなことが論文に書いてある。

統計力学のボルツマンの原理の類似を信じると、ボールがi個入っている壺の個数をn_iとした時
$N=\sum_{i=1}^B n_i$
$E=\sum_{j=1}^B j n_j$
という条件下で、
$W = \frac{N!}{n_1! \cdots n_c!}$
を最大化するような[n_1,...n_E]が、平衡状態では実現される。この条件から、Nが大きい時は、
$n_j/N \propto e^{-\beta j}$
であることが示せるというのは統計力学の教科書に書いてある通りで、dynamicsの話がなければ、まぁそういうもんかという感じであるものの、[1]の論文のように時間発展の仕方が定まれば、定常分布で巨視的状態[n_1,...n_E]が実現される確率は(原理的には)厳密に計算できる量となる(NとEが有限であれば、ただの有限Markov連鎖であるので)。その場合でも、(熱力学極限では)平衡状態が教科書通りに実現されるのかは自明でない気がする。で、[1]の論文にはこれが正しいことの説明があるが、数学の論文でなく、厳密な証明にはなってないと思うし、厳密な証明が知られているのかどうかは不明

このモデルを最初に考えたのは誰かと思ったら

非結合的代数によるランダムな衝突モデルの表現
https://ismrepo.ism.ac.jp/?action=pages_view_main&active_action=repository_view_main_item_detail&item_id=32615&item_no=1&page_id=13&block_id=21

によると、Boltzmann自身が、ほぼ同じモデルを考えていたっぽい。ただ、Boltzmannの結論だと $\sqrt{x} exp(-x/T)$ が出るらしい。ミクロカノニカル分布を仮定しても、この形は出ない。これは、本文中の条件(1.2)を仮定したためだと思うのだけど、ボールを一個ずつ交換する単純なモデルでは、これが成立する根拠はない。代わりに
$A_{\kappa \lambda}^{kl} = A_{kl}^{\kappa \lambda}$
を仮定して、Boltzmannと同様の議論を行えば、C exp(-x/T)は出る。これは良さそうな気もするけど、そもそも、このような計算が厳密には正しくないので微妙。例えば、N=4,E=10として、4つの壺に入っているボールの数(微視的状態)を(n0,n1,n2,n3)で表すとして、(1,2,7,0)->(0,3,7,0)と(1,2,3,4)->(0,3,3,4)の2種類の遷移を考える。1つめと2つ目の壺では、いずれも(1,2)->(0,3)という遷移が起きているけども、(1,2,7,0)から(0,3,7,0)へ遷移する確率と(1,2,3,4)から(0,3,3,4)へ遷移する確率は等しくない。従って、上のAがきちんと決まらない。

#元の論文では、ボールの数をmoneyと解釈しているが、現実の取引のモデルとしての妥当性は疑問(例えば、このモデルでは一旦金持ちとなっても、比較的容易に没落する)

一応、簡単な数値実験。

import numpy as np
from scipy.special import gammaln


N = 4  #--number of urns
urns = np.ones(N , dtype=np.int)*3
E = sum(urns)
p = {}
Q = {}          #-- N matrix
Nstep = 1000000

for i in range(Nstep):
    n0 = np.random.randint(0 , len(urns))
    if urns[n0]==0:continue
    n1 = np.random.randint(0 , len(urns))
    if n0==n1:continue
    bef = (urns[n0],urns[n1])
    urns[n0]-=1
    urns[n1]+=1
    aft = (urns[n0],urns[n1])
    #-- 本当は十分定常になってからcountすべき
    Q[(bef,aft)] = Q.get((bef,aft),0) + 1
    macro_state = tuple([sum(urns==b) for b in range(E+1)])
    p[macro_state] = p.get(macro_state,0) + 1



vtot = sum(p.values())
ctot = np.exp(np.sum(np.log(np.arange(N,E+N))) - np.sum(np.log(np.arange(1,E+1))))
w = np.zeros(E+1)
for k,v in p.items():
   t = np.array(k)
   c0 = np.exp(np.sum(np.log(np.arange(1,N+1))) - np.sum(np.where(t>0 , gammaln(t+1) , 0)))
   #--macrostate , observed probablitity , expected probability by microcanonical ensemble
   print( k , v/vtot , c0/ctot)
   w += t*(v/vtot)


for (s0,s1),c in Q.items():
    print (s0,s1,Q.get(((s1[1],s1[0]),(s0[1],s0[0])),0),c)

#AlderとWainwrightが1957年に始めた分子動力学計算(hard-sphere MD)では、衝突が起きるまで直進運動し、衝突時に保存則を用いて衝突後の速度を計算するというものであったらしい。MaxwellやBoltzmannが仮定していた分子運動論のイメージに近いのだと思う。この系は、決定論的でエネルギーも連続的であるので、Boltzmannの離散的なモデルよりは、現実寄りのモデルではあると思うけど、衝突によるエネルギー交換という単純なルールは保たれている

『Boltzmannの壺』モデルで、ある特定の壺に入っているボールの数の時間変化を追うと、時間平均はE/Nで、沢山サンプルを集めると、分布は、やはりC exp(-x/T)で近似できるようになるっぽい(数値実験の結果)。勿論、ボールの数はEが最大なので、NやEが有限の時は、C exp(-x/T)という分布は近似にしかならない。ところで、最初の分布は、集団に関するもので、後者は時間に関するものなので、この2つが同分布で近似できるというのは、"エルゴード仮説"っぽい(※)。これもやはり、E/Nを一定にして、NとEを無限大に飛ばした極限を考えると、この2つの分布のずれは0に収束していくと予想される。これも自明ではないと思うけど、厳密な証明は考えてないし、誰かが考えているかどうかも知らない

※)普通のエルゴード仮説は、任意の物理量の相空間上の空間平均と、相空間上の時間発展を追って時間平均を取ると等しいというものだけど、分子動力学屋は、よく、一粒子の物理量の時間方向の分布と、集団に対する分布が一致するという意味で、エルゴード仮説という用語を使っている気がする。MD屋のエルゴード仮説は、通常のエルゴード仮説よりも、強い主張になってる気がするけど、数値計算では成立することが多々あるらしい。Boltzmannのモデルでも、MD屋のエルゴード仮説が(近似的に)成立しているっぽい

#特定の壺のボール数の定常分布は、理屈上は、ボール数がk->k+1,k->k-1に変化する確率を求めれば厳密に計算できる。が、この計算は意外と面倒そうである(注目している壺にk個のボールが入っている時、残りの壺のいくつかは空であるかもしれないというのが面倒の原因)。別の直感的な考え方としては、Nが十分大きければ、2つの壺のボール数にほぼ相関はないと考えられる(例えば、N=2の時は、2つの壺のボール数には強い逆相関がある。Nが増えるにつれて、この相関が弱くなるのは、それほど非直感的ではないと思う)。一方、壺は、どれも特別でないので、各壺のボール数の定常分布は、どれも同じ分布に従うと考えられる。なので、集団としての分布は、独立同分布からサンプリングしたものとほぼ同じになるはず。定常分布に於いて、各壺にボールが(x_1,x_2,...,x_c)個入ってる確率をp(x_1,x_2,...,x_c)とすると、この確率分布は、exchangeable(x_1,...x_cの任意の置換によって値が不変)だから、De Finettiの定理によって、Nを無限大に飛ばした時、独立同分布に従うかのように振る舞い、この直感は正当化できそうな気がする(うさんくさい)

#等重率仮定を使った、よくあるタイプの計算「ある壺に入っているボールの数がxである確率は、残りの壺にE-x個のボールの数を分配する仕方の数W(E-x)に比例する」。例として、N=3,E=6の時を考えると、W(n)=n+1なので、p(x)=(7-x)/28が等重率原理から導かれる結果(平均はE/N=2になっている)。一方、定常分布を、厳密な遷移確率から数値的に計算すると、[p(0),p(1),p(2),p(3),p(4),p(5),p(6)] = [ 0.19047619, 0.25396825, 0.20634921, 0.15873016, 0.11111111, 0.06349206, 0.01587302]となって(検算として平均が2.0であることを確認)、当たり前だけど等重率仮定は全然成立していない。まあ、この例は、NやEが小さすぎて、そもそもBoltzmann-Gibbs分布で全然近似されてないので、問題とは言えないけど(ちなみに数値的に実験した限りでは、あるNまでは、p(0)p(1)と逆転するようである)

論文[1]には、エネルギー交換ルールの時間反転対称性を破ると、Boltzmann-Gibbsでない分布が出ると書いてある。私が最初同じモデルと誤解した、類似のモデルとして、multi-urn Ehrenfest modelがある。こっちは、毎回、ボールをランダムに一つ選んで取り出し、ランダムに選んだ壺に入れるというモデル。Boltzmannのモデルと違って、ボールをランダムに選ぶので、ボールが沢山入っている壺からはボールが取り出される確率が高くなり、少ししかボールのない壺からは、ボールが取り出されにくくなる。結果として、Boltzmannのモデルと違って、極端にボールが増えることもなければ減ることもない。十分時間が経った後は、特定の壺に入っているボールの数の時間分布は、B/Nの周りで正規分布(?)するっぽい。状態の集合は、Boltzmannのモデルと同じで、状態遷移の確率がBoltzmannのモデルと異なる。

multi-urn Ehrenfest modelは、元々は、N=2の時に、粒子の拡散の単純化したモデルとして提案されたらしい(N=2の時は単にEhrenfest modelという)。適当な容器内に沢山粒子があって、容器を2つの領域に分割し、各粒子がどっちの領域にいるかという問題を考えると、Ehrenfest modelとの類似性はイメージできる。

『Boltzmannの壺』とmulti-urn Ehrenfest modelの2つのモデルの違いはコードで書くと割と明確。

import numpy as np
import matplotlib.pyplot as plt

model_type = 0
N = 60  #--number of urns
urns = np.ones(N , dtype=np.int)*20
dist = []

for i in range(50000000):
    if model_type==0:
       n0 = np.random.randint(0 , len(urns))
    else:   #--multi-urn Ehrenfest model
       n0 = np.random.choice(np.arange(0,len(urns)) , p=urns/sum(urns))
    if urns[n0]==0:continue
    n1 = np.random.randint(0 , len(urns))
    if n0==n1:continue
    urns[n0]-=1
    urns[n1]+=1
    dist.append( urns[N//2] )  #--適当な壺の時間発展を追う


cnts = np.bincount(dist)
print( np.sum(np.arange(0,len(cnts))*cnts)/np.sum(cnts) )   #--print average

plt.plot(np.arange(0 , len(cnts)) , cnts)
plt.show()

Ehrenfest modelは、カットオフ現象というものが起きることが証明されてる。カットオフ現象は、多分1980年代に、Diaconisという数学者とその周辺の人々が発見した話らしい。multi-urn Ehrenfest modelにしても、最初のモデルにしても、有限Markov連鎖として定式化できる。有限Markov連鎖では、定常分布に到達する(定常分布に収束するとしてだが)のにかかる時間(混合時間)というものを考えることができるけど、定常分布に一様に近づいていくのかというと、必ずしもそうではなく、ある時間で急に定常分布に近付くということが起きるというのが、Diaconisらの発見らしく、カットオフ現象と呼ばれている(分布間の近さは全変動距離というものを使って測る)。

尤も、日常的な感覚からすると、ボールは区別することなく、また直接目にするのは、それぞれの壺に入っているボール数なので、それは、予想される通り、概ね単調に平均に収束していく。分布そのものの近さを測るということに意味があるのかどうかは、よく分からない

カットオフ現象を示すことが分かっているモデルはいくつかあって、最初に示されたのは多分、リフルシャッフルの数学的モデルと思う

Shuflling cards and stopping times
http://statweb.stanford.edu/~cgates/PERSI/year.html#86

マルコフ連鎖と混合時間ーカード・シャッフルの数理ー
http://www.kurims.kyoto-u.ac.jp/~kenkyubu/kokai-koza/H23-kumagai.pdf

2つ目の文献には、「一般のマルコフ連鎖について、カット・オフ現象が起こるための満足のいく必要十分条件は知られていないと言ってよい」と書いてある(数年前の話だけど、多分状況は変わってなさそう?)

Ehrenfest modelでカットオフ現象を示したのも、Diaconisらしい。ちなみに、どっちの解析でも、有限群の表現論が使われている。multi-urnだと、どうなるのか知らない
Asymptotic analysis of a random walk on a hypercube with many dimensions
http://onlinelibrary.wiley.com/doi/10.1002/rsa.3240010105/abstract

Ehrenfestモデルでは、"エントロピー"(正確には、定常分布との相対エントロピー)が単調増加することが厳密に証明されている。これ自体は、もっと広いクラスの有限Markov連鎖で成立するので、特筆すべき性質でもないけど、これはこれで全変動距離とは違う尺度で、定常分布への近さを測っていると見ることもできる。尺度を変えると、カットオフ現象とかがどう見えるのか(あるいは見えないのか)は知らない。ランダムに生成した連結なグラフ上のsimple random walkで実験した限り、KL divergence(負の相対エントロピー)の対数は、かなりキレイに直線的に減少していく(つまり、KL divergenceそのものは指数関数的に減衰する)(数値実験の詳細は下記コード参照)。こういう振る舞いをする理由が十分説明されているのかは不明

import networkx as nx
import numpy as np
import matplotlib.pyplot as plt


G=nx.fast_gnp_random_graph(950,0.02)
P=sorted(nx.connected_components(G), key = len, reverse=True)[0]

#--最大連結成分だけ残す
for c  in G.nodes():
   if not c in P:
       G.remove_node(c)


#-- transition matrix for simple random walk on G
W=np.zeros( (len(G.nodes()) , len(G.nodes())) )
for n , src in enumerate(G.nodes()):
   adjs = G.adj[src].keys()
   cnt = len(adjs)
   for tgt in adjs:
     m = G.nodes().index(tgt)
     W[n,m] = 1.0/cnt


#-- 定常分布の計算
eval,evec = np.linalg.eig(W.T)
idx = np.argmin(np.abs(eval-np.ones(len(G.nodes()))))  #--1.0に最も近い固有値のindex
pi = evec.T[idx]
pi /= np.sum(pi)

#--初期状態(任意)
p = np.zeros( len(G.nodes()) )
p[0] = np.random.random()
p[-1] += 1.0 - p[0]

#--適当なステップ数追跡
maxstep = 50
entropies=[]
for istep in range(maxstep):
    entropy = np.sum(p*np.where(p>0, np.log(p/pi), 0)) #-- KL divergence
    entropies.append(entropy)
    p = np.dot(W.T , p)



#--片対数グラフ
plt.plot(np.arange(0 , len(entropies)) , np.log(entropies) )
plt.show()