pHomパッケージでドーナッツの穴を数えてみた
CRANに登録されているpersistent homologyを計算できるpHomパッケージで3次元トーラスの穴を数えてみました。
ソースはこちら。
データ生成
まず与える人工データの生成ですが、極座標表示したときの回転パラメータ2つでそれぞれ20間隔にして20x20の400点作り、ドーナッツの直径が50に対して標準偏差1のガウスノイズを各点ごとに独立に乗っけてみました。下で呼ぶpHom関数はデフォルトでは1000以上の点がある場合は1000だけ適当にピックアップして動作するようです。なので400なら特に何も言わなくても全部を考慮してくれる。
実行結果
pHomパッケージに与えるのはフィルトレーション(点たちを見るときの解像度)の最大値と、調べたいホモロジーの次元の最大値で、それぞれ5, 2としています。トーラスにおいては、0次元ホモロジー(点のあつまりを点とみなすときの点)、1次元ホモロジー(点のあつまりを線とみなすときの輪)、2次元ホモロジー(点のあつまりを面とみなすときの球)がそれぞれ1,2,1となるのが自然です。こういう輪っかの数のことを各次元のベッチ数と呼ぶそうです。ベッチ、ベッチ、いそべっち。いやすみません。
その条件でホモロジー計算させた結果。バーコードとパーシステントダイアグラム。
b0の図が0次元ベッチ数のバーコードで、横軸のフィルトレーションが1.5あたりで一気に一個に集約されてます。
b1のほうはその切り替わる1.5あたりで輪っかを一つ認識していて、2.3あたりからもう一つの輪っかがリレー式に出来てます。これは多分解像度の調整がうまくいってないせいで、ドーナツを切り分けた断片に相当するパイプチューブみたいな形の組み合わせが次第に集約されていってると認識してるんじゃないかと思います。
そして残念な事にb2つまりドーナツの中の空洞は認識できませんでした。パラメータ設定を色々変えると、解像度がある一定を超えると一気に計算量が増えるのか、pHomの呼び出しが帰ってこなくなります。この辺りのアルゴリズムの計算量のことをもうちょっと学ばないと。
下のパーシステントダイアグラムのほうは上の各バーコードの両端を2次元にプロットしたものなので意味するところは一緒です。