2012-12-11

Fedora17でAmazon Elastic MapReduce Ruby Clientを使った際のメモ

fedora

Amazon Elastic MapReduce Ruby Clientがruby 1.8用であるため、Fedora17では使えない。
ruby1.8をインストールしてAmazon Elastic MapReduce Ruby ClientをFedora17で使えるようにした。
なお、https://github.com/rbotzer/elastic-mapreduce-rubyにruby1.9用のクライアントがあるが、
このクライアントを使ってもEMRを使う事が出来なかった。

ruby1.8のインストール

Fedora17でruby 1.8.7のインストール作業メモ

1. ruby 1.8.7のソースをダウンロードする

2. パッチを当てる

ruby1.8.7はglibc2.14が原因でコンパイルが通らないとのこと。
http://sourceforge.jp/forum/message.php?msg_id=59221

以下のパッチをダウンロードして、
http://bugs.ruby-lang.org/attachments/1931/stdout-rouge-fix.patch

以下の手順でパッチを適用する。

mkdir ruby
mv ruby-1.8.7-p371.zip ruby
mv stdout-rouge-fix.patch ruby
cd ruby
unzip ruby-1.8.7-p371.zip
cd ruby-1.8.7-p371
patch -p1 < ../stdout-rouge-fix.patch

3. makeしてインストールする

今回は/usr/local/ruby/1.8.7にインストールする。
以下の手順でコンパイルとインストールを実行する。

./configure --prefix=/usr/local/ruby/1.8.7
make
sudo make install

4. rubyをruby1.8.7に切り替える

以下のコマンドで、ruby1.8.7に切り替える。

sudo alternatives --install /usr/bin/ruby ruby /usr/local/ruby/1.8.7/bin/ruby 30

環境変数PATHにruby1.8.7へのパスを追加する。
~/.bashrcに以下の行を追加する。

export PATH=/usr/local/ruby/1.8.7/bin:$PATH

5. Amazon Elastic MapReduce Ruby Clientをインストールする
後はマニュアル通りにインストールして使用すればよい。

つまづいた点

credentials.jsonを作成して実行し、ログを見ると以下のようなエラーが出た。

the given ssh key name was invalid

原因は、EC2のキーペアをUS EASTで作成していたため。
TOKYOのEC2でキーペアを作りなおして、TOKYOのEMRを使ったら動いた。実に間抜け。

2012-12-01

Fedora17でCDH4をインストールし、擬似分散モードを起動させた際のメモ

Hadoop

Fedora17でCDH4をインストールした。
また擬似分散モードでの起動を確認した。

インストール

JDK6をインストールする

CDH4がOracleのJDK6に依存しているので、
http://www.oracle.com/technetwork/java/javase/downloads/jdk6u37-downloads-1859587.htmlよりインストーラをダウンロードしてインストールする。

$ cd ダウンロード
$ chmod +x jdk-6u37-linux-x64-rpm.bin
$ sudo ./jdk-6u37-linux-x64-rpm.bin 
$ sudo alternatives --install /usr/bin/java java /usr/java/jdk1.6.0_37/bin/java 20000

~/.bashrcに以下を追加する

$ export JAVA_HOME=/usr/java/jdk1.6.0_37
$ export PATH=$JAVA_HOME/bin:$PATH

以下で有効にする。

$ source ~/.bashrc

javaのバージョンを確認する。

$ java -version
java version "1.6.0_37"
Java(TM) SE Runtime Environment (build 1.6.0_37-b06)
Java HotSpot(TM) 64-Bit Server VM (build 20.12-b01, mixed mode)
$javac -version
javac 1.6.0_37

CDH4.1.2をインストールする

CDHをyumを使ってインストールする。

CDH用にリポジトリを追加する。

https://ccp.cloudera.com/display/CDH4DOC/CDH4+Installation#CDH4Installation-AddingaYumRepositoryのRedHat/CentOS6 64-bit向けrepoファイルをコピーする。

cloudera-cdh4.repo

[cloudera-cdh4]
name=Cloudera's Distribution for Hadoop, Version 4
baseurl=http://archive.cloudera.com/cdh4/redhat/6/x86_64/cdh/4/
gpgkey = http://archive.cloudera.com/cdh4/redhat/6/x86_64/cdh/RPM-GPG-KEY-cloudera    
gpgcheck = 1

上記のテキストファイルcloudera-cdh4.repoを/etc/yum.repos.d以下にコピーする。

リポジトリキーを追加する。

$ sudo rpm --import http://archive.cloudera.com/cdh4/redhat/6/x86_64/cdh/RPM-GPG-KEY-cloudera

yumからHadoopをインストールする。

$ sudo yum clean all
$ sudo yum install hadoop

完全分散モード用の設定ファイルをコピーして、オリジナルをいじらないようにしておく。
以下のコマンドを実行すると、完全分散モード用の設定ファイルがconf.myclusterにコピーされ、
/etc/hadoop/confからconf.myclusterのファイルを修正出来るよになる。

$ cd /etc/hadoop
$ sudo cp -rf conf.empty conf.mycluster

Hadoopを動かしてみる

スタンドアローンモードで動かす

スタンドアローンモードは特に設定ファイルは不要。
必要なパッケージをインストールして、hadoopのgrepを使ってみる。

$ sudo yum install hadoop-mapreduce hadoop-0.20-mapreduce

作業用ディレクトリ~/hadoop/test0/を作成して、そこに適当なテキストファイルを放り込み、hadoopのgrepを実行する。

$ mkdir -p ~/hadoop/test0/
$ cd ~/hadoop/test0/
$ mkdir input
$ cp ~/sample.txt input
$ hadoop jar /usr/lib/hadoop-mapreduce/hadoop-*-examples.jar grep ./input/ output 'hoge'

このコマンドを実行すると、outputディレクトリが作成されて、'hoge'の単語数が./output/part-00000に出力される。

擬似分散モードで動かす

擬似分散モード用のパッケージをインストールして、
擬似分散モードでhadoopのgrepを使ってみる。

$ sudo yum install hadoop-conf-pseudo

このコマンドを実行すると、alternativesで/etc/hadoop/confに/etc/hadoop/conf.pseudoがインストールされる。

$ alternatives --display hadoop-conf
hadoop-conf -ステータスは自動です。
リンクは現在 /etc/hadoop/conf.mycluster を指しています。
/etc/hadoop/conf.empty - 優先項目 10
/etc/hadoop/conf.mycluster - 優先項目 60
/etc/hadoop/conf.pseudo - 優先項目 30
現在の「最適」バージョンは /etc/hadoop/conf.mycluster です。

設定ファイルを擬似分散モードに切り替える。

$ sudo alternatives --config hadoop-conf 

3 プログラムがあり 'hadoop-conf' を提供します。

  選択       コマンド
-----------------------------------------------
   1           /etc/hadoop/conf.empty
*+ 2           /etc/hadoop/conf.mycluster
   3           /etc/hadoop/conf.pseudo

Enter を押して現在の選択 [+] を保持するか、選択番号を入力します:3

*** HDFSのフォーマットと起動

formatしてRe-format filesystem in Storage Directory?でYを選択するclusterIDが更新される。

$ sudo -u hdfs hdfs namenode -format
$ sudo service hadoop-hdfs-namenode start
$ sudo service hadoop-hdfs-datanode start
$ sudo service hadoop-hdfs-secondarynamenode start

