■
球面上の測地流は、古典力学系として見ると、球面の余接空間を相空間としている。その関数環は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
俺たちは何回コインを投げるべきなのか
表と裏が同じ確率で出ることが期待されているコインがあり、コインにイカサマがないことを確認したいとする。なんやかんやあって、統計的手法で、不正の有無(つまり、コインの表が出る確率が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. 具体的な計算は面倒なので省略
ノンパラメトリックベイズ(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とかは多分正例しか食べ(れ)ないので、そのへんどうなっているか理解したい)
