2011-10-10

Haskellで蟻本を解く1

Haskell

最近Haskellを勉強していて、練習がてら、蟻本こと「プログラミングコンテストチャレンジブック」の問題をできるだけ解いてみようと思います。
最初は練習問題の

n本の棒から周長が最大となる三角形を作る

という問題。これをO(nlogn)で解くにはどうしたらいいかというものです。(しらみつぶしでO(n^3)の方法については本に書いてあります)

qsort :: [Int] -> [Int]
qsort [] = []
qsort (x:xs) =
  qsort [x' | x' <- xs, x'>x] ++ [x] ++ qsort [x' | x' <- xs, x'<=x]

maxTriangleSorted :: [Int] -> [Int]
maxTriangleSorted (x:xs)
  | length xs <	2 = []
  | x <	sum xss = [x] ++ xss
  | otherwise = maxTriangleSorted xs
                where xss = take 2 xs

maxTriangle :: [Int] ->	[Int]
maxTriangle = maxTriangleSorted.qsort

maxTriangle [15,2,7,3,8,6]
[8,7,6]

qsortはO(nlogn)で比較は高々n回なので計算時間はO(nlogn)です。
あまりHaskell的にはよろしくないコードだと思いますが…。

プログラミングコンテストチャレンジブック

作者: 秋葉拓哉,岩田陽一,北川宜稔
出版社/メーカー: 毎日コミュニケーションズ
発売日: 2010/09/11
メディア: 単行本（ソフトカバー）
購入: 52人クリック: 1,538回
この商品を含むブログ (83件) を見る

2011-09-10

homebrewからopencvの最新版(2.3.1)をインストールする

cv mac

そのうち必要なくなるとは思いますが、現在のhomebrewのデフォルトではopencvのformulaは2.2をインストールするように記述されているため、最新版を入れるためには少々作業が必要でした。
2.3以降を入れる大きなメリットとして、新しいpython バインディングcv2が使えるようになるというものがあります。pyhtonからopencvを使いたい場合は、なるべく2.3以上にしておきたいところだと思います。

インストールのための作業

インストール作業は以下のようにしました。

opencvのFormulaを書き換える
monoが悪さをするので(インストールされている場合は)アンインストールする
brew upgradeする(またはbrew installする)

opencvのFormulaを書き換える

https://github.com/wjwwood/homebrew/commit/9317f433927a3a7ef9dc40d147830c366ec9b0a7
この通りにopencv.rbを修正するか、あるいはgithubにcommitされた最新版を落としてきます。
（brew updateで反映されれば良いのですが、僕の環境では反映されませんでした）
ローカル環境のFormulaはおそらく/usr/local/Library/Formula以下にあります。

monoが悪さをするので(インストールされている場合は)アンインストールする

この状態でいきなりbrew upgradeをすると失敗することがあると思います。試行錯誤した結果、monoがインストールされていることが原因であるらしいことがわかりました。とりあえずアンインストールして続行します（opencvのインストール後にmonoを入れ直せばおそらく大丈夫だと思います）。
monoのアンインストール方法は、
http://www.mono-project.com/Mono:OSX#Uninstalling_Mono_on_Mac_OS_X
にあるUninstalling Mono on Mac OS Xに従ってシェルスクリプトを実行すればOKです。

brew upgradeする(またはbrew installする)

これでbrew upgradeすれば最新版に更新されると思いますが、僕の環境ではバージョンの管理がうまく行かなかったようで、そのままでは旧バージョンとの共存状態になってしまいました。(homebrewをきちんと理解してないためかもしれません…)
というわけで、まずは

brew remove opencv

で既存のopencvをアンインストールし、改めて

brew install opencv

を実行します。

以上で最新版のopencvがインストールされると思います。

2011-08-14

pythonでパーティクルフィルタによる物体トラッキング

機械学習・統計 python cv

