HFM for R!

やっと少し時間が取れたので、RでHFMをやってみようと思いまして。関数を作りました。次の通りです。

hfm<span class="synStatement">&lt;-</span><span class="synType">function</span><span class="synSpecial">(</span>data<span class="synSpecial">)</span>
<span class="synSpecial">{</span>
<span class="synStatement">if</span><span class="synSpecial">(</span>!is.matrix<span class="synSpecial">(</span>data<span class="synSpecial">))</span>            <span class="synComment">#行列形式でなければ行列形式にしてしまう</span>
data <span class="synStatement">&lt;-</span> as.matrix<span class="synSpecial">(</span>data<span class="synSpecial">)</span>
<span class="synStatement">if</span><span class="synSpecial">(</span>ncol<span class="synSpecial">(</span>data<span class="synSpecial">)</span>!=nrow<span class="synSpecial">(</span>data<span class="synSpecial">))</span>      <span class="synComment">#正方行列かどうかのチェック</span>
stop<span class="synSpecial">(</span><span class="synConstant">&quot;data is not a square matrix&quot;</span><span class="synSpecial">)</span>
<span class="synComment">#エルミート化</span>
s <span class="synStatement">&lt;-</span> <span class="synSpecial">(</span>data+t<span class="synSpecial">(</span>data<span class="synSpecial">))</span>/<span class="synConstant">2</span>
sk <span class="synStatement">&lt;-</span> <span class="synSpecial">(</span>data-t<span class="synSpecial">(</span>data<span class="synSpecial">))</span>/<span class="synConstant">2</span>
Her <span class="synStatement">&lt;-</span> <span class="synType">matrix</span><span class="synSpecial">(</span><span class="synType">complex</span><span class="synSpecial">(</span>re=s<span class="synSpecial">,</span>im=sk<span class="synSpecial">),</span>ncol=ncol<span class="synSpecial">(</span>s<span class="synSpecial">))</span>
rownames<span class="synSpecial">(</span>Her<span class="synSpecial">)</span> <span class="synStatement">&lt;-</span> colnames<span class="synSpecial">(</span>data<span class="synSpecial">)</span>
colnames<span class="synSpecial">(</span>Her<span class="synSpecial">)</span> <span class="synStatement">&lt;-</span> colnames<span class="synSpecial">(</span>data<span class="synSpecial">)</span>
<span class="synComment">#<a class="keyword" href="http://d.hatena.ne.jp/keyword/%B8%C7%CD%AD%C3%CD">固有値</a>分解</span>
eval <span class="synStatement">&lt;-</span> eigen<span class="synSpecial">(</span>Her<span class="synSpecial">)</span>$values
evec <span class="synStatement">&lt;-</span>eigen<span class="synSpecial">(</span>Her<span class="synSpecial">)</span>$vectors
rownames<span class="synSpecial">(</span>evec<span class="synSpecial">)</span> <span class="synStatement">&lt;-</span> colnames<span class="synSpecial">(</span>data<span class="synSpecial">)</span>
<span class="synComment">#絶対値の順に並べ替え</span>
od <span class="synStatement">&lt;-</span> abs<span class="synSpecial">(</span>eval<span class="synSpecial">)</span>
eval <span class="synStatement">&lt;-</span> eval<span class="synSpecial">[</span>order<span class="synSpecial">(</span>od<span class="synSpecial">,</span>decreasing=<span class="synConstant">TRUE</span><span class="synSpecial">)]</span>
evec <span class="synStatement">&lt;-</span> evec<span class="synSpecial">[,</span>order<span class="synSpecial">(</span>od<span class="synSpecial">,</span>decreasing=<span class="synConstant">TRUE</span><span class="synSpecial">)]</span>
<span class="synComment">#寄与率の算出</span>
<a class="keyword" href="http://d.hatena.ne.jp/keyword/gof">gof</a>   <span class="synStatement">&lt;-</span>eval*eval / <span class="synSpecial">(</span>t<span class="synSpecial">(</span>eval<span class="synSpecial">)</span>%*%eval<span class="synSpecial">)</span>
<span class="synStatement">return</span><span class="synSpecial">(</span><span class="synType">list</span><span class="synSpecial">(</span>Hermitian =Her <span class="synSpecial">,</span>Eigen=eval<span class="synSpecial">,</span><a class="keyword" href="http://d.hatena.ne.jp/keyword/GOF">GOF</a>=<a class="keyword" href="http://d.hatena.ne.jp/keyword/gof">gof</a><span class="synSpecial">,</span>Vecs=evec<span class="synSpecial">))</span>
<span class="synError">}</span>
<a class="keyword" href="http://d.hatena.ne.jp/keyword/hfm">hfm</a>.plot <span class="synStatement">&lt;-</span><span class="synType">function</span><span class="synSpecial">(</span>data<span class="synSpecial">,</span>dim<span class="synSpecial">,</span>Xlim=c<span class="synSpecial">(</span>-<span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">),</span>Ylim=c<span class="synSpecial">(</span>-<span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">))</span>
<span class="synSpecial">{</span>
plot<span class="synSpecial">(</span>data$Vecs<span class="synSpecial">[,</span>dim<span class="synSpecial">],</span>
main=paste<span class="synSpecial">(</span><span class="synConstant">&quot;Dim&quot;</span><span class="synSpecial">,</span> dim<span class="synSpecial">,</span> <span class="synConstant">&quot;with Eigenvalue&quot;</span> <span class="synSpecial">,</span>round<span class="synSpecial">(</span>data$Eigen<span class="synSpecial">[</span>dim<span class="synSpecial">],</span><span class="synConstant">3</span><span class="synSpecial">)),</span>
xlab=<span class="synConstant">&quot;real&quot;</span><span class="synSpecial">,</span>    <span class="synComment">#Xラベル</span>
ylab=<span class="synConstant">&quot;imag&quot;</span><span class="synSpecial">,</span>  <span class="synComment">#Yラベル</span>
xlim=Xlim<span class="synSpecial">,</span>  <span class="synComment">#X範囲</span>
ylim=Ylim<span class="synSpecial">,</span>  <span class="synComment">#Y範囲</span>
axes=F<span class="synSpecial">)</span>
axis<span class="synSpecial">(</span><span class="synConstant">1</span><span class="synSpecial">,</span> pos = <span class="synConstant">0</span><span class="synSpecial">,</span> at = -<span class="synConstant">3</span><span class="synSpecial">:</span><span class="synConstant">3</span><span class="synSpecial">,</span> adj = <span class="synConstant">0</span><span class="synSpecial">,</span> col = <span class="synConstant">1</span><span class="synSpecial">)</span> <span class="synComment">#  X 軸を描く</span>
axis<span class="synSpecial">(</span><span class="synConstant">2</span><span class="synSpecial">,</span> pos = <span class="synConstant">0</span><span class="synSpecial">,</span> at = -<span class="synConstant">3</span><span class="synSpecial">:</span><span class="synConstant">3</span><span class="synSpecial">,</span> adj = <span class="synConstant">1</span><span class="synSpecial">,</span> las = <span class="synConstant">2</span><span class="synSpecial">)</span> <span class="synComment">#  Y 軸を描く</span>
<span class="synStatement">for</span><span class="synSpecial">(</span> i <span class="synStatement">in</span> seq<span class="synSpecial">(</span>along = data$Vecs<span class="synSpecial">[,</span>dim<span class="synSpecial">]))</span>
arrows<span class="synSpecial">(</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span>Re<span class="synSpecial">(</span>data$Vecs<span class="synSpecial">[</span>i<span class="synSpecial">,</span>dim<span class="synSpecial">]),</span>Im<span class="synSpecial">(</span>data$Vecs<span class="synSpecial">[</span>i<span class="synSpecial">,</span>dim<span class="synSpecial">]))</span>
text<span class="synSpecial">(</span>data$Vecs<span class="synSpecial">[,</span>dim<span class="synSpecial">],</span>rownames<span class="synSpecial">(</span>data$Vecs<span class="synSpecial">))</span>
<span class="synSpecial">}</span>

