nonparametric
Folders and files
| Name | Name | Last commit date | ||
|---|---|---|---|---|
parent directory.. | ||||
nonparameteric: ノンパラメトリックベイズ
=======================================
statistics.lisp
ノンパラベイズ用の統計＆行列ユーリティティ
（主に行列計算のために）hjs.util.*を使っている
todo: 確率系関数をstatistics/statistics.lispとマージ
重複しているいくつかの乱数はこちらの方が高速にスケールするアルゴリズムとなっている
=======================================
dpm.lisp
Dirichlet Process Mixture(DPM): クラスタ数(k)を推定しながら与えられた事前分布から生成されたクラスタへのボトムアップクラスタリングを行う
Chinese Restaurant Process(CRP)を用いる他のノンパラベイズ手法にライブラリとして提供される
todo: cluster-validationを使えるように結果を整形する？
class:
- cluster
-- gaussian-cluster
クラスタパラメタやそのサンプルに必要な計算の途中経過を保持する目的のクラス
- dpm
-- gaussian-cluster
data: pointのvector
- dp-distribution
-- dp-gaussian
DPM用の分布定義クラス
struct:
- point
任意のデータ型とそのクラスタ割り当ての組
methods:
- (initialize dpm)
初期解をサンプルする
戻り値は更新されたhyper-parameter
- (sampling dpm)
sliced Gibbs samplingによって解を更新する
戻り値は更新されたhyper-parameter
- (make-cluster-result dpm)
クラスタリング結果を入力データのベクトルで返す
Accessor:
- (dpm-clusters dpm)
現在のクラスタリング結果を返す
現状では割り当てゼロのクラスタも見えてしまう
- (dpm-k dpm)
割り当てのあるクラスタ数の合計kを返す
- (dpm-hyper dpm)
hyper-parameter
- (dpm-base dpm)
DPMが仮定している「事前分布」
- (dpm-data dpm)
クラスタリング対象 pointのベクトル
- (estimate-base? dpm)
DPMが事前分布のパラメタをサンプルから推測するべきか(defaultはnil)
tにセットすると学習がうまくいかなくなる恐れもある
- (cluster-size cluster)
- (point-data point)
- (point-cluster point)
読んで字の如し
自分で事前分布を定義するには、独自のdp-distributionのサブクラス、必要であればclusterのサブクラスを定義した上で、その元で
- density-to-cluster
- add-to-cluster
- remove-from-cluster
- make-new-cluster
- base-distribution
- sample-cluster-parameters
- sample-distribution (optional)
などのメソッドを定義する必要がある。
作った事前分布のinstanceを(dpm-base dpm)にセットするか、defaultがそうであるようなdpmのサブクラスを定義して完成
- density-to-cluster (cluster data &rest other-data &key &allow-other-keys)
clusterにデータが所属する確率を返すメソッド
other-dataは条件付き確率を計算したい場合に使う
後述のbase-distributionと合わせて、スケールがあっているなら必ずしも和が1である必要は無い
（むしろ直截的に事後確率を用いれば和は1にはならないが、これを利用して計算をサボれる場合がある）
todo: 引数にdistributionを取る必要がありそうか？
- add-to-cluster (cluster data &rest other-data &key &allow-other-keys)
clusterにデータを追加するメソッド
dataはpoint構造体ではなく、生のデータ(point-data data)である(下記すべて同じ)
defaultメソッドはsizeをインクリメントするだけ
クラスタ重心を求めるために全データの合計を計算するなどのステップはここに入れてしまうとよい
- remove-from-cluster (cluster data &rest other-data &key &allow-other-keys)
clusterからデータを削除するメソッド
defaultメソッドはsizeをデクリメントするだけ
通常はadd-to-clusterと対称に定義するべき
- make-new-cluster (distribution data &optional discarded-cluster)
事前分布と新規客の元で、新しいクラスタをサンプルするメソッド
discarded-clusterは着席ゼロの空クラスタが渡されるか、もしくはnilである（ゴミ抑制）
- base-distribution (dpm distribution data &rest args &key &allow-other-keys)
事前分布とあるデータに対して、既存のクラスタではなくて新しいクラスタに割り当てられる確率を返すメソッド
\int f(x|\theta)G(\theta)d\theta に対するx=dataの値
- sample-cluster-parameters (cluster distribution dpm)
現在のクラスタ割り当てと事前分布からクラスタに存在するパラメタを再サンプリングするメソッド
クラスタ割り当て結果が必要な場合にはclusterに保持するとよい
データ点全体の分布が必要な可能性があるので、引数にdpmを取るが、単純な例では必要ない
- sample-distribution (dpm distribution)
現在のデータ点、およびクラスタパラメタの分布から事前分布のパラメタをサンプルするメソッド
事前分布に関する情報が不完全な場合には、得られたMCMCサンプルから事前分布のパラメタもサンプルする方がよい場合もある
しかし基本的には人間が事前知識を入力する方が望ましい
estimate-base?をtにセットしないなら定義する必要はない
========================================
multi-dpm.lisp
多変量混合正規分布に対するDPM
DPMに対する事前分布の定義の仕方の例
実行例：
DPM(128): (setf data (make-multivar-gauss-dataset 3 4 200)) ;; 人工データ生成用関数
distributions ;; 作成した分布のリストを印字する
(#(3.0227228181608674 -4.917328883250097 9.416279840427581) ;; 平均ベクトル
#2A((0.6149354366098384 0.0 0.0) ;; 分散
(-0.4055310942489165 1.1814636317087717 0.0)
(0.39123857033292375 0.656671363189244 1.6453273683222376))
61) ;; データの数
(#(4.725178282295572 -10.80592572131948 10.935760522066051)
#2A((1.6794388919977667 0.0 0.0)
(0.10077996946225677 2.8513112342592026 0.0)
(0.01318929234028677 -1.5334430908070291 3.911133356965887))
55)
(#(-12.750888916648064 2.9784339490902108 2.2432568569974394)
#2A((1.8236861940336782 0.0 0.0)
(1.26978464750881 2.074938458977714 0.0)
(-0.5586229203010648 0.2864853832500676 0.513250979823732))
43)
(#(3.706127055835794 -2.100318468340544 -5.936122188279907)
#2A((1.1387192901988026 0.0 0.0)
(-0.8338853318826855 1.4984944510243 -1.1102230246251565e-16)
(0.15422119671257709 1.755283189530601 0.6445149218041165))
41)
;; 上の分布から得られたランダムなデータ列
#(#S(POINT :DATA #(3.774485958242961 -10.237199724529862 10.79342433539454) :CLUSTER NIL)
#S(POINT :DATA #(7.969900415944767 2.0638137822323888 -8.574501407746817) :CLUSTER NIL)
#S(POINT :DATA #(4.138140190077053 -11.044781799519692 10.894942699510448) :CLUSTER NIL)
#S(POINT :DATA #(-11.390991006937876 2.352468388156428 4.752483799387133) :CLUSTER NIL)
#S(POINT :DATA #(1.0505170887995634 -4.895694882479956 -3.5828556362300024) :CLUSTER NIL)
#S(POINT :DATA #(4.793301447804544 -4.241785499366994 10.085430501196928) :CLUSTER NIL)
#S(POINT :DATA #(4.492124617269321 -10.835153936892327 10.690929797020202) :CLUSTER NIL)
#S(POINT :DATA #(1.7952202606807408 -4.144550190592244 9.615672379598411) :CLUSTER NIL)
#S(POINT :DATA #(3.380615067547286 -2.638131867036481 -6.454694304762895) :CLUSTER NIL)
#S(POINT :DATA #(-0.13180219410625105 -5.074561211442568 10.582243664171607) :CLUSTER NIL) ...)
DPM(179): (setf m (make-instance 'multivar-gauss-dpm :data data :dim 3))
#<MULTIVAR-GAUSS-DPM @ #x1003a25912>
DPM(180): (initialize m)
13.303837839977257
DPM(181): (time (bench m)) ;; 50回のサンプリング後に現在の状態を印字するループを5回回す
; While compiling (:INTERNAL (:ANONYMOUS-LAMBDA 87) 0):
Warning: Free reference to undeclared variable M assumed special.
---------- iteration 1 ----------
0:
#<size: 58
ave: #(2.927175216693459 -4.846338033224291 9.372464414573544)
std: #2A((0.6976461313410682 0.0 0.0)
(-0.08930776419439851 1.0480008891842754 0.0)
(0.588940897301858 0.5899858608031886 1.4216948269886558))>
1:
#<size: 55
ave: #(4.732399640682938 -10.873523996372105 10.914769284974337)
std: #2A((1.547503724624603 0.0 0.0)
(-0.1097597986394086 3.322611243972814 0.0)
(0.3406420972373237 -2.6330490450676383 4.11058979065462))>
2:
#<size: 40
ave: #(-12.490480051959969 2.949746546416199 2.6247548519158133)
std: #2A((1.7404103077909476 0.0 0.0)
(1.4259513030215993 2.18495039535426 0.0)
(-0.6225890647922969 0.17772685591186194 0.40685610316602205))>
3:
#<size: 23
ave: #(4.296988793983612 -1.5894232076272998 -6.459805813434724)
std: #2A((1.1208661554054586 0.0 0.0)
(-0.821776357188609 2.014179347709142 0.0)
(0.008628758127304521 2.3890564597279313 0.72050609086323))>
4:
#<size: 6
ave: #(3.4358108912056053 -2.8448623968955693 -5.165353635451156)
std: #2A((2.971596050986194 0.0 0.0)
(2.99008187866468 0.6235136724562264 0.0)
(2.391225089915498 0.675428781062474 0.9421551364116058))>
5:
#<size: 5
ave: #(-0.4513647324971928 -6.391492366210508 -2.292455705810453)
std: #2A((1.43527226343607 0.0 0.0)
(-2.641563630102838 1.6727404120994602 0.0)
(0.33934554407951123 2.1928082108878892 1.0764248574171718))>
6:
#<size: 3
ave: #(1.0034301894984483 -5.085339667836983 -3.330031890325851)
std: #2A((3.2110462030178297 0.0 0.0)
(-2.9163343463569853 2.2515586838086175 0.0)
(0.6918864984753642 3.104504271126779 0.8884490791378827))>
7:
#<size: 2
ave: #(0.41893594239215315 -5.399690497228074 -3.1880772959963566)
std: #2A((5.67278702826423 0.0 0.0)
(-7.6894867320249505 2.7862114540806484 0.0)
(-3.142478203436282 2.5856170048595812 0.690369311980487))>
8:
#<size: 2
ave: #(4.446019307102985 1.143722868803642 -3.7303434596422234)
std: #2A((0.7554112250919517 0.0 0.0)
(0.21641014430376873 0.9894094631123735 0.0)
(-0.2868437790614707 0.7761570881896019 0.23888263504264032))>
9:
#<size: 2
ave: #(1.744335794578009 -7.325473363459821 10.345797113840007)
std: #2A((2.301450914896525 0.0 0.0)
(1.5072182101739973 0.5758809098775008 0.0)
(0.8474592693513199 0.1905266656569395 0.7217452287025308))>
---------- iteration 2 ----------
0:
#<size: 61
ave: #(2.8069676457707793 -4.942451267565749 9.512528548611881)
std: #2A((0.6799038400369993 0.0 0.0)
(-0.40444857239877 1.1679362196076075 0.0)
(0.40236064539042604 0.6552619132187644 1.384538255002976))>
1:
#<size: 55
ave: #(4.627570603871044 -10.758847539443064 10.987401543690158)
std: #2A((1.8976322689144043 0.0 0.0)
(0.5424146707856068 3.7452200275288616 0.0)
(-0.43909417691543023 -1.3136461786448672 4.0712577217836))>
2:
#<size: 40
ave: #(-12.662360167739097 3.140428740406515 2.6591768840908143)
std: #2A((1.7305310638280822 0.0 0.0)
(1.937022653476898 1.8616300697326469 0.0)
(-0.5151964991296756 0.21224076226555388 0.38817659242052754))>
3:
#<size: 28
ave: #(3.653125072232231 -1.9515764056098193 -5.903537689044709)
std: #2A((1.1151617556151736 0.0 0.0)
(-0.6392129447096362 1.7026883413996612 0.0)
(-0.030834646488631073 1.741199258817217 0.6300100695389105))>
4:
#<size: 3
ave: #(1.1606708442642515 -2.1746006081025095 -2.9196733052496864)
std: #2A((1.9989818293306156 0.0 0.0)
(0.8988751053200703 1.3749888422191847 0.0)
(-0.19445615911290143 -0.49249121417460423 0.5496154174487236))>
5:
#<size: 3
ave: #(0.03161005359511471 -3.792702923198943 -2.095354083244233)
std: #2A((1.3426485722996682 0.0 0.0)
(0.5923353404445896 0.838655048763351 0.0)
(1.0635970698812685 0.3964052948822824 0.3791505462151047))>
6:
#<size: 1
ave: #(1.348920186647181 -4.623310375775923 -3.969615297803295)
std: #2A((5.3682663028031286 0.0 0.0)
(-4.754291682242515 4.9774230665772325 0.0)
(0.9929242871292591 7.264166195534891 0.5689067096383374))>
7:
#<size: 1
ave: #(3.9376321931818024 -0.011855704787635535 -3.3531719635718833)
std: #2A((1.5365838747680391 0.0 0.0)
(-0.4062840229838975 0.898778363829382 0.0)
(-0.37277701671464936 -0.21374268553867856 0.9770554778017979))>
8:
#<size: 1
ave: #(0.015234369193675323 -4.023344212622708 -0.7009045462957213)
std: #2A((1.091490854883612 0.0 0.0)
(0.4368016964350129 1.270459905309591 0.0)
(0.17503019667099795 0.7450348942986575 0.44615356161913916))>
9:
#<size: 1
ave: #(-6.533570745919246 0.29180235730027926 -0.6386960341817608)
std: #2A((1.5418293577759248 0.0 0.0)
(0.5189984692264599 1.0235381517651103 0.0)
(-0.4077448561915826 -0.6284787474836747 1.0672864342039536))>
---------- iteration 3 ----------
0:
#<size: 61
ave: #(2.5097125632892836 -5.291741497832986 9.498809917234665)
std: #2A((0.5411140518760066 0.0 0.0)
(-0.36802848299336105 1.0317747403435378 0.0)
(0.36385618896617833 0.6248554673055129 1.1271977542069336))>
1:
#<size: 55
ave: #(4.763751577143368 -10.719943140826631 10.978814891747739)
std: #2A((2.046111424308298 0.0 0.0)
(0.20043508360623385 3.632567468239127 0.0)
(0.903543862203009 -1.5535051461195748 3.7379276477205003))>
2:
#<size: 43
ave: #(-12.918686348812814 3.0652096087178986 2.0123801176515306)
std: #2A((1.5969821847124224 0.0 0.0)
(1.9435922489689714 2.0843249795196273 0.0)
(-0.479766130113108 0.21063788017107146 0.449176752593068))>
3:
#<size: 34
ave: #(3.077347661016097 -2.606910019735932 -5.624127453466127)
std: #2A((1.04870266591221 0.0 0.0)
(-1.3834672031726198 1.742183003796471 0.0)
(-0.6285767145495247 1.667946748656608 0.633102715722574))>
4:
#<size: 5
ave: #(0.7859632764612294 -5.466241760115667 -3.5690975618424807)
std: #2A((2.280148594034162 0.0 0.0)
(-1.2861622900571033 2.5947223043251757 0.0)
(0.2620340401231253 2.035078886956046 1.6060393252745941))>
5:
#<size: 1
ave: #(4.339814977740541 0.9040311920161834 -3.802637222638452)
std: #2A((1.3215627583281466 0.0 0.0)
(0.455161781442808 0.5595770487678385 0.0)
(-0.3420897484465336 -0.05069328847789344 0.796814007648923))>
6:
#<size: 1
ave: #(3.8987749205860007 0.5066169984636779 -2.808698511732912)
std: #2A((0.8224163188330691 0.0 0.0)
(0.2795938798948154 1.6314200452124068 0.0)
(-0.18571885683120082 0.4886647129861812 0.5656240050271633))>
7:
#<size: 0
ave: #(-5.446072219664096 0.4646337166114336 2.235052000138044)
std: #2A((1.1637049172868104 0.0 0.0)
(0.18990724546266577 1.043639337426633 0.0)
(0.4817629147795495 -0.35455894505533386 0.35026423607301477))>
8:
#<size: 0
ave: #(3.306418211625072 -3.013710665743015 -5.23985889020107)
std: #2A((1.873169594874314 0.0 0.0)
(-1.2382008061127414 2.062795557190719 0.0)
(0.8276950766672347 2.5778583221868097 0.824809776330044))>
9:
#<size: 0
ave: #(2.278517597481957 0.08409290340986808 -2.5791591668763263)
std: #2A((0.71702697483601 0.0 0.0)
(-0.047064641966757234 1.4651513859784986 0.0)
(0.24340368507424803 -0.22294575610795062 0.49565870876549945))>
---------- iteration 4 ----------
0:
#<size: 61
ave: #(2.4828425826092313 -4.92890727148019 9.427873678109767)
std: #2A((0.6097810269637248 0.0 0.0)
(-0.35610449860921695 0.9780816042038336 0.0)
(0.4228667948828795 0.5911619856689561 1.528804150997672))>
1:
#<size: 55
ave: #(4.854148372513157 -10.742942106189865 10.970786309253226)
std: #2A((1.6859348698841365 0.0 0.0)
(-1.342041885373115 3.999346017012868 0.0)
(1.2916162069428707 -2.0651591555442987 4.110448955846518))>
2:
#<size: 42
ave: #(-12.986199989412432 3.1278997054574313 1.915325691312002)
std: #2A((1.764385199340075 0.0 0.0)
(1.9842071227384714 2.436568532152146 0.0)
(-0.4438355481802078 0.23890665670147898 0.38612303726843816))>
3:
#<size: 33
ave: #(4.239500836173111 -1.7523117038055183 -6.329168160595539)
std: #2A((0.9813093554656801 0.0 0.0)
(-1.0928237778558418 1.2544176642142835 0.0)
(-0.5039803152015857 1.3693810995205458 0.6444713986125342))>
4:
#<size: 4
ave: #(0.27418047898082953 -5.57314972250134 -2.9125553045205725)
std: #2A((4.346262536846539 0.0 0.0)
(-1.9468342722390304 3.9830605487034987 0.0)
(1.460254311129598 3.4422406868617217 1.6539971513080065))>
5:
#<size: 3
ave: #(1.9498946121429213 -3.6286175508176015 -4.5552549551748305)
std: #2A((4.083883248527007 0.0 0.0)
(-2.037659804708729 1.8951425217612128 0.0)
(1.9568393874225325 2.7024612670327777 0.645747404795075))>
6:
#<size: 1
ave: #(-0.5619418957685309 -3.834290836510284 -0.2717122069053709)
std: #2A((0.8996071936021539 0.0 0.0)
(-0.5390163706495277 0.4726829786356856 0.0)
(0.11965520465420415 0.7339110130410332 0.8917265397166774))>
7:
#<size: 1
ave: #(-4.215484934002681 2.5229863325600803 1.2119453465860446)
std: #2A((1.5243030335176517 0.0 0.0)
(-0.7243290671359278 0.45622839520097 0.0)
(-0.3153071733951639 0.39968587585734205 1.275269681537959))>
8:
#<size: 0
ave: #(0.2815886503953172 -3.8676202722828004 -0.4978001680451093)
std: #2A((0.8211466872346592 0.0 0.0)
(0.017360251670413003 0.33905705331625 0.0)
(0.3242917921606121 0.4450200925426773 0.08241750694008637))>
9:
#<size: 0
ave: #(4.0652723644354145 0.13662564352278272 -3.9585615306646265)
std: #2A((1.4439147642784866 0.0 0.0)
(0.10777104032621353 0.7764016473714035 0.0)
(-0.4752170187041059 -0.4284109898229721 0.6360579626925299))>
---------- iteration 5 ----------
0:
#<size: 61
ave: #(2.8389494105390347 -4.847297277243273 9.449730517683742)
std: #2A((0.7338414185028997 0.0 0.0)
(-0.22663066295009468 1.0912295961011884 0.0)
(0.4350865305606222 0.7784648620249829 1.1351873763467162))>
1:
#<size: 55
ave: #(4.832141218234674 -10.729679314056742 10.95703758791401)
std: #2A((1.5135098695308806 0.0 0.0)
(-0.3278868323282931 3.3439134501606333 0.0)
(0.3892161360349809 -0.7996195071882126 3.821915713387091))>
2:
#<size: 42
ave: #(-13.018193397100104 3.1377153452573947 1.754553157044611)
std: #2A((1.7320434023168374 0.0 0.0)
(1.817833708039712 1.9728925311897771 0.0)
(-0.5872114464718041 0.1765474573821052 0.5554686339321179))>
3:
#<size: 19
ave: #(3.8869659653433524 -1.8130105563013674 -6.115404170093998)
std: #2A((1.2727935041129552 0.0 0.0)
(0.2121701490094205 1.7667883953625454 0.0)
(0.39024121211934826 0.9255893082499055 1.1340772957831695))>
4:
#<size: 9
ave: #(2.0891122887721423 -4.286089167117592 -4.181378343344008)
std: #2A((1.2234257004280005 0.0 0.0)
(0.8126831105876011 3.770124644158256 0.0)
(-0.16095803539573883 3.355381122405471 1.6265434764655375))>
5:
#<size: 3
ave: #(0.5623597330105898 -5.258067550663203 -3.4228310177505143)
std: #2A((4.746711548715865 0.0 0.0)
(0.051149480910345976 6.110965690870317 0.0)
(1.7219389426928786 4.154449469768805 1.122266016458616))>
6:
#<size: 3
ave: #(0.8956775118811646 -1.4201650209349728 -1.6122538269908047)
std: #2A((1.493788623366322 0.0 0.0)
(-1.0842073850408864 0.5212017129303019 0.0)
(-0.0603674264384702 -0.011539750998482421 0.4925660868410117))>
7:
#<size: 3
ave: #(5.083108995741341 0.8389967015934188 -3.6196109709273347)
std: #2A((0.7952092465301981 0.0 0.0)
(-0.9765566485358852 0.6314489770374453 0.0)
(0.522986213867192 0.03536102846262183 0.8158863972930634))>
8:
#<size: 3
ave: #(5.229121525838935 -1.2777985422290163 -6.371216946124767)
std: #2A((2.9394585926746832 0.0 0.0)
(-0.7962469428473649 0.3415891910957655 0.0)
(1.4124306225879892 0.6782299557074762 0.9371629428713526))>
9:
#<size: 1
ave: #(-6.883331817688198 1.4950709172012338 3.3741067241587666)
std: #2A((0.7051887265707968 0.0 0.0)
(0.14539966262541817 1.2481210478667843 0.0)
(-0.31211910402732584 0.2711811196946303 1.2076541299769292))>
; cpu time (non-gc) 1.930000 sec user, 0.000000 sec system
; cpu time (gc) 1.480000 sec user, 0.000000 sec system
; cpu time (total) 3.410000 sec user, 0.000000 sec system
; real time 3.406327 sec
; space allocation:
; 4,353,779 cons cells, 526,686,512 other bytes, 0 static bytes
NIL
DPM(182): (dpm-hyper m)
2.343955996656361
DPM(183): (make-cluster-result m)
#(#(#(3.9823169099353253 -5.7538341155936195 9.533561089197962)
#(5.115651920711911 -4.88264142063046 9.122139238338159)
#(1.8776834245347507 -6.433568516594285 10.309460581627109)
#(3.5962489689683084 -4.637103542285794 9.532627083501511)
#(1.3400221559910137 -5.120247183311154 8.741750536406178)
#(4.102951740725052 -5.03882892508714 9.584405072336567)
#(1.4065757173228888 -7.10700132763071 10.558959618467927)
#(5.6100812402325335 -4.251787062174653 9.8578925415673)
#(1.6616682467100166 -4.391931453737993 10.059900574337714)
#(2.826214165794492 -4.549044606522388 9.412875645301698) ...)
#(#(4.569534672126397 -10.99770588007982 10.963515411778348)
#(4.446225550085373 -10.627826669753874 11.209657324492937)
#(5.059401983449978 -10.875220705854394 11.097027547127139)
#(3.774485958242961 -10.237199724529862 10.79342433539454)
#(4.805620159453174 -10.781424097277288 11.341697409559604)
#(4.492124617269321 -10.835153936892327 10.690929797020202)
#(5.066830269657376 -10.83840176135832 10.627420248371244)
#(4.828668092030972 -10.61614405126655 11.290629268822814)
#(5.091464774953241 -11.192689412544237 10.932270927572107)
#(4.7100255242843305 -10.641389270998047 10.65560545598822) ...)
#(#(-14.270841906056676 3.3791097630208267 2.5107616155787182)
#(-14.125488825924549 3.5242427620186447 -0.7217144548380698)
#(-14.26795451978224 3.316809855710163 -1.1928446471102845)
#(-11.681657831788186 3.0631516771560694 3.8662436706183088)
#(-13.379349530383099 2.8185233020395652 1.7508885911619503)
#(-12.844954470503376 2.698838177730742 1.4803285815966345)
#(-11.349353833583601 2.5171214604387964 3.222134775923579)
#(-13.521944292410254 2.5820067742123953 1.326359835873558)
#(-12.52938400148098 3.253060088686058 1.6547957934325) #(-10.96477091823868 2.948306749559725 6.07268720276314)
...)
#(#(4.035736076276529 -2.785943240870092 -5.561292237361581)
#(5.525263004657802 -0.7932348203208375 -7.402114421641164)
#(4.208814838416864 -2.257533315154677 -5.484791402964806)
#(3.993777930194894 -1.7698047546734597 -7.043994496920872)
#(2.6625531551899764 -2.50031272444097 -5.2359242214135735)
#(3.1612568006436055 -0.7121410527396044 -6.158077430120189)
#(3.448651702356854 -1.432867680128318 -6.377380648568867)
#(3.258527666960709 -2.9969096586814423 -4.850104473299948)
#(2.4539480572499683 -2.5474273338692672 -4.539388633150054)
#(4.317989899272528 -0.640867782623791 -7.321263913551618) ...)
#(#(1.12207771303666 -4.691118103798798 -3.830757992107001)
#(2.529954551285324 -4.157809175459205 -3.863175766611737)
#(1.8641471829555067 -3.647419630390355 -4.774547477725869)
#(1.5174351989759236 -3.061963170216118 -5.377208595762013)
#(3.4083304010130924 -4.8351350042145675 -3.861630938964044)
#(1.0505170887995634 -4.895694882479956 -3.5828556362300024)
#(2.802912795060661 -4.987111199235921 -3.0854817817179705)
#(1.3735923062069961 -4.125456973017023 -4.362775565106839)
#(1.7987815622391872 -4.217078197586187 -5.0076417981930055))
#(#(0.35230830640898514 -5.71983718318557 -2.5173806136127896)
#(0.36303962647870236 -5.155839967916859 -3.6002446291731744)
#(1.0361230585246095 -5.47157299728093 -3.6970440283426176))
#(#(-0.5591565012796771 -6.539237663213181 -2.5448159314421495)
#(0.21301664228906336 -6.05830454759155 -2.270224218174471)
#(-0.8675911178880398 -6.366410127153614 -2.314992343081113))
#(#(7.319176827407145 1.4086197926096182 -7.596042392086137)
#(5.725483006553445 -0.7008901859272663 -6.623625676629864)
#(6.334199278394568 0.5019862035471325 -7.317804977024041))
#(#(4.669415782342794 -0.9852399194597998 -6.258813500170089)
#(5.313765231930116 1.2602496804129917 -8.822966973690836)
#(5.431431720772352 -0.8757367528990956 -7.06611786833505))
#(#(-10.600144721004565 2.367613612996009 5.6164543818718595)) ...)
DPM(184):
備考：
毎回のクラスタリング結果は確率的で単一の正解とは言えないが、概ね4つの大きなクラスタが抽出されたことはMCMC列の平均を取ることで分かる
細かい結果を見ると、確かに分割されてしまっている小クラスタは別の固まりとみなせそうなものが得られていることも見て取れる
ここで表示されているstdは分散共分散行列の逆行列のコレスキー分解（1変数でいえば、標準偏差の逆数）にあたる
参考文献：
R. M. Neal, “Markov chain sampling methods for Dirichlet process mixture models,” University of Toronto, Toronto, ON, Canada, Tech. Rep. 9815, Sep. 1998
一変数正規分布で等分散を仮定した場合の丁寧なチュートリアル：
http://biocomp.bioen.uiuc.edu/journal_club_web/dirichlet.pdf
========================================
hdp-lda.lisp
HDP-LDA: topic数(k)を推定しながらトピック推定を行い文書-単語モデルを次元圧縮する
class:
- document
initargs:
id: 識別/印字用。計算には関係ない
words: stringのvector。一つの文書に複数回現れる単語は複数個入れる。hdp-ldaによって破壊的に変更される(下記wordの項参照)
- hdp-lda
initargs:
data: documentのvector
k: optional -- topic数 計算用配列の確保の目安にするだけ。kが可変なのでadjustable-arrayでのallocateになるが、その初期値
default -- 0
alpha: optional -- hyperparameterの初期値。最終的には推定されるが初期解に影響する。大きい程topic数は増える確率が高い
default -- (gamma-random *alpha-base-a* *alpha-base-b*)
gamma: optional -- おおむね同上。alphaとgammaの違いについては参考文献を参照されたい
default -- (gamma-random *gamma-base-a* *gamma-base-b*)
beta: optional -- 各トピックの事前分布を決めるパラメータ
default -- *default-beta*
struct:
- word
各単語は辞書からword-idを引き、関連付けられたCRF(Chienese Restaurant Franchise)のtableへの参照を持つ構造体へと変換される
methods:
- (initialize hdp-lda)
単語を数え上げたり配列を確保した後、初期解をサンプルする
戻り値は更新されたalpha,gamma
- (sampling hdp-lda)
sliced Gibbs samplingによって解を更新する
戻り値は更新されたalpha,gamma
- (assign-theta hdp-lda)
dataとして渡した各documentにthetaを割り当てる
thetaはトピック毎の混合比(トピックとの関連性・この文書内に生成されるトピックの事前確率)を表す
- (get-phi hdp-lda)
行列phiを返す
phiは語彙数V×kの行列であり、各トピックから生成される語の確率を示している
- (revert-word hdp-lda word-id)
word-id(=行列phiの添字)を元の単語表現に戻す
- (get-top-n-words hdp-lda n)
各トピックについて特徴的な頻出語上位n件を出力する
phiをトピック毎にソートしたものを語に戻しただけ
Accessor:
- (topic-count hdp-lda)
現在のkを返す
- (hdp-lda-data hdp-lda)
元データ(documentのvector)を返す
- (vocabulary hdp-lda)
語彙数Vを返す
- (document-id document)
- (document-words document)
- (document-thetas document)
- (word-id word)
読んで字の如し
下位APIについてはソース参照のこと
実行例:
HDP-LDA(5): (setf nips (with-open-file (in "nips.sexp") (read in)))
#(#("aaa" "abstract" "abu" "abu" "abu" "account" "achieve" "achieved" "acknowledgement" "acm" ...)
#("ability" "ability" "ability" "able" "abstract" "acad" "according" "acknowledgement" "activity" "activity" ...)
#("able" "abstract" "accuracy" "accurate" "achieved" "achieved" "achieved" "achieved" "achieved" "achieved" ...)
#("ability" "ability" "ability" "ability" "ability" "able" "absence" "abstract" "access" "accomplished" ...)
#("ability" "able" "able" "abstract" "acceptable" "according" "acknowledgment" "adding" "addison" "addition" ...)
#("accommodate" "accomplish" "acti" "activated" "activation" "activation" "activation" "activation" "activation"
"activation" ...)
#("able" "able" "abstract" "academic" "accept" "accept" "accept" "acceptance" "accepted" "accepted" ...)
#("aaa" "absence" "absence" "achieve" "actually" "addition" "addition" "additional" "additional" "adjust" ...)
#("abort" "absence" "abstract" "accepted" "accepting" "accomplished" "according" "according" "according" "account"
...)
#("ability" "ability" "abound" "abstract" "according" "according" "according" "accu" "accuracy" "accuracy" ...) ...)
HDP-LDA(6): (defun make-nips-dataset (nips)
(let ((id -1))
(map 'vector #'(lambda (x)
(make-instance 'document :id (incf id) :words (copy-seq x))) nips))) ;; wordsはcopyしておく
MAKE-NIPS-DATASET
HDP-LDA(7): (setf data (make-nips-dataset nips))
#(#<DOCUMENT 0 @ #x100d6356a2> #<DOCUMENT 1 @ #x100d6356d2> #<DOCUMENT 2 @ #x100d635702> #<DOCUMENT 3 @ #x100d635732>
#<DOCUMENT 4 @ #x100d635762> #<DOCUMENT 5 @ #x100d635792> #<DOCUMENT 6 @ #x100d6357c2> #<DOCUMENT 7 @ #x100d6357f2>
#<DOCUMENT 8 @ #x100d635822> #<DOCUMENT 9 @ #x100d635852> ...)
HDP-LDA(8): (setf m (make-instance 'hdp-lda :data data))
#<HDP-LDA @ #x10136fb632>
HDP-LDA(9): (time (initialize m))
; While compiling (:INTERNAL (:ANONYMOUS-LAMBDA 0) 0):
Warning: Free reference to undeclared variable M assumed special.
59063568 bytes have been tenured, next gc will be global.
See the documentation for variable EXCL:*GLOBAL-GC-BEHAVIOR* for more information.
87226336 bytes have been tenured, next gc will be global.
See the documentation for variable EXCL:*GLOBAL-GC-BEHAVIOR* for more information.
55039984 bytes have been tenured, next gc will be global.
See the documentation for variable EXCL:*GLOBAL-GC-BEHAVIOR* for more information.
; cpu time (non-gc) 102.840000 sec (00:01:42.840000) user, 0.140000 sec system
; cpu time (gc) 49.570000 sec user, 0.020000 sec system
; cpu time (total) 152.410000 sec (00:02:32.410000) user, 0.160000 sec system
; real time 152.611516 sec (00:02:32.611516)
; space allocation:
; 870,573 cons cells, 82,959,059,408 other bytes, 0 static bytes
150.27874398522025d0
27.410756602369037d0
HDP-LDA(10): (topic-count m)
192
HDP-LDA(11): (time (sampling m))
; While compiling (:INTERNAL (:ANONYMOUS-LAMBDA 1) 0):
Warning: Free reference to undeclared variable M assumed special.
; cpu time (non-gc) 112.670000 sec (00:01:52.670000) user, 0.500000 sec system
; cpu time (gc) 44.310000 sec user, 0.260000 sec system
; cpu time (total) 156.980000 sec (00:02:36.980000) user, 0.760000 sec system
; real time 157.756031 sec (00:02:37.756031)
; space allocation:
; 3,306 cons cells, 91,564,610,016 other bytes, 0 static bytes
367.6535745173031d0
24.501190724062724d0
HDP-LDA(12): (topic-count m)
204
HDP-LDA(13): (get-top-n-words m 10)
#(#("recurrent" "motion" "recognition" "neural" "affine" "euclidean" "mlp" "ica" "cholinergic" "top")
#("network" "set" "result" "output" "error" "number" "weight" "method" "order" "layer")
#("model" "unit" "optimal" "task" "probability" "system" "vector" "high" "net" "condition")
#("noise" "class" "prediction" "separation" "form" "student" "natural" "source" "linear" "tangent")
#("vector" "point" "function" "unit" "weight" "test" "training" "space" "error" "term")
#("problem" "word" "local" "temporal" "trained" "rule" "respect" "update" "transition" "active")
#("feature" "features" "density" "object" "rbf" "size" "search" "hidden" "large" "covariance")
#("neural" "input" "algorithm" "function" "training" "parameter" "image" "control" "information" "linear")
#("neuron" "function" "gaussian" "field" "mean" "small" "case" "correlation" "iteration" "energy")
#("system" "function" "response" "distribution" "images" "vector" "hmm" "order" "sample" "application") ...)
HDP-LDA(14):
参考文献:
Hierarchical Dirichlet Processes, Yee Whye Teh, Michael I Jordan, Matthew J Beal, David M Blei. Journal of the American Statistical Association. December 1, 2006, 101(476): 1566-1581.
備考:
参考文献内の"Chinese Restaurant Franchise"が実装されている
理論上はDPMを使えば実装できるが、効率のためにそのような実装はしていない
収束性を判断するのは難しいが、k、hyperparameters、および(hdp-lda-topic-occurs hdp-lda)の値が収束するのがもっとも全体の収束に近いと思われる
しかし、毎回サンプリングされるため、収束といっても振動は残ることになる
get-top-n-wordsで単語を得る際、そのトピックに属すと推定された語がn未満の場合、単語のid順に語が埋められてしまう
実際には残りの単語がすべて均等な確率なので、必ずしもバグとは言えないが混乱の元とは言える
========================================
hdp.lisp
Hierachical Dirichlet Process(HDP): 階層Dirichlet過程のライブラリ（DPMの拡張）
Direct Assignment表現を用いることでChinese Restaurant Franchiseより省スペース・高速に動作する（はず）
class:
- hdp-cluster << cluster
- hdp << dpm
- hdp-distribution << dp-distribution
Accessor:
- (hdp-gamma hdp)
HDPのgamma（上位DPのハイパーパラメータ）を返す
(dpm-hyper dpm)は下位DPのハイパーパラメータを返す
- (cluster-beta hdp-cluster)
Dirichlet分布からサンプルされたクラスタの事前分布\beta の現在値を返す
- (hdp-beta hdp)
Dirichlet分布からサンプルされた新規クラスタの事前分布\beta の現在値を返す
参考文献:
Hierarchical Dirichlet Processes, Yee Whye Teh, Michael I Jordan, Matthew J Beal, David M Blei. Journal of the American Statistical Association. December 1, 2006, 101(476): 1566-1581.
========================================
hdp-hmm.lisp
HDP-HMM (a.k.a. infinite-HMM): 隠れ状態の数kを推定しながら教師なし隠れマルコフモデルを学習する
離散データに対するHMMであり、事前分布は一様分布を仮定している
全体を一つの大きな系列とみなし、各データ点を一つずつサンプルするGibbsを構成する
class:
- hidden-state << hdp-cluster
- hdp-hmm << hdp
- state-uniform << hdp-distribution
methods:
- trans-prob (before after &rest args &key &allow-other-keys)
隠れ状態beforeから隠れ状態afterへ遷移する確率を返す
defaultは多項分布による最尤推定
- emission-prob (state data &rest args &key &allow-other-keys)
隠れ状態stateからdataが生成される確率を返す
defaultは多項分布による最尤推定
- show-hidden-state (hdp-hmm)
与えたデータ列に対して推定された状態列を頻度順idで置換して返す
参考文献:
Hierarchical Dirichlet Processes, Yee Whye Teh, Michael I Jordan, Matthew J Beal, David M Blei. Journal of the American Statistical Association. December 1, 2006, 101(476): 1566-1581.
========================================
gauss-hmm.lisp
データ生成分布に正規分布を仮定した連続データ用HDP-HMM
HDP-HMMの拡張例
gauss-dpmとhdp-hmmを多重継承することでコードを大部分再利用
emission-probとinitializeの書き換えのみ
========================================
sticky-hdp-hmm.lisp
Sticky HDP-HMM: （特に連続データに対する）HDP-HMMが状態数を過剰に生成してまうことを抑制し、sticky-parameterを推定する
class:
- sticky-hdp-hmm << hdp-hmm
- sticky-hidden-state << hidden-state
- sticky-state-unifrom << state-uniform
Accessor:
- sticky-kappa (sticky-hdp-hmm)
sticky-parameterを返す
HDP-HMMの拡張に対してmix-inとして組み込むことでsticky-parameterを導入できる
詳細は参考文献に譲る
参考文献：
An HDP-HMM for systems with state persistence, Emily B. Fox, Erik B. Sudderth, Michael I. Jordan and Alan S. Willsky. Proceedings of the 25th international conference on Machine learning. 2008, 312-319
========================================
blocked-hdp-hmm.lisp
Blocked HDP-HMM: 系列の前後依存性を考慮し効率的なMCMCを行うため、一つの列を単位としてBlocked Gibbs Samplingを行うHDP-HMM
class
- seq-point << point
系列内の各点を表す
- point-sequence
「系列」を表す
initargs:
seq: seq-pointのvector
- blocked-hdp-hmm << hdp-hmm
initargs:
data: point-sequenceのvector
L: 現在の状態数の上限（Blocked Gibbsを行うために必要）
初期値を与えればあとは必要に合わせて増やすが、極端に小さい値を与えると局所解に落ちやすい傾向がある
- blocked-hidden-state << hidden-state
- block-uniform << state-uniform
Accessor:
- sequence-data (point-sequence)
seq-pointのvectorを返す
- hdp-hmm-l (blocked-hdp-hmm)
現在の状態数の上限Lを返す
実際の状態数kとは異なるが、常にkより大きい
参考文献：
Beam Sampling for the Infinite Hidden Markov Model, Jurgen Van Gael, Yunus Saatci, Yee Whye Teh, and Zoubin Ghahramani.
ICML 2008
========================================
ihmm.lisp
blocked & sticky HDP-HMM: 基礎的に使うべき万能？HDP-HMM これをここではiHMMと呼ぶことにする
class:
- ihmm << blocked-hdp-hmm sticky-hdp-hmm
- ihmm-state << blocked-hidden-state sticky-hidden-state
- ihmm-state-uniform << block-uniform sticky-state-uniform
sticky featureをBlocked Gibbsに対応させるために上書きすべきmethodがあるため、単純なmix-inとはならない
- gauss-ihmm << ihmm gauss-hdp-hmm
- gauss-ihmm-state << ihmm-state gaussian-state
- ihmm-gaussian << ihmm-state-uniform state-gaussian
実行例:
IHMM(14): (setf (values trans1 data1) (make-cyclic-test 500 #(0d0 50d0 0d0 -50d0)))
#(0 0 0 0 0 0 0 0 0 0 ...)
#<POINT-SEQUENCE @ #x1000bd9d72>
IHMM(15): (setf (values trans2 data2) (make-cyclic-test 500 #(0d0 50d0 0d0 -50d0)))
#(0 0 0 0 0 0 0 0 0 0 ...)
#<POINT-SEQUENCE @ #x1000c3f9f2>
IHMM(16): (setf trans (coerce (list trans1 trans2) 'vector))
#(#(0 0 0 0 0 0 0 0 0 0 ...) #(0 0 0 0 0 0 0 0 0 0 ...))
IHMM(17): (setf data (coerce (list data1 data2) 'vector))
#(#<POINT-SEQUENCE @ #x1000bd9d72> #<POINT-SEQUENCE @ #x1000c3f9f2>)
IHMM(18): (map 'vector #'(lambda (x) (map 'vector #'(lambda (y) (point-data y)) (sequence-data x))) data) ;; 生データ確認
#(#(0.7774628524351572d0 -1.317263393017809d0 0.1499122993576713d0 1.760687062844291d0 0.5890700040616544d0
0.7620681604876274d0 0.20820480252339166d0 0.008305963440786393d0 0.6280792673787637d0 -0.3690057889528494d0 ...)
#(-1.2589573109025902d0 -0.6680537631049001d0 -0.977360542281873d0 -0.9204817444451175d0 0.1444950521775598d0
0.7507011449925263d0 -0.3055258644799465d0 -0.6574616916707324d0 -0.9588044829834967d0 0.45389443148718195d0 ...))
IHMM(26): (setf m (make-instance 'gauss-ihmm :data data :l 10)) ;; 状態数上限は真の状態4より多めに設定
#<GAUSS-IHMM @ #x10011b0592>
IHMM(27): (initialize m)
0.041013719237274926d0 ;; alpha (下層DPのhyper-parameter)
4.577657553892466d0 ;; gamma (上層DPのhyper-parameter)
0.8890828011820707d0 ;; kappa (sticky-parameter)
IHMM(28): (dotimes (i 100) (sampling m))
NIL
IHMM(29): (pprint (show-hidden-states m)) ;; 状態推定結果を表示
#(#(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 2 2 2 4 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 1 1 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
#(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 3 3 3 4 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 4 1 1 1 1 1 1 1 1 1 1 1 5 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 6 6 6 6 6))
IHMM(30): (pprint trans) ;; データ生成時に使った真の状態
#(#(0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2)
#(0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0))
;;; 状態遷移時に中間状態が生成されていることを除き、ほぼ正しく推定されている
;;; 真の状態で0と2は共に平均0d0の正規分布であるがHDP-HMMが二つを正しく分離していることに注意
IHMM(31): (dpm-k m)
7
IHMM(32): (hdp-hmm-l m)
11
IHMM(33):