Tokyo.SciPyという勉強会に参加するので、pythonの練習に最近少し勉強しているコンピュータビジョン関連のアルゴリズムを実装してみました。
パーティクルフィルタによる物体トラッキングです。OpenCVを用いてカメラからキャプチャした画像の色情報を使い、追跡したい物体にパーティクルが集まるように尤度関数を設定しています。よく例題として出てくるものです。一応、システムモデル（パーティクルの時間発展を記述）、尤度モデル（パーティクルの重みを記述）を表す関数オブジェクトを入れ替えることで色々な問題に対応できるようにはなっています。練習ということであまりこだわらずにアドホックに書いてしまったところもあり、特にOpenCV周りはいい加減です…。

# -*- coding:utf-8 -*-

import random
import numpy as np
from scipy import stats
import cv

#キャプチャを管理
class Image:
    def __init__(self):
        self.capture = cv.CreateCameraCapture(0)
        self.image = cv.QueryFrame(self.capture)
        cv.ShowImage("Capture",self.image)
        self.size = (self.image.width,self.image.height)

    def create(self):
        self.image = cv.QueryFrame(self.capture)

    def getCol(self,sv):
        x = sv[0]
        y = sv[1]
        #パーティクルが画面外にはみ出した場合は黒として返す
        if((x<0 or x>self.size[0]) or (y<0 or y>self.size[1])):
            return((0,0,0,0))
        else:
            #OpenCVでは座標指定が(y,x)の順
            return(cv.Get2D(self.image,int(sv[1]),int(sv[0])))

#システムモデルを管理
class SystemModel:
    def __init__(self,model):
        self.model = model

    def generate(self,sv,w):
        return(self.model(sv,w))

#尤度モデルを管理
class Likelihood:
    def __init__(self,model):
        self.model = model

    def generate(self,sv,mv):
        return(self.model(sv,mv))

    def normalization(self,svs,mv):
        return(sum([self.generate(sv,mv) for sv in svs]))

#システムモデルの関数オブジェクト
def model_s(sv,w):
    #等速直線運動モデル
    #状態ベクトルは(x,y,vx,vy)を仮定
    F = np.matrix([[1,0,1,0],
                   [0,1,0,1],
                   [0,0,1,0],
                   [0,0,0,1]])
    return(np.array(np.dot(F,sv))[0]+w)

#尤度モデルの関数オブジェクト
def model_l(sv,mv):
    #mv(cvイメージ)からsvの座標値(第1,2成分)に
    #対応する点の色情報を取得する
    mv_col = img.getCol(sv)
    mv_col = mv_col[0:3]
    target_col = (43,22,55) #追跡したい物体の色(BGR)

    #尤度は色情報と指定する色の差に対するガウスモデル
    delta = np.array(mv_col)-np.array(target_col)
    dist_sqr = sum(delta*delta)
    sigma = 10000.0
    gauss = np.exp(-0.5*dist_sqr/(sigma*sigma)) / (np.sqrt(2*np.pi)*sigma)
    return(gauss)

#パーティクルのリサンプリング
def resampling(svs,weights):
    N = len(svs)
    #重みの大きい順にパーティクルをソート
    sorted_particle = sorted([list(x) for x in zip(svs,weights)],key=lambda x:x[1],reverse=True)
    #重みに従ってパーティクルをリサンプリング
    resampled_particle = []
    while(len(resampled_particle)<N):
        for sp in sorted_particle:
            resampled_particle += [sp[0]]*(sp[1]*N)
    resampled_particle = resampled_particle[0:N]

    return(resampled_particle)

#1ステップ分のフィルタ
def filtration(svs,mv,systemModel,likelihood):
    #システムモデルに用いる乱数を生成
    dim = len(svs[1])
    N = len(svs)
    sigma = 5.0
    rnorm = stats.norm.rvs(0,sigma,size=N*dim)
    ranges = zip([N*i for i in range(dim)],[N*i for i in (range(dim+1)[1:])])
    ws = np.array([rnorm[p:q] for p,q in ranges])
    ws = ws.transpose()

    #予測サンプルを生成
    svs_predict = [systemModel.generate(sv,w) for sv,w in zip(svs,ws)]
    
    #尤度重みを計算
    normalization_factor = likelihood.normalization(svs_predict,mv)
    likelihood_weights = [likelihood.generate(sv,mv)/normalization_factor for sv in svs_predict]
    #重みによってリサンプリング
    svs_resampled = resampling(svs_predict,likelihood_weights)
    return(svs_resampled)

