■
Peter Woitという弦理論の批判で知られた人が、数年前に
Use the Moment Map, not Noether's Theorem
http://www.math.columbia.edu/~woit/wordpress/?p=7146
という記事を書いていた(私には要点がイマイチ分からなかったが)
ネーターの定理は、「保存量と対称性が対応する」という標語で知られる有名な定理で、普通、Lagrangianを使って定式化される。moment mapは、ネーターの定理のHamilton形式版というかsymplectic版といえる。数学では、Lagrangianそのものが使われない傾向にあるので、ネーターの定理を見かけることは、少ない。一方、物理学者が、moment mapを使っているのを見たことは、記憶にある限りではない。
物理学者は(数学者とは逆に)Lagrange形式を好んで使うこと、また、彼らが考えたい系の多くは可積分でなく、エネルギー、運動量、角運動量以外の"自明でない"保存量は、"some kind of charge"以外出てこないので、ネーターの定理で得られる保存量を重視する必要性が低い(そのような自明でない保存量は、Runge-Lenz vectorやKowalevski topで出てくるものがあるけど、教科書に載ってないことも多い)ことなどが理由として考えられる
moment mapを使うべき理由として、私が思いつくのは、以下のようなもの
(1)余接空間以外の相空間(≒一般のsymplectic多様体)でも適用可能(このような力学系は、特に数学では、しばしば扱われ、Lagrangianが不明である)
(2)個別の保存量や"対称性"ではなく、保存量や対称性の集合に存在する構造(つまり、Lie群やLie環)について記述する
(3)運動量写像は、"力学には依存しない"ことが定義から分かる
(4)運動量写像の量子化も、数学的に厳密に定義できる
#moment mapでは、時間並進対称性に対応する物理量がエネルギーだという事実が出ないように見えるけど、これは、moment mapの問題ではなく、相空間の取り方の問題だと思う。空間並進対称性に付随する物理量が運動量であることと、時間並進対称性に付随する物理量がエネルギーであることは、明らかに同じ構造を持っているので、通常の扱いでは、相空間が時間とエネルギーに対応する座標を持ってないことがトラブルの原因と思われる。
#時間・エネルギー座標を持つ相空間として、extended phase spaceというものが使われることがある。extended phase spaceは、ガリレイ群のmoment mapを定義しようと思えば必要だし、相対論を考えても不自然ではない。時間/エネルギー変数を含む正準変換を行いたいケースもある(以下の文献などで使われている)
Stackel系の全ての保存量を保つ離散化 (可積分系研究の新展開 : 連続・離散・超離散)
https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/42747
(3)について。ネーターの定理は、(少なくとも、見かけ上は)Lagrangianに依存するが、Lagrangianは、力学系ごとに異なる。一方、運動量写像は、symplectic多様体と、多様体への(symplectic形式を保つ)群作用のみで決まり、Hamiltonianには依存していない。従って、Hamiltonianを取り替えても、相空間と群作用が同じであれば、自明に同じ運動量写像を得る。
群作用が、Hamiltonianを保つなら、HamiltonianとPoisson可換な保存量を得ることができるけど、運動量写像の"型"は
で(Mはsymplectic多様体でetc.)、保存量はM上の関数でないと困る。"型"だけから類推してもpullbackを考えればよいというのは想像に難くなく、これはcomeoment mapと呼ばれることがある。comoment mapの"型"は
となる。は、Lie環の対称代数であるけど、Poisson代数になり、また、
に注意する。comoment mapは、群作用と可換で、Poisson括弧を保つ環準同型。環準同型なので、一次の元
の行き先が決まれば、残りも全部決まる。例えば、symplectic形式を保つ回転対称性SO(3)があれば、M上の関数L1,L2,L3で、{L1,L2}=L3,{L2,L3}=L1,{L3,L1}=L2を満たすものが得られ、これは、so(3)の基底のcomoment mapによる像である
#moment mapの定義によっては、保存量に定数を足す不定性が残ることがある。角運動量の例だと、L1,L2,L3をL1+c1,L2+c2,L3+c3に置き換えても(c1,c2,c3は異なる定数)、それぞれから定まる時間発展は変化しない。このような場合、moment mapのpullbackはPoisson括弧を保たない。moment mapの定義で、comoment mapがPoisson準同型になることを要請する場合もある。このへんは、上の(2)と関係する話
(4)について。いわゆるquantum moment mapも、comoment mapと対比すると分かりやすい。quantum moment mapの概念や名前を誰が最初に明文化したのかは分からない(概念自体は、誰でも思いつくようなものなので、多くの人が同時期に知っていたに違いない)。特に難しい概念ではないけど、日本語の教科書で見かけたことはないので、一応文献をあげておく
Lectures on Calogero-Moser systems
https://arxiv.org/abs/math/0606233
の25ページDefinition4.1を参照
#私は、昔誰かが"量子化運動量写像"と訳していたのを見て、その呼称を使っていたけど、今検索したら、このサイトしか出なかった。そもそも、普通に訳すと、"量子運動量写像"の気がする(これは、一件あった)。
quantum moment mapの"型"は、comoment mapの型に於いて、対称代数を普遍展開環に変え、関数環を非可換環に変えたものとなる(quantum comoment mapと呼んだほうが適切でないかとも思うが、非可換の場合の空間概念がないので、量子の場合、comomentの方しか存在しない)
これは環準同型で、群作用と可換になるという条件が付く。この場合も、一次の元(つまり、Lie環の元)の行き先だけ決まれば、残りは全部決まる。上のEtingofによるLecture Noteでは、大域的な変換ではなく、無限小変換を考えている。この定義だと、ネーターの定理は、ほぼ自明である。
#P,XをAの元として、ad(P)(x) = [P,x]と置くと、ad(P)は線形かつ非可換ライプニッツ則を満たす(ad(P)(QR) = Q(ad(P)(R)) + (ad(P)(Q))Rとなる)。 この条件を満たす元XをDer(A)としているが、X=ad(P)となるPがいつでも存在するとは限らないけど、存在すれば、正に保存量となる(逆に、保存量Pを与えれば、ad(P)によって、Der(A)の元が決まるので、非可換環Aへの無限小変換も自動的に定まることになる。従って、普遍展開環からの環準同型は、いつでも何らかの無限小変換に付随するquantum moment mapとなっている)。一方、群作用を微分して直接得られるのはDer(A)の元である
ここでの量子化の意味は、"filtered quantization"と呼ばれるものを考えるのがいい。一般に、フィルター環A(filtered quntizationではincreasing filterを考える)に対して、associated graded algebra(日本語訳は次数化環?)gr(A)を定義することができて、gr(A)は可換とは限らないけど、可換代数である時には、自然なPoisson構造が入る。このような時、Aはgr(A)のfiltered quantizationと呼ばれる(gr(A)のPoisson構造まで込みで一致する必要がある)。filtered quantizationは、多分、最近使われるようになった用語で、あまり浸透していないっぽい(検索すると、arXiv:1212.0914やarXiv:1704.05144で見つかる)。
名前が付いてなかったものの、かなり昔から、界隈の数学者は、この意味で(暗黙のうちに？漠然と？)量子化を捉えていたっぽい(上記EtingofのLecture Noteでも、filtered quantizationという用語は使われていないが、正にこの意味での量子化になっていることが分かる)。この用語に従えば、Lie環の普遍展開環は対称代数のfiltered quantizationであり、Weyl代数は、上の多項式関数環のfiltered quantizationである
filtered quantizationで得たフィルター環から完備なRees代数を構成すれば、変形量子化になる。1990年以降に見つかってきた新しい例も多く、例えば「有限W代数はSlodowy sliceの量子化である」という時も、filtered quantizationの意味での量子化と理解できる(arXiv:math/0105225)。また別の例としては、(ある種の?)3-Calabi-Yau代数などがあるらしい(よく知らない)
Noncommutative del Pezzo surfaces and Calabi-Yau algebras
https://arxiv.org/abs/0709.3593
Calabi-Yau algebras viewed as deformations of Poisson algebras
https://arxiv.org/abs/1107.4472
associated graded algebraは、可換にならない場合でも、例えばClifford代数からは外積代数が得られる。外積代数には、Poisson superalgebraの構造が入り、『Clifford代数は外積代数の量子化である』と言われるのも、filtered quantizationの一種として理解できる。
#量子化すると、非自明な中心拡大が生じることがある(ガリレイ代数の質量項や、Virasoro代数のcentral chargeなど)。これらは、filtered quantizationとは別の枠組みで理解されるべきようにも見えるけど、よく分からない
(3)に戻ると、運動量写像にHamiltonianは必要ないので、例えば、相空間がsymplectic等質空間であれば、変換がHamiltonianを不変に保つか否かに限らず、運動量写像を考えることができる。そして、量子論では、このような運動量写像の量子化によって、spectrum generating algebraが作られ、波動関数に作用する。このような観察から得るべき教訓は、Hamiltonianを不変にするような"力学的対称性"だけでなく、相空間の幾何構造を不変にする幾何学的対称性に付随する物理量にも目を向けたほうがいいということだと思う
私は、普通のmoment mapで考えると混乱するので、comoment mapしか使わない(computer algebraを利用しやすいように、物事を代数的に考えておきたいというのもある)けど、comoment mapは、保存量が直接得られる点とquantum moment mapとの自明な類似性の二点に於いて、moment mapより優れていると思う。なので、「ネーターの定理より、moment mapより、comoment mapを使おう」といえる(数学では、moment mapの像を見たいこともあるけど)
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_1,...n_E]が、平衡状態では実現される。この条件から、Nが大きい時は、
であることが示せるというのは統計力学の教科書に書いてある通りで、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の結論だとが出るらしい。ミクロカノニカル分布を仮定しても、この形は出ない。これは、本文中の条件(1.2)を仮定したためだと思うのだけど、ボールを一個ずつ交換する単純なモデルでは、これが成立する根拠はない。代わりに
を仮定して、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)
論文[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()
体積粘性係数の分子動力学計算
等方的な流体の場合、粘性係数には、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公式で、体積粘性係数を計算するために
を計算する。
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.
とか書いてある(
よくやるように、
で計算する。これが得られたら、論文の式(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年3月の『Cell』の論文
Targeted Apoptosis of Senescent Cells Restores Tissue Homeostasis in Response to Chemotoxicity and Aging
http://www.cell.com/fulltext/S0092-8674(17)30246-5
大雑把には、(マウスの)"老化細胞"(Senescent cells)でアポトーシスを誘導してやって除去すると、若返るよという話(例えば、Fig7(D)では、毛が明らかに増量している)。ここ数年、老化細胞の除去による"若返り"の報告はいくつかあるが、いずれも老化細胞の除去に遺伝子改変を必要とするものばかりだったのに対して、この論文の方法では、遺伝子改変を必要としない(マウスの実験であるから、ヒトなどでは、別途検討が必要ではあるが)
cf)細胞老化が生体機能に及ぼす影響とその機構に関する研究
http://www.ncgg.go.jp/ncgg-kenkyu/documents/25-18.pdf
http://www.ncgg.go.jp/research/news/20160810.html
肺組織で特異的に老化細胞を除去できるように、遺伝子改変を行ったマウスでの類似実験を行っている
cf)Naturally occurring p16Ink4a-positive cells shorten healthy lifespan
https://www.nature.com/nature/journal/v530/n7589/full/nature16932.html
背景。正常な細胞は、老化刺激(テロメア短小化、DNA損傷、酸化ストレスetc)を受けると、恒久的に細胞周期を停止した細胞老化と呼ばれる状態になることが知られている。こうしてできた老化細胞は、増殖を再開することはないと考えられているが、ただ機能を失って存在するだけでなく、炎症性サイトカイン、増殖因子、細胞外マトリックス分解酵素など、様々な分泌因子を分泌するSASPと呼ばれる現象を起こし、周辺組織に炎症状態を誘発する。SASPは、周囲の細胞で細胞老化や癌を誘導するなど生体にとって有害な作用をもたらすと現在考えられている。
cf)老化した細胞ががん化を促進する仕組みをハエで解明(プレスリリース)
http://research.kyoto-u.ac.jp/research/141011_1/
論文では、老化細胞でFOXO4の発現が亢進することを確認し(Fig1(D))、また、FOXO4を阻害すると、老化細胞の生存率が下がる(Fig(1(H))ことから、FOXO4が老化細胞のアポトーシスを抑制し、老化細胞の生存・維持に寄与していると推測し、FOXO4のp53結合部位を含む部分ペプチドのD-retro inverso (DRI)-isoform(FOXO4-DRI)を設計し、マウスに腹腔投与(i.p. injection)したという感じ(5mg/kgを一日一回x3回)。
#D-retro inverso-isoformというのは初耳だったのだけど、配列とアミノ酸のキラリティ(つまりD型アミノ酸から構成される)が逆になってるような異性体のことらしい。なんだかよく分からないけど、DRI-isoformはしばしば、元のペプチドと同じ活性を保持し、かつ、安定であるらしい(多分プロテアーゼなどで分解されにくい?)。今の場合は、p53とFOXO4の結合を競合阻害するのが目的だろうと思う(部分ペプチドなので、p53を抑制する機能はない)
投与を行った"老いたマウス"の年齢は104〜110週齢程度のよう(Fig7 captionなど)。20ヶ月齢のマウスが人間でいうと60歳相当程度らしいので、人間換算で30歳弱と若干若い感じがする。マウスを直接扱ったことがないので、何とも言えないけど、老化研究では、もう少し高齢のマウスを対象にすることが多いように思う。一応、Fig6(B)の腎臓組織のSA-β-GAL染色では、130week miceでも、結構、老化細胞が増えているらしきことが見て取れるし、最初の方のcfであげた
細胞老化が生体機能に及ぼす影響とその機構に関する研究
http://www.ncgg.go.jp/ncgg-kenkyu/documents/25-18.pdf
でも、12ヶ月齢のマウスが使われている(これは、肺組織に於ける老化細胞の除去実験)。
毒性について。FOXO4は、FOXO1とFOXO3に比べると、注目されることが少なく、非老化細胞ではあまり重要な役割を持っていないらしい。実際にFOXO4をノックアウトしたマウスは、目立った表現型を示さない("FOXO4 knockout mice do not show a striking phenotype")らしいので、FOXO4の阻害は、非老化細胞に対する影響は少ないと予想される。
[疑問]FOXO4の阻害が老化細胞のアポトーシスを引き起こすなら、FOXO4ノックアウトマウスでは、老化細胞の数が少なく(論文の結果が正しければ)見た目の老化も抑制されそうに思えるが、目立った表現型は示さないという報告とは矛盾する気もする。老化が遅いなどの表現型は注意して観察しなければ気付かない微妙なもの(致死とかに比べればstriking phenotypeでない)で単に見落とされているという可能性もあるけど
[補足]癌抑制遺伝子として有名な転写因子p53は、細胞周期を停止し細胞増殖の抑制を行うだけでなく、細胞老化を引き起こし、アポトーシスを誘導することも知られている。これらの機能は、独立しているわけでなく、細胞やDNAがダメージを受けた時に細胞周期を停止し修復を試みると同時に、ダメージが重度で(癌化しそう)であればアポトーシスを起こし、中度であれば細胞老化を引き起こすのだろうと思われる(回復可能な程度に軽度であれば、細胞周期を再開する)。というわけで細胞老化を誘導する薬剤を高濃度で与えると、アポトーシスを誘導するっぽい。培養細胞に於いて、p53の発現を亢進すれば、癌は減るが老化は進行し(寿命も縮むらしい)、一方、抑制すれば、老化は抑えられるが癌は増えるらしい(マウス個体でp53をノックアウトすると早期に死亡するので解析は難しいが、培養細胞での実験と矛盾しない結果が知られている)。このような癌化と細胞老化のトレードオフは、テロメラーゼの高発現/ノックアウトでも観察することができる(マウスはテロメアが長く、テロメラーゼノックアウトマウスで個体老化が観察できるようになるまで数世代かかる)らしい
cf)癌抑制遺伝子p53 による老化の制御
http://www.pharmacol.or.jp/fpj/topic/topic121-2-130.htm
#細胞周期の一時停止、細胞老化の誘導、アポトーシスの選択がどのようにされているのか詳細は不明。p53には、リン酸化部位が沢山あり、中でもSer46がリン酸化される(このリン酸化を行う酵素はHIPK2,ATK kinase,DYRK2など複数ある?)と、アポトーシス関連遺伝子(AIP1という遺伝子が知られる)の発現が誘導されるらしい。一方、Ser15のリン酸化(AMPKなどによる?)は細胞老化関連遺伝子の発現と関連するらしい。DNA損傷が起きた時、Ser46のリン酸化はSer15のリン酸化に比べると、少し後になってから起こるらしい
The p53 response to DNA damage.
https://www.ncbi.nlm.nih.gov/pubmed/15279792
Ser46 phosphorylation of p53 is not always sufficient to induce apoptosis: multiple mechanisms of regulation of p53-dependent apoptosis.
https://www.ncbi.nlm.nih.gov/pubmed/17584297
Palmdelphin, a novel target of p53 with Ser46 phosphorylation, controls cell death in response to DNA damage
http://www.nature.com/cddis/journal/v5/n5/full/cddis2014176a.html?foxtrotcallback=true
素朴な疑問として、何のために細胞老化などという機構があるのか。癌抑制のためだけであれば、とっととアポトーシスを起こすなりして除去されてもよさそうなものだけど、わざわざアポトーシスを抑制してまで生存を図る意味が分からない。最近では、創傷治癒に於いて細胞老化の果たす役割が報告されたりしていて、他にもまだ知られてない生物学的機能が存在する可能性はある
cf)Cellular senescence controls fibrosis in wound healing
https://www.ncbi.nlm.nih.gov/pubmed/20930261
それらの話とはまた別に組織幹細胞の減少と機能低下も個体老化の一因として挙げられることがある。細胞老化と幹細胞は、共に創傷治癒に寄与する、という共通点があるっぽい。7月のNatureに
Hypothalamic stem cells control ageing speed partly through exosomal miRNAs
http://www.nature.com/nature/journal/v548/n7665/full/nature23282.html
という論文が出ていた。これは、視床下部神経幹細胞(htNSC)の減少が老化のトリガーとなっており、特に、exosomal miRNAの放出量が減ってしまうことが主要因なのでないかとのこと(生物種はマウス)。
ちなみに、論文中では"mid-aged mice (11–16 months old)" "aged mice (22 months and older)"と書かれていて、やっぱCell論文の110wk miceとかは中年手前なんだな、という感じ。こっちの論文のFig1(a)などを見ると、16 months-old miceでも、視床下部幹細胞の減少が見て取れるし、高齢マウスでは殆ど消失しかかっている。Cell論文に比べると、視床下部幹細胞の注入("炎症"のせいか、単純に移植しても、すぐに死んでしまうので、ドミナント・ネガティブIκB-αを発現する細胞を作ったと書いてある)などにより、老化の進行が抑制されて寿命が延長してる(16ヶ月齢雄マウスに移植して4ヶ月後の生理機能(Fig3(c))や18ヶ月齢雄マウスに移植して寿命計測(Fig3(d))
この論文では、視床下部幹細胞に注目しているが、
Muscle-derived stem/progenitor cell dysfunction limits healthspan and lifespan in a murine progeria model
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3272577/
のような論文も以前に出ている(タイトルを直訳すると、『筋由来幹/前駆細胞の機能不全が、早老症モデルマウスの健康寿命と寿命を制限する』)。若いマウスのmuscle-derived stem/progenitor cellを移植することにより、早老症モデルマウスの寿命を大幅に延長できたようである(Fig4(a))。通常の老齢マウスに対して、同様の実験を行って寿命延長などが起こるかは書いていない。まぁ、視床下部幹細胞が特別というわけでもなさそうではある。また、移植後、比較的短期間で効果が出ていることから、移植した幹細胞そのものが分化しているわけではなく、何らかの分泌因子が重要らしい(が、分泌因子の正体は不明)というところで、終わっている。視床下部幹細胞の方も事情は似ていて、移植した細胞が直接分化するわけではなく、こちらは、何らかのmiRNAが重要らしいということになっている
物凄く単純化していうと、老化細胞除去は、"悪いもの"を減らすという考え方で、一方、若い幹細胞を移植するというのは、"良い物"を増やすという考え方なので、割と真逆ではある。もし後者が重要であるならば、単純に老化細胞を除去しても、機能低下した幹細胞が回復するとは限らず、高齢になってからの老化細胞除去は、効果が薄いか、手遅れである可能性もある。老化細胞除去が比較的若いマウスに対して行われている理由かもしれないけど、どうだろう。
ノンパラメトリックベイズ(1)infinite GMM
元論文
[1]The Infinite Gaussian Mixture Model
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.9111
論文中に"The likelihood for α may be derived from eq. (12)"とあるけど、何度読んでも、式(15)が式(12)から出るというのは無理がある気がするし、そもそも式(15)の1つ目の式は多分正確ではない。
おそらく、式(16)のChinese Restaurant Process(CRP)から出すのが正しい
具体的には
で、j番目のクラスに割り当てられた要素の数を$n_j$として、そのindexを昇順にI(1,j),...,I(n_j,j)とする。
と略記すると、クラスの数を$k$として
と積を並び替えられる。そして、CRPの定義(論文中の式(16))から
と計算できる。最終的に
という分布を得る。これは、Ewens分布と呼ばれるものになっている。
で、この分布に従って、n個の要素をk個のクラスに分割する確率は
となる。但し、S(n,k)は第一種Stirling数で、Γ(n+x)/Γ(x)をxのn次多項式と見た時のk次の係数に等しい。これが、式(15)の1つ目の式の正しい形と思う。第一種Stirling数の定義から、
が容易に分かる。nとkを固定してしまえば、式(15)の2つ目の式は、正しいので、論文通り、実装しても問題は生じない
逆に、Ewens分布からCRPを導出するのも、条件付き確率分布の定義通り計算すればいいので容易。
Ewens分布は、exchangeableなランダム分割のモデルの一つで、他の例として、Ewens分布の1パラメータ変形であるEwens-Pitaman分布は、Pitman-Yor過程とかで使われるらしい。exchangeableという性質は、上のEwens分布で、c_1,...,c_nの順番を任意に置換しても結合確率が変わらないことを指している。exchangebilityは、今はGibbsサンプリングが使えるために要求される性質と思っておけばいいっぽい。
cf)Exchangeable random variables
https://en.wikipedia.org/wiki/Exchangeable_random_variables
原理的には、exchangebleなランダム分割のモデルに対してinfinite GMMができるのだと思う(というわけで、Pitman-Yor mixture modelというものもあるらしい)。infinite GMMに於いて、Ewens分布を選ぶ理由は特になさそうに見えるけど
Exchangeable Gibbs partitions and Stirling triangles
https://arxiv.org/abs/math/0412494
によれば、Ewens-Pitman分布は、exchangeableなランダム分割としては、かなり一般的なモデルとなっているらしい
#論文[1]と同じようにして、Pitman-Yor mixture modelを実装してみようかと思ったけど、hyperparameterの更新式が複雑になった。これをどう乗り越えるべきか確信がなかった(力技で計算する以外の方法が思いつかなかった)ので、今回はpending
[余談]ちゃんと読んだわけではないけど、以下の論文2.4.2節によれば、Ewens分布は、Jack測度というものと同一視できるらしい(Jack測度は、自然数nの分割上で定められた測度)。論文で扱われている測度は、Macdonald測度とでも呼ぶべきだろうか(一般的に通りのいい名前は付いていないっぽい？)。Schur多項式->Jack多項式->Macdonald多項式という(1パラメータ、2パラメータ)変形に対応して、Plancherel測度->Jack測度->"Macdonald測度"という変形があるっぽい。
A probabilistic interpretation of the Macdonald polynomials
https://arxiv.org/abs/1007.4779
[余談2]arXiv:hep-th/0306238によれば"In many aspects, the set of partitions equipped with the Plancherel measure is the proper discretization of the Gaussian Unitary Ensemble (GUE) of random matrices"だそうである。Jack多項式のαは任意パラメータであるが、α=1,2,1/2の時は、ある対称空間上の球関数となる。一方、ランダム行列側では、ガウス型アンサンブルとして、GUE,GOE,GSEが知られている。α=1のJack測度がPlancherel測度で、その適当な極限は、GUEランダム行列の固有値の極限分布に一致するらしい(arXiv:math/9903176,arXiv:1402.4615)。というわけで、Ewens分布は数学的に見ても、非常に素性のよさそうな分布である
#ランダム行列について何も知らないけど、"hermitian random matrix ensembles"については、ある種の対称空間のCartan classとの対応があって、GUE/GOE/GSEは、その一例であるらしい(arXiv:cond-mat/0304363,arXiv:0707.0418)。今やランダム行列ユーザーの最大勢力は物理屋だと思うけど、これらのrandom matrix ensemblesは、class Bとclass DIII-oddを除けば、トポロジカル絶縁体・超伝導体の分類にも現れる(arXiv:0905.2029)。class Bとclass DIIIが物理で現れないというわけでもないっぽい(arXiv:cond-mat/0103089,arXiv:cond-mat/0103137)けど、ハミられている理由はよく分からない。最近は、"non-hermitian random matrices"も、応用があるそうな。
もうひとつ、論文[1]の式(17)の後半で、新しいクラスを割り当てる確率の計算で、積分が必要になるけど、
[2]Markov chain sampling methods for Dirichlet process mixture models
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.115.3668
のAlgorithm8を使って積分を回避する。この論文には"However, the equilibirum distribution of the Markov chain defined by Algorithm 8 is exactly correct for any value of m"とあるのだけど、何故か分からない(mが大きければ、積分のモンテカルロ評価になるというのは論文に書いてある通りで明らかだけど)。このあとのステップで、同じ基底分布G0からφ_cをサンプリングすることになるので、そうなりそうな気もするけど...
全体の流れとしては、初期化が終わったら、以下の処理を繰り返す感じになる
(1)クラス割当の更新。論文[2]のalgorithm8を使う
(2)各クラスのcomponent parameters(今の場合、個々のクラスは、単一のガウス分布に対応しており、その平均と分散)をサンプリング。論文[1]の式(4)と(8)を使う
(3)hyperparametersの更新。論文[1]の適当な分布に従ってサンプリング。更新するhyperparameterは大きく分けると、基底分布G0のパラメータ(4つある)とDirichlet過程の集中度αの二種類となる
これを十分に繰り返せば、クラス割当/component parameters/hyperparametersの分布が得られるが、一番尤もらしい割当とパラメータが欲しいだけなので、それぞれの尤度を計算して一番いいものを選ぶ。十分多く繰り返せば、MAP推定値に近づくはず。
尤度は、DPGMMで対象データが生成される確率を計算する。ハイパーパラメータを決めれば、DPGMMに従うデータ点を生成していけるが、そのようにして実際に生成される点は、元の学習データとは、あまり似てないものになる公算が高そうな気がするけど、大丈夫なんだろうか。各クラスターの平均(と分散)は、基底分布G0に従って生成するので、特にクラスター数が少ない時には、元の学習データにあるクラスターとずれたところにクラスターを生成しても不思議ではなiい。そのへんは、通常の混合ガウス分布の方が、ずっと有能そうではある
#infinite GMMやDPGMMはノンパラメトリックベイズ的手法の一つと言うけど、ノンパラメトリックという言葉は全く正しくないように思うし、誤解と混乱を招く言い方の気がする。infinite GMMで使われているDPGMMのパラメータ(上でhyperparameterと書いてるもの)はαと基底分布のパラメータからなる。データの次元が一次元の場合は、基底分布のパラメータは4つ(正規分布の分が2つとガンマ分布の分が2つ)なので、合計5つのパラメータを与えれば、モデルを指定できる。一方、(一次元の)混合ガウス分布の場合は混合数をkとすると、平均と分散がk組、混合比がk-1個のパラメータで指定されるので、3k-1個のパラメータで定まる。一次元では、DPGMMは、混合数2の混合ガウス分布と同じ数のパラメータを持つ。
現在では、MCMC以外の、より高速な推定手法も提案されてはいるようだけど、特に調べてはいない。
問題として、βの分布が、k=1の時、log-concaveではない気がする(ので、ARSが失敗する場合がある)。以下のコードで、対数確率密度関数のplotをしてみたけど、k=1の時は、xが2か3あたりを超えると、微妙に凹でない
import matplotlib.pyplot as plt from scipy.special import gammaln s,w=np.array([ 0.69756949 ]) , 1.159287098021641 k = len(s) x = np.arange(1 , 100)*0.1 y = -k*gammaln(x/2) - 1/(2*x) +((k*x-3)/2)*np.log(x/2) + (x/2)*np.sum(np.log(s*w)) - (x*w/2)*np.sum(s) plt.plot(x, y) plt.show()
αも同様に、log-concaveでないと思う(プロットして見ても、かなり分かりにくいが、kが小さい時の方が、"凹度"は高い気がする)
import matplotlib.pyplot as plt from scipy.special import gammaln k,n=1,500 x = np.arange(1,101)*0.5 y = (k-3/2)*np.log(x) - 1/(2*x) + gammaln(x) - gammaln(n+x) plt.plot(x, y) plt.show()
(私が間違ってるのでないとすると)結構イイカゲンな論文だな、という印象なんだけど、どうなんだろう
以下に実装したコードを置いておく。面倒なので、データが一次元空間上に分布している場合のみを考えている。αとβのサンプリングは、論文に従ってARSで行っているものの、上で書いた理由で、正しくはない(し、エラーを起こす場合がある)。今回の目的はサンプリングではなく、MAP推定なので、αとβのサンプリングが多少悪くても、よかろうと思って、そのままにしてある
実際に、DPGMMでデータを生成して、infinite GMMでfittingして見ると、データが一次元のせいもあって、クラス間の分離が悪い場合は、本来存在するより少ないクラスしか出てこない場合があるように思える。対数尤度を計算してみると、infinite GMMで計算した方が高くなっているので、実装が間違ってるわけでもないと思う
まぁ、さしあたって、infinite GMMは、ノンパラメトリックベイズの一番簡単な例題として勉強しただけで、本来の目的ではないので、とりあえず、諸々の問題点・課題(α、βのサンプリング、高次元データ対応etc.)は、そのままにしておいて、次に進む
"""
This is an implementation of the following paper written in python 3.5
(The only univariate case is supported)
The Infinite Gaussian Mixture Model
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.9111
"""
import numpy as np
import math
from scipy.special import digamma
def rnd_from_piecewise_exponential(coeffs , intercepts , points):
"""
density function h(x) in [points[i] , points[i+1]]
h(x) = exp(coeffs[i]*x + intercepts[i])
"""
assert(len(coeffs)==len(intercepts)),(coeffs,intercepts)
assert(len(coeffs)+1==len(points)),(coeffs,points)
u = np.random.random()
CDF = [0.0] #--cumulative density value at each points
for i in range(len(points)-1):
z_cur = points[i]
z_next = points[i+1]
assert(z_cur<z_next),(points)
a = coeffs[i]
b = intercepts[i]
if a!=0.0:
r = np.exp(b)*(np.exp(a*z_next) - np.exp(a*z_cur))/a
else:
r = (znext-z)*np.exp(b)
assert(not np.isinf(r)) , (coeffs , intercepts , points)
r0 = CDF[-1]
CDF.append( r + r0 )
u = CDF[-1]*u #--CDFが正規化されてないので
for i in range(len(points)-1):
if CDF[i] < u and u <=CDF[i+1]:
a = coeffs[i]
b = intercepts[i]
z = points[i]
"""
compute t such that
( exp(a*t+b) - exp(a*z+b) )/a == u - CDF[i]
"""
if a!=0.0:
t0 = a*(u - CDF[i]) + np.exp(a*z+b)
assert(t0 > 0.0),t0
t = ( np.log(t0) - b )/a
else:
t = (u - CDF[i])/np.exp(b) + z
assert(t >= z and t <= points[i+1]),(t,z,points[i+1])
return t
assert(False),("rnd_from_piecewise_exponential",u,CDF,points,coeffs,intercepts)
def ars(h , hprime , points , support=(0 , np.inf)):
xs = [x for x in points]
xs.sort()
assert(len(xs)>0)
assert(xs[-1]>0)
while hprime(xs[-1])>=0.0:
xs.append( 2*xs[-1] )
while 1:
xs.sort()
x0 = xs[0]
xN = xs[-1]
zs = [support[0]]
coeffs = []
intercepts = []
for i in range(len(xs)-1):
x = xs[i]
xnext = xs[i+1]
assert(xnext>x)
zi = (h(xnext) - h(x) - xnext*hprime(xnext)+x*hprime(x))/(hprime(x) - hprime(xnext))
assert(zi>zs[-1]) , (zi,zs,xs , xnext,h(x) , h(xnext) , hprime(x) , hprime(xnext))
zs.append( zi )
coeffs.append( hprime(x) )
intercepts.append( h(x) - hprime(x)*x )
zs.append( support[1] )
coeffs.append( hprime(xN) )
intercepts.append( h(xN) - hprime(xN)*xN )
y = rnd_from_piecewise_exponential(coeffs , intercepts , zs)
w = np.random.random()
ukval = 0.0
assert(not np.isinf(y)),(coeffs,intercepts,zs,y)
for i in range(len(zs)-1):
if zs[i]<=y and y<=zs[i+1]:
ukval = h(xs[i]) + (y-xs[i])*hprime(xs[i])
break
if w<=np.exp(h(y) - ukval):
return y
else:
xs.append( y )
def gen_alpha(k , n):
def pdfln(x):
C = sum([np.log(i) for i in range(1,n+1)])
return C + (k-3/2)*np.log(x) - 1.0/(2*x) + math.lgamma(x) - math.lgamma(n+x)
def pdfln_prime(x):
return (k-3/2)/x + 1.0/(2*x*x) + digamma(x) - digamma(n+x)
while 1:
try:
return ars(pdfln , pdfln_prime , [0.1 , 1.0 , 3.0])
except:
continue
def gen_beta(s , w):
def pdfln(x):
k = len(s)
return -3*math.log(x)/2 - 1.0/(2*x) - k*math.lgamma(x/2)+(k*x/2)*math.log(x*w/2) + np.sum((x/2-1)*np.log(s) - (x*w/2)*s)
def pdfln_prime(x):
k = len(s)
return -(k/2)*digamma(x/2) + 1.0/(2*x*x) + (k/2)*math.log(x/2) + (k*x-3)/(2*x) + np.sum(np.log(s) - w*s)/2 + k*math.log(w)/2
while 1:
try:
return ars(pdfln , pdfln_prime , [0.1 , 1.0 , 3.0])
except:
continue
def rnd_gamma(alpha,theta):
return np.random.gamma(alpha/2 , 2*theta/alpha)
def pdf_normal(x , mean , precision):
return math.sqrt(precision/(2*np.pi))*np.exp(-precision*(x-mean)**2/2)
def pdf_gamma(x , alpha,theta):
lprob = (alpha/2-1)*np.log(x) - alpha*x/(2*theta) - math.lgamma(alpha/2) - (alpha/2)*(math.log(2*theta/alpha))
return math.exp(lprob)
#-- normal gamma distribution pdf
#-- for multivariate case, this should be normal wishart distribution pdf
def pdf_G0(mu , s , lam , r , beta , w):
return pdf_normal(mu ,lam , r)*pdf_gamma(s,beta,1.0/w)
def pdf_lewens(k,n,alpha,ns):
return (k*np.log(alpha)) + sum([math.lgamma(n) for n in ns]) - (math.lgamma(n+alpha) - math.lgamma(alpha))
class DPGMM:
def __init__(self , alpha , lam , r , beta, w):
self.alpha = alpha
self.lam = lam
self.r = r
self.beta = beta
self.w = w
self.cnts = []
self.params = [] #-- component parameters
self.indicators = [] #-- for debug
def __iter__(self):
return self
def __next__(self):
N = sum(self.cnts)
k = len(self.cnts)
probs = [n/(N+self.alpha) for n in self.cnts]
probs.append( (self.alpha)/(N+self.alpha) )
c = np.random.choice( list(range(k+1)) , p=probs)
self.indicators.append(c)
if c<k:
mu,s = self.params[c]
self.cnts[c] += 1
else:
mu = np.random.normal(self.lam , math.sqrt(1.0/self.r))
s = rnd_gamma(self.beta , 1.0/self.w)
self.params.append( (mu,s) )
self.cnts.append( 1 )
return np.random.normal(mu , math.sqrt(1.0/s) )
def likelihood(self,xs):
k_rep = len(self.cnts)
N = len(self.indicators)
Ptmp = pdf_lewens(k_rep,N, self.alpha, np.array(self.cnts))
for j,(m,s) in enumerate(self.params):
mask = (np.array(self.indicators)==j).astype(int)
Ptmp += np.sum( mask*np.log(pdf_normal(xs,m,s)) )
Ptmp += np.log( pdf_G0(m,s,self.lam , self.r, self.beta, self.w) )
return Ptmp
def iGMM(xs , iter=200):
assert(type(xs)==np.ndarray)
N = len(xs)
mean_all = np.sum(xs)/N
var_all = np.sum(xs**2)/N - mean_all*mean_all
stdev_all = math.sqrt(var_all)
best_indicators = np.zeros(len(xs) , dtype=np.int)
best_alpha = 1.0/rnd_gamma(1.0 , 1.0)
best_beta = 1.0/rnd_gamma(1.0 , 1.0)
best_lambda = np.random.normal(mean_all , stdev_all)
best_r = rnd_gamma(1 , 1.0/var_all)
best_w = rnd_gamma(1 , var_all)
best_mparams = [mean_all]
best_sparams = [1.0/var_all]
k_rep = 1
Ptmp = 0.0
for j,(m,s) in enumerate(zip(best_mparams,best_sparams)):
mask = (best_indicators==j).astype(int)
Ptmp += np.sum( mask*np.log(pdf_normal(xs,m,s)) )
Ptmp += np.log( pdf_G0(m,s,best_lambda,best_r,best_beta,best_w) )
Pmax = Ptmp
#print(Pmax)
for it in range(iter):
#-- step(1) update class indicators
h = N
k_rep = np.max( best_indicators ) + 1
assert(h >k_rep)
new_indicators = np.copy( best_indicators )
top_mparams = np.zeros(h)
top_sparams = np.zeros(h)
top_mparams[:k_rep] = np.copy(best_mparams)
top_sparams[:k_rep] = np.copy(best_sparams)
top_mparams[k_rep:] = [np.random.normal(best_lambda , math.sqrt(1.0/best_r)) for _ in range(0,h-k_rep)]
top_sparams[k_rep:] = [rnd_gamma(best_beta , 1.0/best_w) for _ in range(0,h-k_rep)]
for (n,x) in enumerate(xs):
top_nums = np.bincount( new_indicators )
k_rep = len(top_nums)
c = new_indicators[n]
top_nums[c] -= 1
Fvals = np.array([pdf_normal(x , mu , s) for (mu,s) in zip(top_mparams,top_sparams)][:h])
if top_nums[c]==0:
k_rep -= 1
weights = np.ones(h)*(best_alpha/(h-k_rep))
weights[:k_rep+1] = top_nums
weights[c] = (best_alpha/(h-k_rep))
else:
weights = np.ones(h)*(best_alpha/(h-k_rep))
weights[:k_rep] = top_nums
Fvals *= weights
Ftot = np.sum(Fvals)
new_c = np.random.choice(list(range(h)) , p=Fvals/Ftot)
if top_nums[c]==0:
if (new_c==c or new_c>=k_rep):
new_indicators[n] = c
else:
new_indicators[n] = new_c
mask = (new_indicators>=c).astype(int)
new_indicators -= mask
top_mparams[c:h-1] = top_mparams[c+1:h]
top_sparams[c:h-1] = top_sparams[c+1:h]
top_mparams[h-1] = 0.0
top_sparams[h-1] = 0.0
h -= 1
elif new_c>=k_rep:
#-- swap new_c and k_rep
new_indicators[n] = k_rep
top_mparams[k_rep],top_mparams[new_c] = top_mparams[new_c] , top_mparams[k_rep]
top_sparams[k_rep],top_sparams[new_c] = top_sparams[new_c] , top_sparams[k_rep]
else:
new_indicators[n] = new_c
#-- step(2) update component parameters
k_rep = np.max( new_indicators ) + 1
new_mparams = []
new_sparams = []
for j in range(k_rep):
mask = (new_indicators==j).astype(int)
nj = np.sum(mask)
ybar_j = np.sum(mask*xs)/nj
sj = top_sparams[j]
mu_j = np.random.normal((ybar_j*nj*sj*best_lambda*best_r)/(nj*sj+best_r) , math.sqrt(1.0/(nj*sj+best_r)))
var_j = np.sum( mask*(xs-np.ones(N)*mu_j)**2 )
sj = rnd_gamma(best_beta+nj , (best_beta+nj)/(best_w*best_beta+var_j))
new_mparams.append( mu_j )
new_sparams.append( sj )
#-- step(3) update hyperparameters
new_lambda = np.random.normal( (mean_all/var_all+best_r*sum(new_mparams))/(1.0/var_all+k_rep*best_r) ,
math.sqrt(1.0/(var_all+k_rep*best_r)) )
new_r = rnd_gamma(k_rep+1 , (k_rep+1)/(var_all + sum([(x-new_lambda)**2 for x in new_mparams])))
new_w = rnd_gamma(k_rep*best_beta+1 , (k_rep*best_beta)/(1.0/var_all + best_beta*sum(new_sparams)))
new_beta = gen_beta(np.array(new_sparams) , new_w)
new_alpha = gen_alpha(k_rep , N)
#-- compute P
Ptmp = pdf_lewens(k_rep,N,new_alpha,np.bincount(new_indicators))
for j,(m,s) in enumerate(zip(new_mparams,new_sparams)):
mask = (new_indicators==j).astype(int)
Ptmp += np.sum( mask*np.log(pdf_normal(xs,m,s)) )
Ptmp += np.log( pdf_G0(m,s,new_lambda,new_r,new_beta,new_w) )
P = Ptmp
#print(P , k_rep)
if P>Pmax:
best_indicators = new_indicators
best_mparams = new_mparams
best_sparams = new_sparams
best_alpha = new_alpha
best_beta = new_beta
best_r = new_r
best_w = new_w
best_lambda = new_lambda
Pmax = P
else:
pass
return (best_indicators , best_mparams , best_sparams , best_alpha,Pmax)
#-- simple test
if __name__=="__main__":
m = DPGMM(1.1,0.0 , 15.0 , 1.1,0.5)
dat=[]
for _ in range(200):
dat.append(next(m))
dat=np.array(dat)
print( iGMM(dat) )
今後の目標としては、
・infinite HMMを理解する。階層Dirichlet過程というものを使うらしい
・infinite PCFGを理解する
The Infinite PCFG using Hierarchical Dirichlet Processes
https://cs.stanford.edu/~pliang/papers/hdppcfg-emnlp2007.pdf
あたり。Dirichlet過程を使っているものは、須く、Pitman-Yor過程に一般化できるのだと思うので、必要に応じて理解したい。
やりたいこととしては、CCG(Combinatory Categorial Grammar)の確率版であるPCCGを作れるので、"infinite PCCG"ができないだろうかというあたり。CCGはCFGより強いが、どっちもCYK parsingでparseできるし、infinite PCFGができるなら、infinite PCCGができない道理はなさそうな気がする。まぁ、未知の言語の文法構造を、人出によるタグ付けとかなしで、推定できるようになれば嬉しいのだけど、infinite PCFGが流行ってるようにも見えないので、どうなのか(そもそも、文法推論業界で得られてる知見によれば、正規表現程度でも、"ちゃんと"学習するには、正例のみだけではダメで、負例が必要ということなので、もっと強いCFGやCCGでは、尚更、負例の提示が必要そうだけど、今のinfinite PCFGとかは多分正例しか食べ(れ)ないので、そのへんどうなっているか理解したい)
俺たちは何回コインを投げるべきなのか
表と裏が同じ確率で出ることが期待されているコインがあり、コインにイカサマがないことを確認したいとする。なんやかんやあって、統計的手法で、不正の有無(つまり、コインの表が出る確率が50%か否か)に決着をつけることになった場合、コインを何回投げれば十分といえるか
統計的手法にも色々あるかもしれないけど、不正の有無を統計的にチェックしろと言われた場合、オーソドックスに使われるのは仮説検定と思う。表が出る回数は二項分布で表されるが、コインを投げる回数Nが十分大きい時、これは正規分布で近似できるようになって、実際に表が出た回数をx、有意水準は、例えば95%とすると
であれば、コインにイカサマがないという仮説は棄却される(但し、p=1/2で、不正がない時、表が出る確率)。二項分布を正規分布で近似した時の分散σは
で、おおまかに、イカサマがない時の期待値Npから2σ以上外れていればアウトという判定基準になっている(偏差値でいうと、70以上か30以下)。イカサマを疑われている側からすると、イカサマしてない場合でも、冤罪をかけられる確率が5%あるということなので、もっと厳しくしたいと思うかもしれないが、そこらへんは話し合いで同意が取れたとする。もし、表が出る確率が本当に寸分の狂いもなく1/2だったら、(Nが十分大きい時には)コインを投げる回数が100回だろうが5000兆回だろうが、冤罪発生率が5%であることは変わらない。どんなに公正なコインでも、たまたま偏った出目になる可能性は0でないので、残念ながら、これを0%にすることはできない
逆に、本当にイカサマがあって、表が出る確率が60%くらいあった場合、コインをN=5回しか投げないで、有意差がないという結論になっても、明らかに、Nが小さすぎるので、もっとNを増やす必要がある。こういうのは検出力(statistical power)が足りないと言われる。また、極端な話、表が出る確率が90%であれば、少し試してみれば不正があることを見抜けそうであるけど、表が出る確率が60%などの場合は、それより多くのNを必要とすると予想される。確率的なものだから、表が出る確率が0%でない限り、たまたま全部表が出て、イカサマがあるという結論になってしまう可能性はある。けど、例えば、N=10として、表が出る確率が90%あるコインを投げた時、全部表になる確率は35%くらいあるのに対して、表が出る確率が60%のコインを投げて全部表になる確率は0.6%くらいしかない。
数学ばかり見てないで、現実を見ると、コインは表裏対称でないことが多いし、そもそも、人間が肉眼で区別できないほど表裏対称だったら、コインの表が出たのか裏が出たのか判定できない。表面にマーキングしたり汚れがある場合など、人間にとっては微細に思える差であっても、表が出る確率が、1/2から僅かにずれるには十分と思われる。何にせよ、コインが表が出る確率が完全に1/2であることは多分期待できない。仮説検定は、差があるかないかしか調べないので、表が出る確率が50.0001%とかで、わずかに裏より表が出やすいだけという場合でも、Nが大きければ、有意差ありという結果になる可能性は非常に高くなりうる。実際、最初に書いた仮説検定の式を少し変形すると
のようになって、表が出た割合x/Nが許される範囲は、Nが大きくなるに従って狭くなっていく。一方、Nが大きくなれば、x/Nは、大数の法則によって、表が出る真の確率に収束していく。結果として、表が出る確率が寸分の狂いもなく1/2でない限り、有意差ありという結論になる可能性は、Nが大きくなると、どんどん高まっていくことになる。しかし、表が出る確率が、実は50.0001%で、統計的に有意な差があると言われても、常識的な感覚からすると、意味のある差とは思えない
ということを考えると、不正を検出するために必要なNの大きさは、どのくらいの不公平まで許容範囲と考えるかに依存する。つまり、表が出る確率が50±ε(%)の範囲外であればよくないコインであるが、範囲内であれば問題ない水準と考えられるεの大きさを事前に決めておく必要がある。このεの大きさは、有意水準αと同じく客観的な決め方はない。
更に冤罪発生率を0%にできないのと同様、コインに不正がある場合でも、たまたま、表と裏が同じ回数ずつ出てしまう可能性は0にはならない。つまり、不正がある時に不正を100%検出することはできない。とはいえ、Nを大きく取れば、表が出る確率が真に1/2でない限り、このようなことが起きる確率は減っていくと考えられる。というわけで、不正がある時に、正しく不正が検出できる確率を、有意水準と同様に設定する必要がある。これが検出力の定義で、1-βで表されることが多い。というわけで、α、β、εの3つのパラメータを設定すると、不正検出のためにコインを投げる回数Nを決められそうな気がする。
もう少し定量的に評価することを考える。さっきの式に戻ると、有意水準α=0.05と設定して、コインをN回投げた時、表が出た回数xが
を満たせば、コインは公平であると帰結される(前回同様、p=1/2)。
コインに不正があって、表が出る確率がp+εであったとする。このとき、コインが表が出る回数は、やはり二項分布で記述され、上の区間に入る確率を計算することができる。この確率は、表が出る確率がp+εであるにも関わらず、不正はないと結論づけてしまう確率とみなせる。この二項分布は、Nが十分大きい時、平均N(p+ε)、分散N(p+ε)(1-p-ε)の正規分布で記述されることに注意すれば、不正検出ミスの確率は簡単に計算できる。表が出る確率が、例えばp+2εだとか、p+2.4εだというケースもありえるけど、検出力が最も低くなるのは、本当の表が出る確率がp+εか、p-εの場合なので、この2点のみチェックすればいい。更に、p+εの場合と、p-εの場合で、結果は同じことになるので、片方だけ見ればいい。この理由から、以下、εは正としておく
自然言語で書くと状況は複雑だけど、数式で書くと、表が出る真の確率がp+εの時に
という条件を満たすということになる。但し、p=1/2で
とする(積分変数のtはコインを投げて表が出た"割合"を表すので、平均と分散を、それぞれN、N^2で割っていることに注意)。で、積分の変数変換をすると
のようになる。erf関数
を使うと、最終的に
となる。とりあえず左辺を評価するためにεを決める必要がある。まぁ、ε=0.01としておく。
数値的に、いくつか計算してみる。pythonで以下のようなコードを実行すると
import math for N in [10000,15000,20000,25000,30000,35000,40000,50000]: p,eps = 0.5,0.01 u = 1.96/math.sqrt(N/p/(1-p)) sigma2 = (p+eps)*(1-p-eps)/N print(N , 1.0-(math.erf((u-eps)/(math.sqrt(2*sigma2))) - math.erf((-u-eps)/(math.sqrt(2*sigma2))))*0.5)
以下のような結果を得た。
10000 0.5159939775494847 15000 0.6877923080943419 20000 0.8074680951250419 25000 0.8854187382561318 30000 0.9337611517760266 35000 0.9626265172738802 40000 0.9793451544602965 50000 0.9940083977523118
大体
N=20000回で、検出力0.8ちょっと
N=25000回で、検出力88.5くらい
N=50000回で、検出力0.994...
などとなる。つまり、20000回くらいやれば、コインが表が出る確率が51%の時、80%以上の確率で有意差がある(有意水準は5%)という検定結果を得られる計算。まぁ、数万回くらいはコインを投げることになりそう。一日が86400秒なことを考えると、実行するには、くそだるい数字である。
検出力は、0.8とか0.9にすることが多いらしいが、αやε同様、客観的に決める基準はないので慣習的なものに過ぎない。また、有意水準を厳しくすると、同じ検出力を得るには、より多くの試行回数を必要とする。例えば、3σを要求すると
10000 0.15860713527390746 15000 0.29094698846888134 20000 0.43187317605341713 25000 0.5644691796250857 30000 0.6787457862331097 35000 0.7708974857449786 40000 0.841393149895479 50000 0.9295476650393496
のようになった。εが厳しすぎる気もするので、5%くらいまで緩めると、数千回程度でも十分な検出力を期待できるようになる。大雑把に言って、有意区間の幅は、Nの平方根に反比例するので、(αやβを変えずに)εを1/10にしたければ、Nを100倍にする必要があり、逆にεを10倍にしてよければ、Nは1/100程度まで小さくできる。このおかげで、今の場合、(例えば、濡れ衣を着せてやろうと思って)Nを大きくして小さな差でも有意になるようにしようとしても、実際に実行するのはかなり難しくなる。例えば、εを0.001%にして(上で計算したεが1%の時と)同じ有意水準、検出力で検定しようと思えば、100万倍ものデータ量が必要となる。
元ネタ(検定力と検出力は同じもの):
【翻訳】ダメな統計学 (3) 検定力と検定力の足りない統計
http://id.fnshr.info/2014/12/17/stats-done-wrong-03/
これは本にもなってる。Webページには、コインの例が載っているのだけど、上に書いたような計算は載っていなかった(数式を載せるような本ではないせいだろうと思う)ので、計算してみた。適切な統計学の教科書には例題として載っていそうな気もするけど、私が大学で受けた授業では、検出力/検定力とか聞かなかった気がする。
演習問題:各面が同じ確率で出ることが期待されている(6面)サイコロに不正がないかを統計的に調べたいとすると、俺たちは一体、何回サイコロを投げるべきか？
方針:教科書的には、カイ二乗(適合度)検定を行うことが第一歩となるだろうと思う(他の検定手法もあるかもしれないが)。
次に、コインの場合のεに相当する量を何らかの形で定義する必要がある。これに正しい方法というのがあるわけではないけど、サイコロの各面の期待される出現確率をp_i(i=1,...6)として、実際の出現確率をp_i'とする時、次のような量を定義すると便利
各面が同じ確率出ることを期待されているという条件は、勿論
を意味する。wは効果量(effect size)と呼ばれる。例えば、1〜3が出る確率が1/6+1/100で、4〜6が出る確率は1/6-1/100である(各面1%くらいの誤差)とすると、wは0.06となる。実際の研究では、wは、0.1,0.3,0.5などの値が選ばれることが多いらしい
効果量の定義の仕方は、検定の手法ごとに異なりうるし、検定の手法を決めても、必ずしも標準的なものが一つあるというわけでもないらしいので、Cohen's wという名前で呼ばれることもあるらしい
Cohen's w
https://en.wikipedia.org/wiki/Effect_size#Cohen.27s_w
それで、サイコロを投げる回数をNとして、wを使うと、カイ二乗値は、"平均的には"
となり、Nに比例して大きくなることが分かる(w=0でも、期待値は自由度kに等しくなるので、本当の意味での平均ではない)。カイ二乗値は、帰無仮説が正しい時(つまり、p_i=p_i'が成立する時)には、近似的にカイ二乗分布に従う(今の場合、自由度は5)というのが、カイ二乗検定の要旨だったけど、p_iとp_i'がずれている時には、非心カイ二乗分布で近似されることは直感的に分かる。検出力は、この非心カイ二乗分布を積分して計算できるけど、非心度λを知る必要がある。ところで、非心カイ二乗分布の期待値が自由度kとλの和で書けることに注意すれば、λと、上で計算したカイ二乗値の値Nw^2とが関係することは意外ではない。etc. 具体的な計算は面倒なので省略
■
球面上の測地流は、古典力学系として見ると、球面の余接空間を相空間としている。その関数環はPoisson代数であり、Poisson括弧の具体的な計算も、Diracの方法によってDirac括弧として計算できる。このPoisson代数を眺めると、Euclidean Lie algebraの対称代数からの全射準同型があることに気づく。というわけで、量子化した代数は、Euclidean Lie algebraの普遍展開環の何らかの商であると考えるのが自然。このことに最初に気付いた人が誰なのかは知らないけど、
Fundamental algebra for quantum mechanics on S^D and gauge potentials
http://aip.scitation.org/doi/abs/10.1063/1.530099
には、そういう風に書いてある。
Euclidean Lie algebra $e(N+1)$の既約ユニタリ表現は、Wignerの方法で群の表現から作ることができて、への$Spin(N+1)$作用への軌道は、一般に球面で、little groupは$Spin(N)$となる(原点の軌道は一点からなりlittele groupはSpin(N+1)となるので、この場合は誘導表現はSpin(N+1)の既約表現である)。なので、$e(N+1)$の無限次元既約ユニタリ表現は、球面の半径$r$と$Spin(N)$の既約ユニタリ表現(WeylのユニタリトリックによってLie環$so(N,C)$の有限次元既約表現と言っても同じ)でラベルされる。Spin(N)の自明表現に対応するのは、標準的なL^2(S^N)への作用である。
$e(3)$の場合は、自明でないcoadjoint orbitは全てと微分同相になる(
とすると、pをp-(b/a)xに移すことで、O(a,0)とO(a,b)の微分同相写像が得られる)。$e(3)^{*}$や
上で定義される力学系は色々あると思うけど、(Liouvilleの意味での)可積分系としては、球面上の測地流、spherical pendulum、Lagrange top、Kowalesvki topなどがある。後者2つのintegrable topは$e(3)^{*}$上で定義されているが、coadjoint orbitに制限することで、
上で定義されていると思うこともできる($e(3)^{*}$上のHamilton力学系は非正準ハミルトン力学系であり、coadjoint orbitはいわゆるCasimir不変量で特徴付けられるsymplectic leavesとなっている)。
参考)3次元Euclid代数の無限次元既約ユニタリ表現については、以下に要点をまとめた
http://formalgroup.tumblr.com/post/160390024555/3%E6%AC%A1%E5%85%83euclid%E4%BB%A3%E6%95%B0%E3%81%AE%E6%97%A2%E7%B4%84%E3%83%A6%E3%83%8B%E3%82%BF%E3%83%AA%E8%A1%A8%E7%8F%BE
$e(3)^{*}$上で定義されたHamilton力学系を量子化することを考えると、自然な波動関数の空間/状態空間は、$e(3)$の既約ユニタリー表現である。これは、上にも書いた通り、一つではなく、無数にある。物理の実際の問題では、波動関数の空間はどれか一つに決める必要があるが、今の場合、$e(3)$のユニタリ表現であり、それを既約分解すると(※)、既約ユニタリ表現を得る。どのような既約ユニタリ表現が出るかは実際の物理の問題設定に依存する。
(※)I型の局所コンパクト群の場合、任意のユニタリ表現は、直積分の形に一意に分解できることが知られている(らしい)
#Casimir作用素は、U(e(3))の中心の元を指す。一方、$e(3)^{*}$上の関数環は対称代数S(e(3))=gr(U(e(3)))であり、grによって、Casimir作用素をS(e(3))に持って行った像がCasimir不変量となる。一般には、Poisson代数の元の内、任意の要素とPoisson可換なものをCasimir不変量と呼ぶのだけども、Poisson代数がLie環の対称代数の場合は、Casimir不変量は、このようにして得られる。
数学としては、各既約ユニタリ表現ごとにHamiltonianのスペクトルを計算できれば、問題は解けたことになる。Kepler系や球面上の測地流の時は分岐則から答えがわかるのでほぼ自明な感じだったが、Kowalevski topやLagrange topはそこまで簡単ではない。こうした量子可積分系の解法が理解されてきたのは1980年代以降のことだと思う。例えば、Kowalevski topのLax行列は1980年代後半に複数の人によって発見されている。
Lax representation with a spectral parameter for the Kowalewski top and its generalizations
http://link.springer.com/article/10.1007/BF00403470
The Kowalewski Top 99 Years Later: A Lax Pair, Generalizations and Explicit Solutions
https://projecteuclid.org/euclid.cmp/1104178400
不思議なことに、Kowalevski topのLax行列は、so(3,2)(のループ代数)に値を取る(Lagrange topの場合は、so(3,1)で、一般化して、so(p,q)-topというものが得られる)。この2つのintegrable topはGASDYM(generalized anti self-dual Yang-Mills)方程式の次元簡約によっても得られるが、その時ゲージ群は、Kowalevski topの場合はSO(3,2)、Lagrange topの場合はSO(3,1)とする。このへんの理由は、私にはよく分からない。
やや風変わりなintegrable topとして、Goryachev-Chaplygin(GC) topがある。これは、$e(3)^{*}$全体では可積分系ではないが、これが可積分になるようなcoadjoint orbitが存在する(大体文献見ると、$r=1$に固定しているのだけど、$r$の値は0でなければ何でもよくて"スピン"が0になることのみが可積分性の条件の気がする。半径の値はスケーリングで規格化してしまえば、何でも同じということかもしれないけど、一応同値でない既約表現が出るので、気にはなる)。e(3)のスピン0の表現空間(完備化すると、L^2(S^2)になる)には、so(3,2)の極小表現あるいは、同じことであるがsp(4,R)の振動子表現(の片割れ)が実現できる。この事実は、quantum GC topを"解く"のに利用された(最終的に、2つの無限次三重対角行列の固有値を計算する問題に帰着する)
Exact solution of the Goryachev-Chaplygin problem in quantum mechanics
http://iopscience.iop.org/article/10.1088/0305-4470/15/6/015/meta
など。ただ、この論文の計算は、Lie環の元の平方根を取るなど、やや危うい操作がある
Goryachev-Chaplygin top and the inverse scattering method
http://link.springer.com/article/10.1007/BF02107243
では、この問題は回避されていて、ついでにL-operatorを与えるなど、系統的な感じにもなっているが、最終的に出る式は同じ。ところで、後者の論文のR-matrixの形とRLL関係式を見ると、Yangian Y(sl(2))との関係を想像させる(当時はまだYangianは知られてなかったはず)。detailをcheckしてないので正しいのかどうか保証できないけど、以下の論文は、まさにその問題を考えているように思う
Transfer Matrix and Yangian Related to Chaplygin-Goryachev Top
http://iopscience.iop.org/article/10.1088/0256-307X/13/9/001
Yangian, Truncated Yangian and Quantum Integrable Systems
https://arxiv.org/abs/cond-mat/9509081
上の論文では引用されてないけど、truncated Yangianは、1988年にCherednikが考えたものらしい(当時は"truncated"という言い方はされなかった模様)。T(u)を有限次の展開で打ち切るので"truncated"と言われているが、打ち切る次数の大きさをレベルと呼ぶ。その言い方に従えば、上で考えているのは、レベル3のtruncated Yangianということになると思う
また、2000年頃になって、ある種の有限W代数が、truncated Yangianと同型になるということが発見された(数学的な有限W代数の定義では、gl(Np)のサイズpのJordan blockがN個並んでるような冪零元eを取ってきて、W(gl(Np),e)を考えると、これがY(gl(N))のレベルpのtruncated Yangianと同型になる、ということらしい(多分...))
Yangian realisations from finite W algebras
https://arxiv.org/abs/hep-th/9803243
RTT presentation of finite W-algebras
https://arxiv.org/abs/math/0005111
Yangians and W-algebras
https://arxiv.org/abs/1305.4100
#Higgs algebraは、(gl(N)ではなく)sl(N)の場合であるが、上のタイプの有限W代数であり、Y(sl(2))の(レベル2の)truncated Yangianとして実現されるのだと思う(要確認)
#現在、様々なタイプの"量子群"が知られており、それらの関係を理解するのに、RLL/RTT関係式を使うのは見やすいように思う。"sl(2)シリーズ"の場合のみではあるけど、以下の論文は各種"量子群"の関係が書いてあるので、沢山あるなぁと思って、ゲンナリするのに丁度いい
Cladistics of Double Yangians and Elliptic Algebras
https://arxiv.org/abs/math/9906189
#上の論文に出てくるR-matrixは、全て差法性R(u,v)=R(u-v)を持つが、応用上は、そうでないR-matrixが出てくることもある(ShastryのR-matrixなど)。こうしたR-matrixでも、RLL関係式を通して、何らかの代数構造が得られるらしいけど、謎が多い
その後、一般化として、gl(N)の場合、一般の有限W代数がtruncted shifted Yangianというもので実現されるという話があり、
Shifted Yangians and finite W-algebras
https://arxiv.org/abs/math/0407012
Representations of shifted Yangians and finite W-algebras
https://arxiv.org/abs/math/0508003
有限W代数やtruncated Yangianの表現論は、まだわからないことも多い(私が理解してないことが多いというわけだけでもなく)ので、現状、GC topを解くのに直接役立つわけではないし、今後も役立つのかどうかわからない。まぁでも、動機として、可積分系のLaxペアは、ヒマな人が試行錯誤して勘で見つけてるようにしか私には見えないものがあるので理解できる形にしたいというのもあるし、これらの代数から新しい可積分系を得ることもあるかもしれないし
Kowalevski topやLagrange topの場合は、classical r-matrixの形は分かっている
Classical R-matrices for generalized so(p,q) tops
https://arxiv.org/abs/nlin/0401025
これらの系に於いて、古典論のレベルでは、通常のループ代数ではなく、twisted loop algebraが出てくるので、量子化すると、多分twsited Yangian(のquantum double?)になるのでないかと思っているけど、定かではない(つまり、今までの話を素朴に総合するならtruncated twisted Yangianのようなものを考えるべきなのでないかと思うけど....)。この手の対称対から得られる力学系としては、他にEuler-Manakov topなどがある
#Yangianは、(1)XXX模型や(一次元)Hubbard模型のような格子模型、(2)(massiveな)integrable field theoryなどの対称性として出てきて、色んな定義があるけど、有限自由度の古典可積分系との関係では、half-loop algebraの量子化として導入される(Lie bialgebraの量子化に関する一般論から出るらしい)。古典可積分系で、ループ代数が出るというのは、多分1970年代から知られていた(ループ変数は、Laxペアのスペクトルパラメータ)。有限自由度の可積分系でYangianのような巨大な代数がそのまま対称性として働くわけではなく、truncationが効いてくるのだろうと私は理解している
・half-loop algebra <-> Yangian
・loop-algebra <-> quantum double of Yangian
・twisted half-loop algebra <-> twisted Yangian
・twisted loop algebra <-> ???
#twisted YangianはOlshanskiiが考えたのが始まりで、その後、Mackayらによって一般のsymmetric pairに対するtwisted Yangianが定義された模様
Twisted Yangians and symmetric pairs
https://arxiv.org/abs/math/0305285
Twisted Yangians for symmetric pairs of types B, C, D
https://arxiv.org/abs/1407.5247
Representations of twisted Yangians of types B, C, D: I
https://arxiv.org/abs/1605.06733