hdfs中に必要なディレクトリを作成をする。

$ sudo -u hdfs hadoop fs -mkdir /tmp
$ sudo -u hdfs hadoop fs -chmod -R 1777 /tmp
$ sudo -u hdfs hadoop fs -ls /
Found 1 items
drwxrwxrwt   - hdfs supergroup          0 2012-11-30 21:35 /tmp
$ sudo -u hdfs hadoop fs -mkdir /var/log/hadoop-yarn
$ sudo -u hdfs hadoop fs -chown yarn:mapred /var/log/hadoop-yarn
$ sudo -u hdfs hadoop fs -mkdir /tmp/hadoop-yarn/staging
$ sudo -u hdfs hadoop fs -chmod -R 1777 /tmp/hadoop-yarn/staging
$ sudo -u hdfs hadoop fs -mkdir /tmp/hadoop-yarn/staging/history/done_intermediate
$ sudo -u hdfs hadoop fs -chmod -R 1777 /tmp/hadoop-yarn/staging/history/done_intermediate
$ sudo -u hdfs hadoop fs -chown -R mapred:mapred /tmp/hadoop-yarn/staging
$ sudo -u hdfs hadoop fs -mkdir /user/$USER
$ sudo -u hdfs hadoop fs -chown $USER /user/$USER

リソースマネージャの起動をする。

$ sudo service hadoop-yarn-resourcemanager start
$ sudo service hadoop-yarn-nodemanager start
$ sudo service hadoop-mapreduce-historyserver start

ここまで実行すると、psコマンドで以下のプロセスが立ち上がっている事を確認する事が出来る。