#初期パーティクルを生成する。モデルが(x,y,vx,vy)の4次元モデルであることを仮定
def initStateVectors(imageSize,sampleSize):
    xs = [random.uniform(0,imageSize[0]) for i in range(sampleSize)]
    ys = [random.uniform(0,imageSize[1]) for i in range(sampleSize)]
    vxs = [random.uniform(0,5) for i in range(sampleSize)]
    vys = [random.uniform(0,5) for i in range(sampleSize)]

    return([list(s) for s in zip(xs,ys,vxs,vys)])

#パーティクル付きの画像を表示する
def showImage(svs,img):
    #パーティクル描画用のコピー
    dst = cv.CloneImage(img.image)
    #パーティクルを書き込み
    for sv in svs:
        #パーティクル位置をintに変換
        cv.Circle(dst,(int(sv[0]),int(sv[1])),3,cv.CV_RGB(0,0,255))
    #表示
    cv.ShowImage("Capture",dst)


if(__name__=="__main__"):
    #イメージソースを指定
    img = Image()
    #パーティクル数を指定
    sampleSize = 100
    #モデルオブジェクトを生成
    systemModel = SystemModel(model_s)
    likelihood = Likelihood(model_l)
    #初期パーティクルを生成
    svs = initStateVectors(img.size,sampleSize)

    while(True):
        #描画
        showImage(svs,img)
        #観測
        img.create()
        #フィルタ
        svs = filtration(svs,img,systemModel,likelihood)

せっかくなので第一回Tokyo.ScipyでこれについてLTさせてもらいました。
そのときの資料は以下。

Particle Filter Tracking in Python

View more presentations from Kohta Ishikawa

2011-07-05

多変量対数正規分布の計算 (basic statistics of multivariate lognormal distribution)

機械学習・統計

多変量正規分布については色々な解説がありますが、多変量対数正規分布についての結果が知りたいと思ったら意外と(英語でも)解説が見つからなかったので、仕方なく自分で計算しました。メモ代わりにまとめておきます。

Basic statistics(means and (co)variances) of multivariate lognormal distribution are calculated below. Although explanations are in Japanese, you might see mathematical expressions and get information you want. Comments are welcome if you have any questions.

多変量対数正規分布の確率密度関数

多変量対数正規分布の確率密度関数は
$P(\bf{y})=\left(\frac{1}{\sqrt{2\pi}}\right)^p \frac{1}{|\bf{\Sigma}|}\frac{1}{y_1y_2\cdots y_p}e^{-\frac{1}{2}(\log \bf{y}-\bf{\mu})^T\bf{\Sigma}^{-1}(\log \bf{y}-\bf{\mu})}$
という形になります。ベクトル $\bf{\mu}$ や行列 $\bf{\Sigma}$ は多変量正規分布の平均ベクトル、分散共分散行列に対応しますが、多変量対数正規分布の平均や分散共分散にはなっていないことに注意が必要です。これらのパラメータを用いて、多変量対数正規分布の平均ベクトルと分散共分散行列を求めることが目的になります。

変数変換

