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()