$  ps -ef | grep java
hdfs      2113     1  0 18:06 ?        00:00:10 /usr/java/jdk1.6.0_37/bin/java -Dproc_namenode -Xmx1000m -Dhadoop.log.dir=/var/log/hadoop-hdfs -Dhadoop.log.file=hadoop-hdfs-namenode-localhost.localdomain.log -Dhadoop.home.dir=/usr/lib/hadoop -Dhadoop.id.str=hdfs -Dhadoop.root.logger=INFO,RFA -Djava.library.path=/usr/lib/hadoop/lib/native -Dhadoop.policy.file=hadoop-policy.xml -Djava.net.preferIPv4Stack=true -Dhadoop.security.logger=INFO,RFAS org.apache.hadoop.hdfs.server.namenode.NameNode
hdfs      2248     1  0 18:07 ?        00:00:09 /usr/java/jdk1.6.0_37/bin/java -Dproc_datanode -Xmx1000m -Dhadoop.log.dir=/var/log/hadoop-hdfs -Dhadoop.log.file=hadoop-hdfs-datanode-localhost.localdomain.log -Dhadoop.home.dir=/usr/lib/hadoop -Dhadoop.id.str=hdfs -Dhadoop.root.logger=INFO,RFA -Djava.library.path=/usr/lib/hadoop/lib/native -Dhadoop.policy.file=hadoop-policy.xml -Djava.net.preferIPv4Stack=true -server -Dhadoop.security.logger=INFO,RFAS org.apache.hadoop.hdfs.server.datanode.DataNode
hdfs      2446     1  0 18:09 ?        00:00:05 /usr/java/jdk1.6.0_37/bin/java -Dproc_secondarynamenode -Xmx1000m -Dhadoop.log.dir=/var/log/hadoop-hdfs -Dhadoop.log.file=hadoop-hdfs-secondarynamenode-localhost.localdomain.log -Dhadoop.home.dir=/usr/lib/hadoop -Dhadoop.id.str=hdfs -Dhadoop.root.logger=INFO,RFA -Djava.library.path=/usr/lib/hadoop/lib/native -Dhadoop.policy.file=hadoop-policy.xml -Djava.net.preferIPv4Stack=true -Dhadoop.security.logger=INFO,RFAS org.apache.hadoop.hdfs.server.namenode.SecondaryNameNode
yarn      2952     1  0 18:12 ?        00:00:10 /usr/java/jdk1.6.0_37/bin/java -Dproc_resourcemanager -Xmx1000m -Dhadoop.log.dir=/var/log/hadoop-yarn -Dyarn.log.dir=/var/log/hadoop-yarn -Dhadoop.log.file=yarn-yarn-resourcemanager-localhost.localdomain.log -Dyarn.log.file=yarn-yarn-resourcemanager-localhost.localdomain.log -Dyarn.home.dir=/usr/lib/hadoop-yarn -Dhadoop.home.dir=/usr/lib/hadoop-yarn -Dhadoop.root.logger=INFO,RFA -Dyarn.root.logger=INFO,RFA -Djava.library.path=/usr/lib/hadoop/lib/native -classpath /etc/hadoop/conf:/etc/hadoop/conf:/etc/hadoop/conf:/usr/lib/hadoop/lib/*:/usr/lib/hadoop/.//*:/usr/lib/hadoop-hdfs/./:/usr/lib/hadoop-hdfs/lib/*:/usr/lib/hadoop-hdfs/.//*:/usr/lib/hadoop-yarn/lib/*:/usr/lib/hadoop-yarn/.//*:/usr/lib/hadoop-mapreduce/lib/*:/usr/lib/hadoop-mapreduce/.//*:/usr/lib/hadoop-yarn/.//*:/usr/lib/hadoop-yarn/lib/*:/etc/hadoop/conf/rm-config/log4j.properties org.apache.hadoop.yarn.server.resourcemanager.ResourceManager
yarn      3224     1  1 18:12 ?        00:00:13 /usr/java/jdk1.6.0_37/bin/java -Dproc_nodemanager -Xmx1000m -server -Dhadoop.log.dir=/var/log/hadoop-yarn -Dyarn.log.dir=/var/log/hadoop-yarn -Dhadoop.log.file=yarn-yarn-nodemanager-localhost.localdomain.log -Dyarn.log.file=yarn-yarn-nodemanager-localhost.localdomain.log -Dyarn.home.dir=/usr/lib/hadoop-yarn -Dhadoop.home.dir=/usr/lib/hadoop-yarn -Dhadoop.root.logger=INFO,RFA -Dyarn.root.logger=INFO,RFA -Djava.library.path=/usr/lib/hadoop/lib/native -classpath /etc/hadoop/conf:/etc/hadoop/conf:/etc/hadoop/conf:/usr/lib/hadoop/lib/*:/usr/lib/hadoop/.//*:/usr/lib/hadoop-hdfs/./:/usr/lib/hadoop-hdfs/lib/*:/usr/lib/hadoop-hdfs/.//*:/usr/lib/hadoop-yarn/lib/*:/usr/lib/hadoop-yarn/.//*:/usr/lib/hadoop-mapreduce/lib/*:/usr/lib/hadoop-mapreduce/.//*:/usr/lib/hadoop-yarn/.//*:/usr/lib/hadoop-yarn/lib/*:/etc/hadoop/conf/nm-config/log4j.properties org.apache.hadoop.yarn.server.nodemanager.NodeManager
mapred    3558     1  0 18:17 ?        00:00:06 /usr/java/jdk1.6.0_37/bin/java -Dproc_historyserver -Xmx1000m -Dhadoop.log.dir=/var/log/hadoop-mapreduce -Dhadoop.log.file=yarn-mapred-historyserver-localhost.localdomain.log -Dhadoop.home.dir=/usr/lib/hadoop -Dhadoop.id.str= -Dhadoop.root.logger=INFO,RFA -Djava.library.path=/usr/lib/hadoop/lib/native -Dhadoop.policy.file=hadoop-policy.xml -Djava.net.preferIPv4Stack=true -Dmapred.jobsummary.logger=INFO,JSA -Dhadoop.security.logger=INFO,NullAppender org.apache.hadoop.mapreduce.v2.hs.JobHistoryServer

jpsコマンドを使うと以下のように３つだけ表示される。

$ sudo jps
3558 JobHistoryServer
4867 Jps
2952 ResourceManager
3224 NodeManager

以上で擬似分散モードでのHDFSの起動が完了する。

つまずいた箇所

DataNodeが起動しない

"sudo service hadoop-hdfs-namenode"を実行しても、DataNodeが起動しなかった。
原因は、以下のDataNodeのログの一部からもわかる通り、NameNodeとDataNodeでclusterIDが一致しなかったためだ。

DataNodeの実行ログ/var/log/hadoop-hdfs/hadoop-hdfs-datanode-localhost.localdomain.log の一部抜粋

2012-12-01 17:47:56,233 FATAL org.apache.hadoop.hdfs.server.datanode.DataNode: Initialization failed for block pool Block pool BP-1334899232-127.0.0.1-1354351535263 (storage id DS-776307469-127.0.0.1-50010-1354277718964) service to localhost/127.0.0.1:8020
java.io.IOException: Incompatible clusterIDs in /var/lib/hadoop-hdfs/cache/hdfs/dfs/data: namenode clusterID = CID-51ca4356-b4a2-4646-bea3-9af43adc50a1; datanode clusterID = CID-fa5c4789-7c02-4884-91fd-921d46bbfa7a
        at org.apache.hadoop.hdfs.server.datanode.DataStorage.doTransition(DataStorage.java:390)
        at org.apache.hadoop.hdfs.server.datanode.DataStorage.recoverTransitionRead(DataStorage.java:190)
        at org.apache.hadoop.hdfs.server.datanode.DataStorage.recoverTransitionRead(DataStorage.java:218)
        at org.apache.hadoop.hdfs.server.datanode.DataNode.initStorage(DataNode.java:849)
        at org.apache.hadoop.hdfs.server.datanode.DataNode.initBlockPool(DataNode.java:820)
        at org.apache.hadoop.hdfs.server.datanode.BPOfferService.verifyAndSetNamespaceInfo(BPOfferService.java:308)
        at org.apache.hadoop.hdfs.server.datanode.BPServiceActor.connectToNNAndHandshake(BPServiceActor.java:218)
        at org.apache.hadoop.hdfs.server.datanode.BPServiceActor.run(BPServiceActor.java:661)
        at java.lang.Thread.run(Thread.java:662)

clusterIDを一致させる正当な手順が見つけられなかったため、今回は以下のファイルのclusterIDを直接エディタで編集して問題を解決した。

以下のファイルの"clusterID="以降をエディタで編集した。
/var/lib/hadoop-hdfs/cache/hdfs/dfs/data/current/VERSION
/var/lib/hadoop-hdfs/cache/hdfs/dfs/namesecondary/current/VERSION

NameNodeのclusterIDは以下のファイルに記述されている。
/var/lib/hadoop-hdfs/cache/hdfs/dfs/name/current/VERSION

これらファイルを編集して、擬似分散モードの起動手順を繰り返したところ、全てのサービスが立ち上がり、hadoopのサンプルプログラムが実行出来た。

2012-11-18

パラメータの事後分布計算を実装した

Haskell 機械学習 PRML

PRML3.3.1のパラメータの分布を読んで、パラメータの事後分布計算を実装した。

事前分布、事後分布、尤度関数の関係

観測データ $\bf{D}$ が得られた時に、そのデータがパラメータ $\bf{w}$ に依存していたとする。
この時、観測データ $\bf{D}$ があるパラメータ $\bf{w}$ の条件下で得られる確率は $P(\bf{D}|\bf{w})$ である。
$P(\bf{D}|\bf{w})$ は $\bf{w}$ の値によって変わる関数なので、尤度関数と呼ばれる。
また、パラメータ $\bf{w}$ の確率分布 $P(\bf{w})$ が事前分布であり、
観測データ $\bf{D}$ が得られた後でわかるパラメータ $\bf{w}$ の確からしさは事後分布 $P(\bf{w}|\bf{D})$ である。
これら事前分布、事後分布、尤度関数は、ベイズの定理より以下の関係にある。(PRML 1.43式参照)
$P(\bf{w}|\bf{D}) = \frac {P(\bf{D}|\bf{w}) P(\bf{w})} {p(\bf{D})}$
つまり、事後分布 $P(\bf{w}|\bf{D})$ $\propto$ 尤度関数 $P(\bf{D}|\bf{w})$ $\times$ 事前分布 $P(\bf{w})$ である。

尤度関数

観測データ $\bf{D}$ のそれぞれの値が独立同分布である時、観測データ $\bf{D}$ の得られる確率はは各データが得られる確率の積となる。
例えば、独立同分布で得られる値 $y_1, y_2$ の同時確率は $P(y_1,y_2) = P(y_1)P(y_2)$ である。
同様に、パラメータ $\bf{w}$ の条件下で観測データ $\bf{D}$ の得られる確率 $P(\bf{D}|\bf{w})$ は、
それぞれの観測値が得られる確率の積となるので以下のように表現できる。
$P(\bf{D}|\bf{w}) = \prod_{i=1}^{N} P(t_i|x_i,\bf{w})$

特に、PRML3.1で導入した線形基底関数モデルを採用した場合に、入力 $\bf{X} = \{x_1, ..., x_N\}$ と対応する目標値 $\bf{t} = {t_1, ... , t_N}$ が得られた場合の尤度関数は、
$P(\bf{t}|\bf{X},\bf{w},\beta) = \prod_{n=1}^{N}N(t_n|\bf{w}^{T}\bf{\phi}(x_n),\beta^{-1}) ... (3.10)$
なお、入力変数のモデル化はしないので、以降 $\bf{X}$ は省略する。

パラメータの分布

式(3.10)の尤度関数 $P(\bf{t}|\bf{w})$ の指数部分が $\bf{w}$ の2次形式であるため、
この尤度関数に対応する共役事前分布が以下のようなガウス分布で与えられる事に注目する。
$P(\bf{w}) = N(\bf{w}|\bf{m_0}, \bf{S_0})$
事後分布は尤度関数と事前分布の積に比例するので、以下のように書ける。
$P(\bf{w}|\bf{t}) \propto \exp \Bigl(-\frac {1}{2} (\bf{w} - \bf{m_0})^{T}\bf{S}_{0}^{-1}(\bf{w} - \bf{m_0}) - \sum_{n=1}^{N}(t_{n}-\bf{w}^{T}\bf{\phi}(x_n))^2 \Bigr)$
この式の指数部分を平方完成して、観測データ取得後の事後分布 $P(\bf{w}|\bf{t}) = N(\bf{w}|\bf{m_N},\bf{S_N})$ も指数部分を展開して平方完成して比較すると、
以下のような期待値 $\bf{m_N}$ と共分散 $\bf{S_{N}^{-1}}$ が得られる。
$\bf{m_N} = \bf{S_N}(S_{0}^{-1}\bf{m_{0}} + \beta \bf{\Phi}^{T}\bf{t})$
$\bf{S_{N}^{-1}} = \bf{S_{0}^{-1}} + \beta \bf{\Phi}^{T} \bf{\Phi}$

ここで、 $P(\bf{w}|\alpha) = N(\bf{w}|\bf{0}, \alpha^{-1}\bf{I})$ を事前分布とした場合に、
目標値 $\bf{t}$ と入力値 $\bf{X}$ が与えられた場合の、事後分布 $P(\bf{w}|\bf{t})$ を計算した結果を示す。
目標値 $\bf{t}$ は、 $t = 3.14 * x + 1.23$ で計算したtに-1から1の範囲で得られるランダムな値を加えて作成した。
以下のグラフは、横軸に $w_0$ 、縦軸に[w_1]を取り、色で座標上の確率を表したものである。

$w_0 = 1.23, w_1 = 3.14$ の座標で確立が集中的に高くなっていることがわかる。

この計算を行うプログラムを以下に示す。

import System.Random
import Numeric.LinearAlgebra
import Graphics.Gnuplot.Simple
import Control.Applicative

-- ノイズを付与する
addNoise :: (RandomGen g, Random a, Num a) =>
            (a, a) -> [a] -> g -> [a]
addNoise range xs g = zipWith (+) xs (randomRs range g)

-- 座標を文字列にする
from2dToString :: (Show a) => [a] -> [a] -> String
from2dToString xs ys =
  concat $ zipWith (\x y -> show x ++ " " ++ show y ++ "\n") xs ys

-- 座標を文字列にする
from3dToString :: (Show a1, Show a2, Show a3)
                  => [(a1, a2)] -> [a3] -> String
from3dToString xs ts = concat $
                      zipWith (\(x,y) z ->    (show x) ++ " "
                                           ++ (show y) ++ " "
                                           ++ (show z) ++ "\n")
                      xs ts

-- 多変量ガウス分布
-- mu    : D次元平均ベクトル 実装上はD×1の行列
-- sigma : D×Dの共分散行列
-- x     : D次元ベクトルx   実装上はD×1の行列
-- gaussDist :: Matrix a -> Matrix a -> Matrix a -> a
gaussDist mu sigma x = c * exp a
  where
    d = fromIntegral $ cols x
    c = 1.0 / (((2 * pi) ** (d / 0.5)) * ((det sigma) ** 0.5))
    y = x - mu
    a = -0.5 * mahalanobisDistance mu sigma x
    mahalanobisDistance avrg covariance x
      = ((trans $ x - avrg)
         `multiply` (inv covariance)
         `multiply` y) @@> (0,0)
      where
        y = x - avrg

-- 多項式基底関数
polynomialBasis :: Floating a => a -> a -> a
polynomialBasis n x = x ** n

-- 基底関数のリストと入力変数xからベクトル?(x)を計算する。
-- 実装上は M×1の行列
basisFuncVector basiss x = (m >< 1) $ [f x | f <- basiss]
  where
    m = length basiss

-- 基底関数の配列と入力変数xから計画行列Φを計算する
designMatrix basiss xs = (n >< m) [f x | x <- xs, f <- basiss]
  where
    n = length xs
    m = length basiss

-- 事後分布
postDist alpha beta xs ts = gaussDist m s
  where
    phi = designMatrix (map polynomialBasis [0,1]) xs
    phiTphi = (trans phi) `multiply` phi
    n = (cols phiTphi)
    alphaI = diag . fromList $ take n $ repeat alpha
    betaI  = diag . fromList $ take n $ repeat beta
    invS = alphaI + betaI `multiply` phiTphi
    s = inv invS
    m = betaI `multiply` s `multiply` (trans phi)
        `multiply` ((length ts >< 1) ts)

-- 事後分布を受け取り、関数内で生成した各座標での確率を計算する。
-- 座標と確率を文字列にして返す。
calcDist dist = str
  where
    xy  = [(x,y) :: (Double, Double)
           | x <- linearScale 100 (-5, 5)
           , y <- linearScale 100 (-5, 5)]
    xy' = map ((2 >< 1) . (\(x,y) -> [x,y])) xy
    z   = map dist xy'
    str = from3dToString xy z

-- 事後分布を計算してファイルに出力する。
-- xの個数が少ないと1を超えた確率が出る場合がある
runPostDistCalc = do
  gen <- getStdGen
  let xs = linearScale 1000 (0, 10)
      f  = \x -> 3.14 * x + 1.23
      ts = addNoise (-1.0 ::Double, 1.0::Double) (map f xs) gen
      alpha = 1.0
      beta  = 1.0
      g = postDist alpha beta xs ts
      str = calcDist g
  writeFile "gaussDist.dat" str
  return ()

目標値tと入力値xを一つずつ入力して事後分布を更新していく計算方法を考える事も出来る。
以下のアニメーションでは、更新された事後分布 $P(\bf{w})$ が徐々に $w_0 = 1.23, w_1 = 3.14$ に収束していく事がわかる。

このアニメーションの確率分布は以下のプログラムで計算した。

import System.Random
import Numeric.LinearAlgebra
import Graphics.Gnuplot.Simple
import Control.Applicative
import Control.Monad

-- ガウス分布の共分散と期待値を保持するデータ構造
data GaussDist = GaussDist {
    mu    :: Matrix Double
  , sigma :: Matrix Double
} deriving (Show, Eq)

-- パラメータWとその確率を保持するためのデータ構造
data Pos = Pos Double Double Double deriving Eq 
instance Show Pos where
  show (Pos x y z) = show x ++ " " ++
                     show y ++ " " ++
                     show z ++ "\n"

-- 多変量ガウス分布
-- mu    : D次元平均ベクトル 実装上はD×1の行列
-- sigma : D×Dの共分散行列
-- x     : D次元ベクトルx   実装上はD×1の行列
-- gaussDist :: Matrix a -> Matrix a -> Matrix a -> a
gaussDist mu sigma x = c * exp a
  where
    d = fromIntegral $ cols x
    c = 1.0 / (((2 * pi) ** (d / 0.5)) * ((det sigma) ** 0.5))
    y = x - mu
    a = -0.5 * mahalanobisDistance mu sigma x
    mahalanobisDistance avrg covariance x
      = ((trans $ x - avrg)
         `multiply` (inv covariance)
         `multiply` y) @@> (0,0)
      where
        y = x - avrg

-- 座標wの確率を計算する
-- wはパラメータベクトル 実装上はD×1の行列
probability :: GaussDist -> Matrix Double -> Double
probability dist w = gaussDist m s w
  where
    m = mu dist
    s = sigma dist

-- 多項式基底関数
polynomialBasis :: Floating a => a -> a -> a
polynomialBasis n x = x ** n

-- 確率分布を更新する。演習3.8により共分散と期待値を更新する
-- 第一引数 係数α
-- 第二引数 係数β
-- 第三引数 基底関数にx_n+1を適用した?_n+1縦ベクトル(実装上はM行1列の行列)
-- 第四引数 新しい目標値t_n+1
-- 第五引数 古いガウス分布
update :: Double -> Double -> Matrix Double ->
          Double -> GaussDist -> GaussDist
update alpha beta phi t dist = dist'
  where
    m = mu    dist
    s = sigma dist
    invS' = (inv s) + ((1><1) [beta])*(phi `multiply` (trans phi))
    lhs   = (inv s) `multiply` m
    rhs   = ((1><1) [beta * t]) * phi
    m'    = (inv invS') `multiply` (lhs + rhs)
    dist' = GaussDist m' (inv invS')

-- ノイズを付与する
addNoise :: (RandomGen g, Random a, Num a) =>
            (a, a) -> [a] -> g -> [a]
addNoise range xs g = zipWith (+) xs (randomRs range g)

-- 訓練データを作成する
-- 関数t(x) = 3.14 x + 1.23から 0 <= x < 5の範囲で生成した値に
-- -1.0 から 1.0の範囲の値をランダムに加えた(x,t)のリストを返す。
makeTrainingData :: IO [(Double, Double)]
makeTrainingData = do
  gen <- getStdGen
  let xs = linearScale 1000 (0, 5)
      f  = \x -> 3.14 * x + 1.23
      ts = addNoise (-1.0 ::Double, 1.0::Double) (map f xs) gen
      xts =  zip xs ts
  return xts

-- 事後分布の逐次更新計算を行う
-- 入力された目標値のリストから目標値を一つずつ取り出して事後分布を更新していく。
sequentialUpdateCalc :: Double -> Double ->
                        [(Double -> Double)] ->
                        GaussDist -> [(Double, Double)] -> [[Pos]]
sequentialUpdateCalc alpha beta basiss priorDist xts =
  seqCalc priorDist xts
  where
    update' = update alpha beta
    -- -5 < w0 < 5 , -5 < w1 <5の範囲にある確率P(w1,w2)を全て計算する
    w = [(2 >< 1) [w0, w1] :: Matrix Double
        | w0 <- linearScale 100 (-5, 5)
        , w1 <- linearScale 100 (-5, 5)]
    -- 戻り値で座標wを返す際に行列のままだと表示しづらいので、リストに変換する
    wl = map (concat . toLists) w
    -- 実際の逐次計算は以下の関数で行う。
    -- (x,t)を一つずつ
    seqCalc _ [] = []
    seqCalc priorDist ((x,t):xts) = poss:(seqCalc postDist xts)
      where
        -- 基底関数にxを適用したベクトル?nを計算する
        phi  = (2 >< 1) [f x | f <- basiss]
        -- 全ての(w0,w1)に対して、その座標での確率を計算する
        ps   = map (probability priorDist) w
        -- w0,w1,P(w0,w1)からなるリストを作成する
        poss = zipWith (\ (w0:w1:[]) p -> Pos w0 w1 p) wl ps
        -- 事前分布から事後分布を計算する(演習3.8)
        postDist = update' phi t priorDist

-- 文字列のリストの各要素をファイルに出力する。        
writeSequential name strs = writeSeq 0 strs
  where
    writeSeq _ [] = return ()
    writeSeq n (str:strs) = do
      let fileName = name ++ show n ++ ".dat"
      writeFile fileName str
      writeSeq (n + 1) strs

-- 逐次更新処理を実行する      
runSequentialUpdate = do
  -- 訓練データを作成する
  xts <- makeTrainingData
  let basiss  = map polynomialBasis ([0,1] :: [Double])
      -- 事前分布 P(w) = N(w|m0,s0)
      m0 = (2 >< 1) [0,0]
      s0 = diag . fromList $ [1,1]
      priorDist = GaussDist m0 s0
      -- 更新されていく確率分布を計算する
      posss = sequentialUpdateCalc 1 1 basiss priorDist xts
      -- 計算した確率分布を文字列にする
      strs = map (concat . (map show)) posss
  writeSequential "update" strs
  return ()

2012-11-07

PRML３章の線形回帰を実装した

PRML Haskell 機械学習

PRML３章の線形回帰を実装した。ベイズ線形回帰は未実装。

中谷さんのgihyo.jp連載記事「機械学習　はじめよう」
第11回 "線形回帰を実装してみよう"が非常に参考になりました。

http://gihyo.jp/dev/serial/01/machine-learning/0011?page=1

実装は以下の通り。

import System.Random
import Numeric.LinearAlgebra
import Graphics.Gnuplot.Simple
import Control.Applicative

-- ノイズを付与する
addNoise :: (RandomGen g, Random a, Num a) =>
            (a, a) -> [a] -> g -> [a]
addNoise range xs g = zipWith (+) xs (randomRs range g)

-- 座標を文字列にする
pointsToString :: (Show a) => [a] -> [a] -> String
pointsToString xs ys =
  concat $
  zipWith (\x y -> show x ++ " " ++ show y ++ "\n") xs ys

-- 座標を文字列にする
to3dPosString xs ts = concat $
                      zipWith (\(x,y) z ->    (show x) ++ " "
                                           ++ (show y) ++ " "
                                           ++ (show z) ++ "\n")
                      xs ts

-- ガウス基底関数
gaussianBasis s mu x = exp $ (- (x - mu) ** 2.0) / s

-- 多項式基底関数
polynomialBasis n x = x ** n

-- 二次元用多項式基底関数
-- 多次元の場合にこのような基底関数でよいのかどうかは確認していない。
polynomialBasis' n (x,y) = x ** n + y ** n

-- 基底関数の配列と入力変数xから計画行列Φを計算する
designMatrix basiss xs = (n >< m) [f x | x <- xs, f <- basiss]
  where
    n = length xs
    m = length basiss

-- 正則化係数、計画行列と目標変数の配列から
-- 基底関数の線形結合で用いるパラメータ w を計算する
parameters lambda phi ts = w
  where
    -- 配列を縦ベクトルに変換する
    t       = (length ts >< 1) ts 
    m       = cols phi
    -- (lambda * 単位行列) を計算する
    lambdaI = fromLists $ map (map (*lambda)) (toLists $ ident m)
    -- 計画行列の転置行列を計算する
    phiT = trans phi
    -- w = ((λ I + (Φ^T)Φ)^(-1))(Φ^T)t <= PRML P.143 (3.28)式より
    w = (inv (lambdaI + phiT `multiply` phi))
        `multiply` phiT `multiply` t

-- 線形回帰
-- 正則化係数, パラメータの要素数, 基底関数, 入力値xs,目標値tsから、
-- 線形回帰によって推定された関数fを返す。
linearRegression lambda m basis xs ts = f
  where
    -- 基底関数ベクトルを作る
    basiss = (\_ -> 1.0) : (map (basis . fromIntegral) [1..m])
    -- 計画行列を計算する
    dm = designMatrix basiss xs
    -- パラメータベクトルwを計算してパラメータのリストにする
    ws = concat . toLists $ parameters lambda dm ts
    -- 訓練データから推測される関数を返す
    f  = \x -> sum $ zipWith (*) ws $ map (\g -> g x) basiss

-- xが一次元の場合の計算
test2d = do
  gen <- getStdGen
  let xs = linearScale 1000 (0, 2*pi)
      ts = addNoise (-1.0 ::Double, 1.0::Double) (map sin xs) gen
      f  = linearRegression 0.2 100 (gaussianBasis 0.2) xs ts
      f' = linearRegression 0.2 10 polynomialBasis xs ts
      ys  = map f xs
      ys' = map f' xs
      xtstr  = pointsToString xs ts
      xystr  = pointsToString xs ys
      xystr' = pointsToString xs ys'
  writeFile "train.dat"      xtstr
  writeFile "gaussian.dat"   xystr
  writeFile "polynomial.dat" xystr'


-- xが二次元の場合の計算
test3d = do
  gen <- getStdGen
  let xs = linearScale 50 (0, 2*pi)
      ys = linearScale 50 (0, 2*pi)
      xs' = [(x,y) | x <- xs, y <- ys]
      org = \(x, y) -> cos(0.5 * (x - pi)) * cos(0.5 * (y -pi))
      ts = addNoise (-0.2::Double, 0.2::Double) (map org xs') gen
      f = linearRegression 0.2 10 polynomialBasis' xs' ts
      zs  = map f xs'
      xtstr  = to3dPosString xs' ts
      xzstr  = to3dPosString xs' zs
  writeFile "train3d.dat"      xtstr
  writeFile "gaussian3d.dat"   xzstr
  return ()

main = do
  test2d
  test3d

xが1次元の場合の線形回帰の計算結果は以下のとおり。
ガウス基底関数の分散を大きくすると、もっとスムーズな関数を得る事は出来る。

xが2次元の場合の線形回帰の計算結果は以下の通り。
2次元での多項式基底関数のとり方が正しいかどうか要確認。

2012-11-02

形態素解析エンジンMeCabをFedora17でHaskellから使う際のメモ

Haskell 形態素解析 fedora

Fedora17でのMeCabのインストールとUTF-8用の設定方法、
及びHaskellからMeCabを使う方法を紹介する。

以下のページを参考にさせて頂きました。
http://kane.meta-scheme.jp/article/37183101.html

MeCabのインストールと設定

Fedora17にはmecabパッケージがあるので、
インストールは以下のようにyumからインストールするだけでよい。

yum install *mecab*

上記の方法でインストールしたMeCabは、デフォルトではEUC-JPの辞書ファイルを使うように設定されているため不便である。

MeCabの設定ファイル"/etc/mecabdir"のdicdirで始まる行をコメントアウトして、
以下のようにUTF-8の辞書ファイルを参照するように編集する。

;
; Configuration file of MeCab
;
; $Id: mecabrc.in,v 1.3 2006/05/29 15:36:08 taku-ku Exp $;
;
; dicdir = /usr/lib64/mecab/dic/jumandic-EUCJP
dicdir = /usr/lib64/mecab/dic/jumandic

; userdic = /home/foo/bar/user.dic

; output-format-type = wakati
; input-buffer-size = 8192

; node-format = %m\n
; bos-format = %S\n
; eos-format = EOS\n

UTF-8の辞書ファイルは、上記のyumコマンドを実行した際に/usr/lib64/mecab/dic/jumandicにインストールされたものである。

これで以下のようにmecabが使えるようになった。

$ echo "メカブは美味しくて健康にもいいよ。" | mecab 
メカブ	名詞,人名,*,*,*,*,*
は	助詞,副助詞,*,*,は,は,*
美味しくて	形容詞,*,イ形容詞イ段,タ系連用テ形,美味しい,おいしくて,代表表記:美味しい
健康に	形容詞,*,ナ形容詞,ダ列基本連用形,健康だ,けんこうに,代表表記:健康だ
も	助詞,副助詞,*,*,も,も,*
いい	形容詞,*,イ形容詞アウオ段,基本形,いい,いい,代表表記:いい
よ	助詞,終助詞,*,*,よ,よ,*
。	特殊,句点,*,*,。,。,*
EOS

HaskellのMeCab バインディング Text.MeCab をインストールする

cabalからインストールする。

cabal update
cabal install mecab

HaskellからMeCabを使う。

使い方の概要は以下の通り。

newもしくはnew2を使ってTagger インスタンスを作成する。
Tagger インスタンスと解析対象の文字列を引数に指定してパースする。

以下のMeCab C++サンプルとほぼ同じ動作をするHaskellサンプルを示す。
http://mecab.googlecode.com/svn/trunk/mecab/doc/libmecab.html

サンプルで使用した関数の概要は以下の通り。

new2 :: String -> IO Text.MeCab.MeCab

Tagger インスタンスを作成する関数。
引数には""を指定する。
MeCabの動作を変更するオプションを指定すると思われるが、確認していない。

parse :: MeCabString str => Text.MeCab.MeCab -> str -> IO str

Tagger インスタンスを第一引数に、解析対象の文字列を第二引数に指定すると、
最も精度がよい解析結果の文字列を取得するIOモナドを作る。

parseNBest :: MeCabString str => Text.MeCab.MeCab -> Int -> str -> IO str

精度のよい解析結果を、上位から第二引数で指定した分だけ返す文字列にして返すIOモナドを作る。

parseNBestInit :: MeCabString str => Text.MeCab.MeCab -> str -> IO ()

後述するnext関数で、精度のよい解析結果を順番に取得するための初期化を行う。

next :: MeCabString str => Text.MeCab.MeCab -> IO (Maybe str)

精度のよい解析結果を一つ取り出す。

parseToNode :: MeCabString str => Text.MeCab.MeCab -> str -> IO [Node str]

最も精度のよい解析結果を、文字列ではなく各ノードのリストとして受け取る。

import Text.MeCab
import Control.Monad

toString :: Text.MeCab.Node String -> String
toString node = foldr (\x acc-> x ++ " " ++ acc) [] attrs
  where
    id   = (show $ nodeId node)
    stat = case nodeStat node of
      BOS -> "BOS"
      EOS -> "EOS"
      _   -> (nodeSurface node) ++ " " ++ (show $ nodeRlength node)
    feature = nodeFeature node
    rcAttr  = show $ nodeRcAttr node
    lcAttr  = show $ nodeLcAttr node
    posid   = show $ nodePosid  node
    ctype   = show $ nodeCharType node
    stat'   = show $ nodeStat   node
    isbest  = show $ nodeIsBest node
    alpha   = show $ nodeAlpha  node
    beta    = show $ nodeBeta   node
    prob    = show $ nodeProb   node
    cost    = show $ nodeCost   node
    attrs = [id, stat, feature]
    --attrs = [id, stat, feature, rcAttr, lcAttr, posid, ctype,
    --         stat', isbest, alpha, beta, prob, cost]
    
test = do
  let input = "太郎は次郎が持っている本を花子に渡した。"
  mecab  <- new2 ""
  -- Gets tagged result in string format.
  result <- parse mecab input
  putStrLn $ "INPUT : " ++ input
  putStrLn $ "RESULT: " ++ result

  -- Gets N best results in string format.
  nbresult <- parseNBest mecab 3 input
  putStrLn $ "NBEST: " ++ nbresult

  -- Gets N best results in sequence.
  parseNBestInit mecab input
  forM_ [0..1] $ \n -> do
    nextResult <- next mecab
    case nextResult of
      Just nr -> putStrLn $ (show n) ++ ":" ++ nr
      _       -> return ()

  -- Gets Node object.
  nodes <- parseToNode mecab input
  mapM_ (putStrLn . toString) nodes

  return ()

main = test

2012-10-21

Haskellの線形代数用ライブラリで多変量ガウス分布を計算した際のメモ

Haskell 機械学習 PRML

Haskellで行列計算したい場合にはhmatrixという線形代数用ライブラリが使える。
この記事では、行列計算でよく使いそうな手続きを紹介する。
また、それらの例として多変量ガウス分布を実装する。

hmarixのインストール

例によってcabalでインストールする。

cabal update
cabal install hmatrix

hmatrixの簡単な使い方

例えば3×3の行列を作る場合には、以下のように実行する。

*Main> :m Numeric.LinearAlgebra
Prelude Numeric.LinearAlgebra> (3 >< 3) ([1..9] :: [Double])
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]

2×4の行列を作る場合には、以下のようにする。

Prelude Numeric.LinearAlgebra> (2 >< 4) ([1..8] :: [Double])
(2><4)
 [ 1.0, 2.0, 3.0, 4.0
 , 5.0, 6.0, 7.0, 8.0 ]

fromListsを使えば、aから行列を作ることが出来る。

Prelude Numeric.LinearAlgebra> fromLists ([[1,0], [0,1]] :: [[Double]])
(2><2)
 [ 1.0, 0.0
 , 0.0, 1.0 ]

toListsで行列をaに戻す事も出来る。

Prelude Numeric.LinearAlgebra> let m3x3 = (3 >< 3) ([1..9] :: [Double])
Prelude Numeric.LinearAlgebra> toLists m3x3
[[1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,9.0]]

行列の行数はrowsで求める事が出来る。
行列の列数はcolsで求める事が出来る。

Prelude Numeric.LinearAlgebra> let m2x4 = (2 >< 4) ([1..8] :: [Double])
Prelude Numeric.LinearAlgebra> rows m2x4
2
Prelude Numeric.LinearAlgebra> cols m2x4
4

転置行列はtransで求める。
行列式はdetで求める。
逆行列はinvで求める。

Prelude Numeric.LinearAlgebra> trans m2x4
(4><2)
 [ 1.0, 5.0
 , 2.0, 6.0
 , 3.0, 7.0
 , 4.0, 8.0 ]
Prelude Numeric.LinearAlgebra> det m2x4
*** Exception: det of nonsquare (2><4) matrix
Prelude Numeric.LinearAlgebra> det m3x3
6.661338147750939e-16
Prelude Numeric.LinearAlgebra> inv m3x3
(3><3)
 [ -4.503599627370498e15,   9.007199254740992e15,  -4.503599627370496e15
 ,  9.007199254740996e15, -1.8014398509481984e16,   9.007199254740991e15
 , -4.503599627370498e15,   9.007199254740992e15, -4.5035996273704955e15 ]

行列の積はmultiplyで計算する。

Prelude Numeric.LinearAlgebra> let m1x2 = (1 >< 2) ([1,2]  :: [Double])
Prelude Numeric.LinearAlgebra> let m2x1 = (2 >< 1) ([2,-1] :: [Double])
Prelude Numeric.LinearAlgebra> m1x2 `multiply` m2x1
(1><1)
 [ 0.0 ]

行列から位置を指定して値を取り出す時には@@>演算子を使う。

Prelude Numeric.LinearAlgebra> m3x3 @@> (0,0)
1.0
Prelude Numeric.LinearAlgebra> m3x3 @@> (2,2)
9.0

連立1次方程式も解ける。
以下のような方程式が与えられていた場合には、

$\left\{ 2x_1 + 3x_2 - x_3 = -3\\ -x_1 + 2x_2 + 2x_3 = 1\\ x_1 + x_2 + x_3 = -2\right.$

次のように実行して方程式を解くことが出来る。

Prelude Numeric.LinearAlgebra> let a = (3><3) ([ 2, 3,-1,-1, 2, 2,1, 1,-1] :: [Double])
Prelude Numeric.LinearAlgebra> let b = (3><1) ([-3, 1, -2]:: [Double])
Prelude Numeric.LinearAlgebra> linearSolve a b
(3><1)
 [  1.0
 , -1.0
 ,  2.0 ]

多変量ガウス分布

hmatrixを使って多変量ガウス分布を計算する。
このプログラムをコンパイルして実行すると、"gauss.dat"というファイルに計算結果が出力される。

import Numeric.LinearAlgebra
import Graphics.Gnuplot.Simple

-- mahalanobisDistance :: Matrix a -> Matrix a -> Matrix a -> a
mahalanobisDistance avrg covariance x = ((trans y)
                                         `multiply` (inv covariance)
                                         `multiply` y) @@> (0,0)
  where
    y = x - avrg

-- gaussDist :: Matrix a -> Matrix a -> Matrix a -> a
gaussDist mu sigma x = c * exp a
  where
    d = fromIntegral $ cols x
    c = 1.0 / (((2 * pi) ** (d / 0.5)) * ((det sigma) ** 0.5))
    y = x - mu
    a = -0.5 * mahalanobisDistance mu sigma x

xs =  [(2 >< 1) ([x1, x2] :: [Double]) |
        x1 <- linearScale 200 (-1, 1)
      , x2 <- linearScale 200 (-1, 1)]

averageVector    = (2 >< 1) ([0.0, 0.0]:: [Double])
covarianceMatrix = (2 >< 2) ([0.1, 0.05, 0.05, 0.1] :: [Double])

calcGaussDist avrgVec covMatrix xs = xyz
  where
    values = map (gaussDist averageVector covarianceMatrix) xs
    xl     = map (concat . toLists) xs
    xyz    = zipWith (\(x:y:[]) v -> x:y:v:[]) xl values

toSpaceSepString :: Show a => [a] -> String
toSpaceSepString xs = foldr (\x acc -> (show x) ++ " " ++ acc) "\n" xs

main = do
  let results = calcGaussDist averageVector covarianceMatrix xs
      str     = concat $ map toSpaceSepString results
  writeFile "gauss.dat" str
  return ()

"gauss.dat"が作成されたら、以下のように実行してgnuplotで表示する。

gnuplot> set pm3d
gnuplot> set dgrid3d 51,51
gnuplot> set hidden3d
gnuplot> set terminal png
Terminal type set to 'png'
Options are 'nocrop font /usr/share/fonts/dejavu/DejaVuSans.ttf 12 size 640,480 '
gnuplot>  set output 'gauss1.png'
gnuplot> splot "gauss.dat" u 1:2:3 with lines

画像が作成されたら以下のように画像表示コマンドを実行して画像を表示して確認する。

eog gauss1.png

2012-09-20

多項式曲線フィッティングを実装してみた

Haskell PRML 機械学習

PRMLの"1.1 例:多項式曲線フィッティング"を読んだ。
実際に実装して、科書通りにフィッティング出来ることを確認した。

多項式曲線フィッティングの実装

xを0~1の間で10個選び、それぞれのxに対してsin(2πx)を計算した値にノイズを加えて訓練集合tを作成する。
そして、この訓練集合tを以下の式でフィッティングする。

$y\left( x,w \right) = w_0 + w_1 x + w_2 x^2 + \ldots + w_M x^M = \sum_{j=0}^{M} w_j x^j \hspace{50} (1.1)$

フィッティングは、訓練集合tに対して二乗和誤差が最小になるwを見つける事で行う。
二乗和誤差は以下の式で表される。

$E(w) = \frac{1}{2} \sum_{n=1}^N \{ y(x_n,w)-t_n\}^2 \hspace{50} (1.2)$

(1.2)に(1.1)を代入して、各wで偏微分した結果が0になる時に二乗和誤差が最小である。

$\frac{\partial E(w)}{\partial w_i}=\sum_{n=1}^{N}\{\sum_{j=0}^{M}w_j x_{n}^{j} -t_n\}x_n^i \\ = \sum_{n=1}^{N}\{\sum_{j=0}^{M}w_j x_{n}^{i+j} -t_n x_n^i\} \\ = \sum_{n=1}^{N} \sum_{j=0}^{M} w_j x_n^{i+j} - \sum_{n=1}^{N} t_n x_n^i \\ = \sum_{j=0}^{M}w_j \sum_{n=1}^{N}x_n^{i+j} - \sum_{n=1}^{N}t_n x_n^i \\ = \sum_{j=0}^{M} A_{ij} w_j - T_i = 0$

したがって、
$\sum_{j=0}^{M} A_{ij} w_j = T_i$
$A_{ij} = \sum_{n=1}^{N} x_n^{i+j} \\ T_i = \sum_{n=1}^{N} t_n x_n^i$

のように書く事が出来て(演習1.1)、M+1個の一次方程式から連立一次方程式が得られる。
その連立一次方程式はAw = Tの行列形式にして解く事が出来る。

ただし、このままでは過学習が発生するので、以下のように罰金項を加える。

$\widetilde{E(w)} = \frac{1}{2} \sum_{n=1}^{N} \{y(x_n, w) - t_n\}^2 + \frac{\lambda}{2}\| w\| ^2 \\ \| w \| = w_0^2 + w_1^2 + w_2^2 + \ldots + w_M^2 \hspace{50} (1.4)$

(1.4)も、偏微分して0になる条件から最適なwを求める。
(1.4)に(1.1)を代入して各wで偏微分すると、

$\widetilde{E(w)} = \frac{1}{2} \sum_{n=1}^{N} \{ \sum_{j=0}^{M} w_j x^j -t_n\}^2 + \frac{\lambda}{2} \sum_{m=0}^{M} w_m^2$

$\frac{\partial \widetilde{E(w)}}{\partial w_i} = \sum_{j=0}^M A_{ij} w_j - T_i + \lambda w_i = 0$

となり、
$A_{ij}' = \begin{cases} \sum_{n=1}^{N} x_n^{i+j} & (i \neq j) \\ \lambda + \sum_{n=1}^{N} x_n^(i+i) & (i = j) \end{cases}$
とおくと、
$\sum_{j=0}^{M} A_{ij}' w_j = T_i$
が得られ(演習1.2)、ここから連立一次方程式を解くことでwを得る。

この一連の計算を実装したコードは以下の通りである。

import System.Random
import Control.Monad
import Numeric.LinearAlgebra
import Graphics.Gnuplot.Simple
import Data.List

trainingDataFileName = "training.dat"

-- 多項式の係数を計算 mは次数
polynomialCoefficients :: Field a => Bool -> a -> Int -> [a] -> [a] -> [a]
polynomialCoefficients isRegularized lambda m xn tn = [mtrxW @@> (i, 0) | i <- [0 .. m]]
  where
    m'    = m + 1
    -- Ai,j と Ti は 演習1.1を参照せよ
    a i j = sum $ map (^ (i + j)) xn
    -- 式1.4のように正則化した場合、Ai,jは i == jで係数λが加えられる
    -- 正則化した場合のAi,jの導出は演習1.2を参照せよ
    a' i j | i /= j    =  sum $ map (^ (i + j)) xn
           | otherwise =  lambda + (sum $ map (^ (i + j)) xn)
    t i   = sum $ zipWith (*) tn $ map (^ i) xn
    -- Aw = T の行列として扱い、連立方程式を解く
    fA | isRegularized = a'
       | otherwise     = a
    mtrxA = (m' >< m') [fA i j | i <- [0 .. m], j <- [0 .. m]]
    mtrxT = (m' >< 1 ) [t i| i <- [0 .. m]]
    mtrxW = linearSolve mtrxA mtrxT

-- ノイズを付与する
addNoise :: (RandomGen g, Random a, Num a) => (a, a) -> [a] -> g -> [a]
addNoise range xs g = zipWith (+) xs (randomRs range g)

-- gnuplotでフィッティング曲線, 訓練データ, sin(x)をpng形式で出力する
toFittingCurveScript :: Show a => [[a]] -> String -> String
toFittingCurveScript wss attr = "set xrange [-0.5: 2*pi + 0.5]\n"
                           ++ "set yrange [-1.5: 1.5]\n"
                           ++ "set term png\n"
                           ++ (concat $ map order (zip [0..] wss))
  where
    order (m, ws) = "set output \"order" ++ attr ++ show m ++ ".png\"\n"
                    ++ "set title " ++ "\"order " ++ attr ++ show m ++ "\"\n"
                    ++ "plot " ++ (toPolynomial ws)
                           ++ " w l lt 1 t \"fitting\", " 
                           ++ "sin(x) lt 2, "
                           ++ "\"" ++ trainingDataFileName ++ "\" w p pt 6 lt 3\n"
    toPolynomial xs = concat . map toTerm $ zip [0..] xs
    toTerm (i, coefficient) = "+ (" ++ show coefficient ++ ")" 
                              ++ (concat . take i $ repeat "*x")

erms :: Floating a => [a] -> [a] -> [a] -> a
erms xs ts ws = (2 * e / n) ** 0.5
  where
    f x = (sum . zipWith (*) ws $ scanl (*) 1 $ repeat x)
    e   = sum . map (^2) . zipWith (-) ts $ (map f xs)
    n   = genericLength xs

ermsList :: Floating a => [a] -> [a] -> [[a]] -> [a]
ermsList xs ts wss = map (erms xs ts) wss

valuesToString :: Show a => [a] -> String
valuesToString xs = concat $ map (\x -> show x ++ "\n") xs

main = do
  gen <- getStdGen
  let xs   = linearScale 10 (0,2*pi)
      ts   = addNoise (-0.5 :: Double, 0.5 :: Double) (fmap sin xs) gen
      wss  = map (\m -> polynomialCoefficients False 0   m xs ts) [0 .. 10]
      wss' = map (\m -> polynomialCoefficients True  1.5 m xs ts) [0 .. 10]
      testxs = linearScale 100 (0,2*pi)
      testys = addNoise (-0.5 :: Double, 0.5 :: Double) (fmap sin testxs) gen

  writeFile "plot.plt" $ toFittingCurveScript wss ""
  writeFile "plot_reg.plt" $ toFittingCurveScript wss' "_regularized"
  writeFile "erms_train.dat" . valuesToString $ ermsList xs ts wss
  writeFile "erms_test.dat"  . valuesToString $ ermsList testxs testys wss
  writeFile trainingDataFileName $ concat . map 
                (\(x,y) -> show x ++ " " ++ show y ++ "\n") $ zip xs ts

実行手順は以下の通り。

ghc test.hs
./test
gnuplot -p < plot.plt
gnuplot -p < plot_reg.plt

正則化無しで過学習する(1.1)式は以下のようになる。

正則化ありで、過学習を抑制した結果は以下の通り。

Ermsを計算した結果は、erms_train.datとerms_test.datに書きだすので、
gnuplotで以下のように実行して結果を表示する。

set term png
set output "erms.png"
plot "erms_train.dat" w lp lt 3 pt 6, "erms_test.dat" w lp lt 1 pt 6