1変数の対数正規分布は、 $x=(\log y - \mu)/\sigma$ と変数変換をすると $x$ についての標準正規分布になります。これに対応する変数変換のを多変量対数正規分布に対しても行うことで、各種統計量の計算を行います。
今、 $\bf{\Sigma}$ の固有値 $\lambda_k$ と固有ベクトル $\bf{u}_k$ が分かったとします。このとき $\bf{\Sigma}$ は
$\Sigma_{ij} = \sum_k \lambda_k u_{ki}u_{kj}$
および $\bf{\Sigma}^{-1}$ は
$\Sigma^{-1}_{ij} = \sum_k \frac{1}{\lambda_k} u_{ki}u_{kj}$
と書けます（証明略）。ここで、 $u_{ki}$ はk番目の固有ベクトルの第i成分を表します。これを利用して、変数変換を
$t_i = \sum_j \frac{1}{\sqrt{\lambda_i}}u_{ij}(\log y_j - \mu_j)$
とします。多変量の積分の変数変換にはヤコビアンを計算する必要があり、
$\frac{d(t_1, \cdots, t_p)}{d(y_1, \cdots, y_p)} = \left| \begin{matrix} \frac{dt_1}{dy_1} & \cdots & \frac{dt_1}{dy_p} \\ \vdots & & \vdots \\ \frac{dt_p}{dy_1} & \cdots & \frac{dt_p}{dy_p} \\ \end{matrix} \right| = \left| \begin{matrix} \frac{u_{11}}{\sqrt{\lambda_1}y_1} & \cdots & \frac{u_{1p}}{\sqrt{\lambda_1}y_p} \\ \vdots & & \vdots \\ \frac{u_{p1}}{\sqrt{\lambda_p}y_1} & \cdots & \frac{u_{pp}}{\sqrt{\lambda_p}y_p} \\ \end{matrix} \right|$
となります。行列式の性質を使って計算を進めると
$\frac{d(t_1, \cdots, t_p)}{d(y_1, \cdots, y_p)} = \frac{1}{y_1\cdots y_p}\frac{1}{\sqrt{\lambda_1 \cdots \lambda_p}}\left|U|$
となります( $U_{ij} = u_{ij}$ )。Uは直交行列なので、行列式は1となり、結局ヤコビアンは
$\frac{d(t_1, \cdots, t_p)}{d(y_1, \cdots, y_p)} = \frac{1}{y_1\cdots y_p}\frac{1}{\sqrt{\lambda_1 \cdots \lambda_p}}$
となります。yからtへ変換する際には逆変換となり、その場合、ヤコビアンは逆数となります。

平均