これで、例えば先日の北京五輪、野球一次予選のデータを分析してみよう。8チームの総当たり、非対称行列なのでHFMの例示にもってこいですな。
次のプログラムを実行してください。

H<span class="synStatement">&lt;-</span><span class="synType">matrix</span><span class="synSpecial">(</span>c<span class="synSpecial">(</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">5</span><span class="synSpecial">,</span><span class="synConstant">7</span><span class="synSpecial">,</span><span class="synConstant">6</span><span class="synSpecial">,</span><span class="synConstant">8</span><span class="synSpecial">,</span><span class="synConstant">2</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">6</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">3</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">8</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">5</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">10</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">6</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">9</span><span class="synSpecial">,</span><span class="synConstant">10</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">8</span><span class="synSpecial">,</span><span class="synConstant">7</span><span class="synSpecial">,</span><span class="synConstant">5</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">7</span><span class="synSpecial">,</span><span class="synConstant">9</span><span class="synSpecial">,</span><span class="synConstant">5</span><span class="synSpecial">,</span><span class="synConstant">7</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">14</span><span class="synSpecial">,</span><span class="synConstant">17</span><span class="synSpecial">,</span><span class="synConstant">7</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">5</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">6</span><span class="synSpecial">,</span><span class="synConstant">6</span><span class="synSpecial">,</span><span class="synConstant">10</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">,</span><span class="synConstant">3</span><span class="synSpecial">,</span><span class="synConstant">2</span><span class="synSpecial">,</span><span class="synConstant">2</span><span class="synSpecial">,</span><span class="synConstant">0</span><span class="synSpecial">),</span>nrow=<span class="synConstant">8</span><span class="synSpecial">)</span>
colnames<span class="synSpecial">(</span>H<span class="synSpecial">)</span> <span class="synStatement">&lt;-</span> c<span class="synSpecial">(</span><span class="synConstant">&quot;Tiwan&quot;</span><span class="synSpecial">,</span><span class="synConstant">&quot;Nederland&quot;</span><span class="synSpecial">,</span><span class="synConstant">&quot;China&quot;</span><span class="synSpecial">,</span><span class="synConstant">&quot;Canada&quot;</span><span class="synSpecial">,</span><span class="synConstant">&quot;Korea&quot;</span><span class="synSpecial">,</span><span class="synConstant">&quot;America&quot;</span><span class="synSpecial">,</span><span class="synConstant">&quot;Cuba&quot;</span><span class="synSpecial">,</span><span class="synConstant">&quot;Japan&quot;</span><span class="synSpecial">)</span>
result<span class="synStatement">&lt;-</span><a class="keyword" href="http://d.hatena.ne.jp/keyword/hfm">hfm</a><span class="synSpecial">(</span>H<span class="synSpecial">)</span>
<a class="keyword" href="http://d.hatena.ne.jp/keyword/hfm">hfm</a>.plot<span class="synSpecial">(</span>result<span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">,</span>c<span class="synSpecial">(</span>-<span class="synConstant">0.5</span><span class="synSpecial">,</span><span class="synConstant">0.5</span><span class="synSpecial">),</span>c<span class="synSpecial">(</span>-<span class="synConstant">0.5</span><span class="synSpecial">,</span><span class="synConstant">0.5</span><span class="synSpecial">))</span>

上二行がデータ入力、三行目がhfm関数の実行。四行目はプロット。より詳しい解説は、本家サイトKosugitti Laboにreadme等といっしょにパッケージングしたものを置いてありますので、そちらからDLして使ってみてください。
結果の一部を書きますと、

$Eigen
<span class="synSpecial">[</span><span class="synConstant">1</span><span class="synSpecial">]</span>  <span class="synConstant">35.598043</span> -<span class="synConstant">17.753666</span> -<span class="synConstant">10.926530</span>  -<span class="synConstant">8.851190</span>   <span class="synConstant">6.448190</span>  -<span class="synConstant">5.022845</span>
<span class="synSpecial">[</span><span class="synConstant">7</span><span class="synSpecial">]</span>   <span class="synConstant">2.113688</span>  -<span class="synConstant">1.605690</span>
$<a class="keyword" href="http://d.hatena.ne.jp/keyword/GOF">GOF</a>
<span class="synSpecial">[</span><span class="synConstant">1</span><span class="synSpecial">]</span> <span class="synConstant">0.683506285</span> <span class="synConstant">0.170006836</span> <span class="synConstant">0.064395398</span> <span class="synConstant">0.042256505</span> <span class="synConstant">0.022426731</span> <span class="synConstant">0.013607858</span>
<span class="synSpecial">[</span><span class="synConstant">7</span><span class="synSpecial">]</span> <span class="synConstant">0.002409751</span> <span class="synConstant">0.001390637</span>

となる。第一次元の寄与率が68%、第二次元が17%。二次元目までで85%だから、どこまで使うかが考え処なんだけど、それよりも注目すべきは固有値の符号。正と負が入り乱れてます。これはつまり、データ全体にある単位（ノルム）を決めるのが難しい*1ということ。語弊を恐れずに簡略化して言えば、一つのモノサシで決められるほどはっきりした構造を保ってないということです。

第一次元をプロットしたのが冒頭の図。台湾を軸に、上と下に分かれています。非対称関係なので、上が台湾より強い、下が台湾より弱い、と考えてもらえばいいでしょう。原点から一番遠いのがキューバ。キューバは8試合で52点とって23点取られてます。原点からの距離で言えば、次は韓国かアメリカか。同じぐらいのポジションです。中身を見てみると、韓国は41点取って22点取られてる。アメリカは40点とって23点取られてる。しかも、韓国vsアメリカは8-7で韓国の勝ち。この辺が図の微妙な位置関係に、ちゃんと反映されているのが見て取れます。

星野さんは「勝ったものが強い」という迷言を吐いてましたが、これを見ると日本は金メダルがとれようはずもなく、予選の結果だけで行けばキューバが金なのです。実際は韓国が金、キューバが銀、アメリカ銅でした。上位２チームが入れ替わったのは、いつかも述べたように、また分析結果が綺麗な距離空間に落とし込めなかったように、尺度が安定しない、あるいは、勝負のアヤというのが効いてくる競技だからでしょう。