$y_i$ の平均値は
$E(y_i)=\int dy_1 \cdots dy_p \left(\frac{1}{\sqrt{2\pi}}\right)^p \frac{1}{|\bf{\Sigma}|}\frac{y_i}{y_1\cdots y_p}e^{-\frac{1}{2}(\log \bf{y}-\bf{\mu})^T\bf{\Sigma}^{-1}(\log \bf{y}-\bf{\mu})}$
となります。先ほどの変数変換を行うと、確率密度関数の部分は無相関な標準正規分布となり、 $y_i$ は
$y_i = \e\left({\mu_i + \sum_j u_{ji}\sqrt{\lambda_j}t_j}\right)$
となります。ここで $U$ は直交行列であり、逆行列が $U^{-1} = U^T$ であるという性質を使っています。従って、
$E(y_i)=\int dt_1 \cdots dt_p \left(\frac{1}{\sqrt{2\pi}}\right)^p \e\left({\mu_i + \sum_j u_{ji}\sqrt{\lambda_j}t_j}\right) e^{-\frac{1}{2}\bf{t}^2$
となります。積分を計算するためには各 $t_j$ について平方完成すればよく、eの肩に乗っている項のうち関連するもののみを見ると、
$-\frac{1}{2}t_j^2 + u_{ji}\sqrt{\lambda_j}t_j = -\frac{1}{2}(t_j - u_{ji}\sqrt{\lambda_j})^2 + \frac{1}{2}u_{ji}^2 \lambda_j$
とすることができます。 $-\frac{1}{2}(t_j - u_{ji}\sqrt{\lambda_j})^2$ の項は $t_j$ で積分すると規格化定数 $\frac{1}{\sqrt{2\pi}}$ とキャンセルして消えるので、結局残るのは
$E(y_i)= \e \left(\mu_i + \frac{1}{2} \sum_j u_{ji}^2 \lambda_j \right)$
になります。ここで、
$\Sigma_{ij} = \sum_k \lambda_k u_{ki}u_{kj}$
という式を思い出すと、
$\frac{1}{2} \sum_j u_{ji}^2 \lambda_j = \frac{1}{2} \Sigma_{ii} = \frac{1}{2} \sigma_i^2$
であることが分かります。結局、多変量対数正規分布の平均値は
$E(y_i)= \e \left(\mu_i +\frac{1}{2} \sigma_i^2)$
と計算できることが分かりました。

分散共分散

平均の次は分散(共分散)を計算します。
共分散の計算はやや面倒ですが、考え方は平均のときと同じものを使うことができます。
$Cov(y_i,y_j)=\int dy_1 \cdots dy_p \left(\frac{1}{\sqrt{2\pi}}\right)^p \frac{1}{|\bf{\Sigma}|}\frac{(y_i-\bar{y}_i)(y_j-\bar{y}_j)}{y_1\cdots y_p}e^{-\frac{1}{2}(\log \bf{y}-\bf{\mu})^T\bf{\Sigma}^{-1}(\log \bf{y}-\bf{\mu})}$
平均値 $\bar{y}_i, \bar{y}_j$ は既に求まっているのでそれを代入し、平均の計算と同様に標準正規分布への変数変換を行うと
$Cov(y_i,y_j)= \int dt_1 \cdots dt_p \left(\frac{1}{\sqrt{2\pi}}\right)^p \left( e^{\mu_i + u_{ki}\sqrt{\lambda_k}t_k} - e^{\mu_i + \frac{1}{2}\sigma_i^2} \right) \left( e^{\mu_j + u_{kj}\sqrt{\lambda_k}t_k} - e^{\mu_j + \frac{1}{2}\sigma_j^2} \right \big) e^{-\frac{1}{2}\bf{t}^2}$
ただし、ここでスペース節約のため、アインシュタインの縮約記法（単一の項内で重複する添字については和を取る）を用いています。すなわち、eの肩に載っている項のkについては和を取るという意味になります。
これらの項のうち、積分後に生き残る（ゼロにならない）のは(eの肩に)tの1乗が含まれる項です。また生き残る項についての積分の結果は、平均の計算のときと同様に、元々の確率密度関数から来る $t^2$ 項と合わせた平方完成による「おつり」の項です。この計算を地道に行うと
$Cov(y_i,y_j) = e^{\mu_i + \mu_j} \left( e^{\frac{1}{2}\lambda_k u_{ki}^2 + \frac{1}{2}\lambda_k u_{kj}^2 + \lambda_k u_{ki}u_{kj}} - e^{\frac{1}{2}(\sigma_i^2 + \sigma_j^2)}\right)$
となり、平均値のときと同様、 $\lambda_k u_{ki}^2 = \sigma_i^2$ 、 $\lambda_k u_{ki}u_{kj} = \sigma_{ij}$ となるので、結局
$Cov(y_i,y_j) = e^{\mu_i + \frac{1}{2}\sigma_i^2 + \mu_j + \frac{1}{2}\sigma_j^2}\left( e^{\sigma_{ij}} - 1 \right) = E(y_i)E(y_j) \big( e^{\sigma_{ij}} - 1\big)$
となります。

2011-05-28

Tokyo.R #14で発表しました

R C++ 機械学習・統計

気がついたら結構参加回数を重ねていたTokyo.Rにて、プレッシャーに耐え切れなくなり発表することになりました。
isomapという多様体学習の手法が使いたくて実装したので、それの周りについての発表としました。先日の多次元尺度構成法もこの手法の一部として使われています。
ただ実装についてはまだ少々バグがあるっぽいのと、isomapだけならveganというパッケージで使えることがわかったため、公開はもう少しブラッシュアップしてからにしようかなと思っています。
少しパンチが弱いなあと思っていますが、思ったより好評だったようでよかったです。

Rでisomap（多様体学習のはなし）

View more presentations from Kohta Ishikawa

ちなみに、slideshareには別の資料も少し置いてあって、統計学入門は「統計とか興味あるけどよくわからないなあ」という人におすすめです。

2011-05-14

多次元尺度構成法イントロダクション

R 機械学習・統計

ちょっと多次元尺度構成法について理解する必要があったので簡単なまとめです。

多次元尺度構成法とは、多数の多次元点間の距離データのみが与えられたときに、その距離を再現するような空間(座標系)を逆算する手法です。(多次元)座標値が与えられれば距離は自由に計算できますが、距離だけが与えられた場合に座標値を計算するのは直感的には中々できません。

今、n個の点がありそれぞれの点間の $n^2-n$ 個の距離 $d_{ij} \quad i,j\in \{1,\cdots,n\}$ がデータとして与えられたとします( $d_{ii}$ はゼロです)。空間の次元をqとし、未知のn個の座標ベクトルを $\bf{x}_i \in \mathcal{R}^q, \quad i \in {1,\cdots,n}$ 、それらを縦に並べた行列を

$\bf{X} = \left( \begin{array} \bf{x}_1^t \\ \vdots \\ \bf{x}_n^t \end{array}\right)$

とします(tは転地記号)。また、点間の距離はユークリッド距離

$d_{ij}^2 = (\bf{x_i} - \bf{x_j})^2 = \sum_k \left(x_{ik}x_{ik} + x_{jk}x_{jk} - 2*x_{ik}x_{jk}\right)$

とします。
ここで

$\bf{B} = \bf{X}\bf{X}^t$

という行列を考えます。もし何らかの形で距離データからBを求めることができれば、Bを対角化し固有値の平方根をとることで、 $\bf{X}$ を求めることができます。 $\bf{B}$ の成分は

$b_{ij} = \sum_k x_{ik}x^t_{kj} = \sum_kx_{ik}x_{jk}$ 　　　　　　　　(a)

であり、距離行列の成分と比べて次のような関係にあることがわかります。

$d_{ij}^2 = b_{ii} + b_{jj} - 2b_{ij}$ 　　　　　　　　(b)

(少々ややこしいですが、ここではi,jについては和を取りません。左辺にも現れるもので、ダミー変数ではないとわかります。)
この関係式によって、距離行列(手元にあるデータ)と $\bf{B}$ を結びつけることができます。

ところが、実はこの関係式だけでは距離行列から $\bf{B}$ を(すなわち座標値を)決めることができません。考えてみれば当然で、座標軸は原点をずらしたり回転したりしても距離は変わりません。逆に、距離だけがわかっている場合は原点や回転の自由度があるため座標軸を一意に決めることができないのです。
回転については直交変換の自由度が残ることになりますが、原点の位置を決めればまずは一つの座標系を書き下すことができるようになります。普通、原点は計算したデータ座標値の平均がゼロとなるようにとられます。

$\sum_{i=1}^n x_{ik} = 0, \quad k\in \{1,\cdots,q\}$

この関係式を(a)式に適用すると、 $b_{ij}$ もまたiまたはjの一方について和をとるとゼロとなることがわかります。これにより、(b)についてiで和をとると

$\sum_i d_{ij}^2 = \text{Tr}(\bf{B}) + nb_{jj}$

となります。ここで $\text{Tr}(\bf{B})$ は $\bf{B}$ の対角成分の和(トレース)で、nは $\bf{B}$ の次元です。同様に、jについての和、i,j両方についての和をとると

$\sum_j d_{ij}^2 = nb_{ii} + \text{Tr}(\bf{B})$

$\sum_{i,j} d_{ij}^2 = 2n\text{Tr}(\bf{B})$

という関係式も得られ、これら3つを $b_{ii},b_{jj},\text{Tr}(\bf{B})$ について解くことで、(b)式より $b_{ij}$ を $d_{ij}^2$ を用いて書き下すことができます。

$b_{ij} = -\frac{1}{2}\big( d_{ij}^2 - d_{i\cdot}^2 - d_{\cdot j}^2 + d_{\cdot \cdot}^2 \big)$
$d_{i\cdot}^2 = \frac{1}{n}\sum_j d_{ij}^2$
$d_{\cdot j}^2 = \frac{1}{n}\sum_i d_{ij}^2$
$d_{\cdot \cdot}^2 = \frac{1}{n^2}\sum_{ij} d_{ij}^2$
これで $\bf{B}$ がわかったので、ここから $\bf{B}=\bf{X}\bf{X}^t$ を解いて $\bf{X}$ を求めます。そのためには $\bf{X}$ の固有値・固有ベクトルを求めればよいというのが線形代数の常套手段です。実際

$\bf{B}\bf{x}^k = \lambda_k \bf{x}^k$

となる固有値 $\lambda_k$ 、固有ベクトル $\bf{x}^k$ が求まったとすると、

$\bf{B} = \sum_k\lambda_k\bf{x}^k\bf{x}^{kt}$

と書くことができ、これは、固有ベクトルを並べた行列

$\bf{V} = (\bf{x}_1, \cdots , \bf{x}_n)$

と、固有値を対角成分に並べた行列 $\bf{\Lambda} = \text{diag} \{\lambda_1, \cdots , \lambda_n\}$ を用いて、

$\bf{B} = \bf{V}\bf{\Lambda}\bf{V}^t = \bf{V}\bf{\Lambda}^{1/2}\left(\bf{V} \bf{\Lambda}^{1/2}\right)^t$

と書けます。ただし $\bf{\Lambda}^{1/2} = \text{diag} \{\sqrt{\lambda_1}, \cdots , \sqrt{\lambda_n}\}$ です。
従って、求める $\bf{X}$ は

$\bf{X} = \bf{V}\bf{\Lambda}^{1/2}$

であることがわかります。
ここで注意が必要です。固有値の平方根を要素とする行列を作りましたが、固有値は一般に負の数や複素数に成り得るので、この議論が成り立つためには固有値は0または正でなければなりません。これは行列 $\bf{B}$ が半正定値でなければならないということを意味します。実際には必ずしも半正定値ではなく、また欲しいのはなるべくズレの少ない低次元表現ですので、固有値を降順にソートし、大きいものからp個持ってきた行列で $\bf{\Lambda}^{1/2}$ を近似します。これによってp次元で表現された座標値 $\bf{X}$ が求まります。

長々と解説しましたが、Rでこれを実行するのは簡単です。statsパッケージにcmdscaleという関数があり、デフォルトでは2次元の多次元尺度構成法を行ってくれます。

> library(stats)
> iris.dist <- dist(iris[,-5])        #距離行列を作成
> iris.cmd <- cmdscale(iris.dist,k=2) #多次元尺度構成法
> plot(iris.cmd,col=iris$Species)     #結果を色つきでプロット

以下のようなグラフが得られます。

irisは元々4次元のデータですが、多次元尺度構成法によって得られた2次元表現でも、そこそこきちんと種別を反映してくれていることがわかると思います。irisは元データが座標値に対応するため、低次元表現を求めるという意味では主成分分析と同等ですが、データによってはそもそも距離データしか得られないということもあり、そういう場合には得にデータの可視化に役立ちます。

今回紹介した多次元尺度構成法は古典的多次元尺度構成法と呼ばれるもので、実際には距離がユークリッドでない場合や非線形な場合などを主眼に置いた様々な改良手法が存在します。いずれその辺りも書くかもしれません。

参考文献

RとS-PLUSによる多変量解析

作者: B.エヴェリット,石田基広
出版社/メーカー: シュプリンガー・ジャパン株式会社
発売日: 2007/06/28
メディア: 単行本
クリック: 40回
この商品を含むブログ (9件) を見る

Rと多次元尺度法

2011-04-28

Rcppを使うときのメモ

R C++

Rcppを使おうとして苦戦したのでメモ

Rcppのインストール

R上で

install.packages("Rcpp")

すればOK

C++側から読み込むヘッダファイル(Rcpp.h)の設定

C:\Program Files\R\R-2.13.0\include

などにあるので、パスを張るなりg++等のインクルードディレクトリに追加するなりする必要がある。

コンパイル環境の整備

RからGSL(WinXP用)
この辺に書いてあるようなMinGW,MSYS,(Perl?)のインストールとパス設定が必要(GSLは必須ではない)。

C++ソースコードのコンパイル

R用関数を書いたソースをコンパイルする際、色々なサイトに書いてあるのは

R CMD SHLIB hogehoge.cpp

というコマンド。これは単純にg++のコマンドをラップしたものを実行するよう。ところが、Makevarsファイルで

PKG_CPPFLAGS = -I"C:\Program Files\hoge"

とか書いても、ダブルクォーテーションが勝手に省略されてしまうよう。この場合は仕方が無いのでR CMD SHLIB hogehoge.cppでコンソールに出たコマンドを適宜修正して直接打ち込むことでコンパイルできる。

ググって出てくる内容にこのあたりを追加すれば、大体使えるのではないかと思います。