2012-12-09

「充足可能性問題(3-SAT)を解く乱択アルゴリズム」 by Haxe

SAT Haxe

初雪です。さっき吹雪いてました。てか3か月ほどブログ更新滞ってました(>_<)
統計とか機械学習の勉強も進んでないのですが、今日はリハビリに、また例の 3-SAT ネタでも投下しときます。

半年以上前に『Rで「充足可能性問題(3-SAT)を解く乱択アルゴリズム」』て記事でこんなこと書いています：

新しい言語を覚えるときに、数学ガール乱択アルゴリズムの「充足可能性問題（3-SAT）を解く乱択アルゴリズム」（p.353）を実装するという癖がついてしましました。

また新しい言語で書いてみました。
今仕事で使っている、Haxe（wikipedia:Haxe）で。

Haxeというのが元々ActionScriptから派生した言語で、基本構文もECMAScriptにほぼ準拠してるので、結局以前に書いた JavaScript 版とかなり近い実装になっちゃいました。今回も焼き直し感が否めないです。

ただ、今までの言語と違う点が1点。
今までの言語（最期の一覧します）は、所謂「スクリプト言語」だったのですが、今回は一応分類上は「コンパイル言語」。今までは言語の動的な性質を利用（eval 関数またはそれに相当するものを利用して文字列を解析・実行等）してきたのですが、今回はそれが使えない*1。
で。
今回は、「クローズ（Clause）を表す文字列を自前で解析」しています。parseClause 関数参照。
てかこの方法は他の言語でも応用できますね。しかも、今までは言語の癖に合わせて入力を少しずつ変えてきたけど、入力形式を統一することも可能ですね。

ちなみに動作は、JavaScript、C++、Java *2 の各プラットフォームで確認しています。

最後に、今までの Rw3sat の実装を一覧列挙。

Ruby版: https://gist.github.com/960956
JavaScript版: https://gist.github.com/966408
CoffeScript版: https://gist.github.com/1000913
R版: https://gist.github.com/2811277
Python版: https://gist.github.com/2818831
Haxe版: https://gist.github.com/4148331

*1:一応、Haxe のコンパイル先を JavaScript にすれば、条件コンパイルや untyped コードを組み合わせて JavaScript のコードも書ける、つまり eval も使えるのですが、クロスプラットフォームで動くものにしたかったので今回は使用していません。

*2:日本語が文字化けしますが、一応動作します。

2012-09-02

RMeCab と RCaBoCha をインストールしてみた

R MeCab CaboCha

えと、今日は取り敢えずインストールの記録だけ。

前準備

R はずいぶん前にインストールしてあったし、MeCab と CaboCha もインストール済。
なので今回は難しいことはほとんどありませんでした。
一応、インストール環境。

OS: Mac OS X Lion 10.7.4
R: ver.2.15.1 (R64)
MeCab: ver.0.994
CaboCha: ver.0.64

RMeCab のインストール

RMeCab のサイトから Mac Lion 用のバイナリパッケージをDLして、同ページ内「Macintosh 版バイナリのインストール方法」に従ってインストール。
ただ一部うまくいかなかったので、うまくいった方法を以下に記述します。

DLしたアーカイブ（RMeCab_0.995.tgz）を適当なフォルダに解凍
R64(GUI)のメニューから [パッケージとデータ]→[パッケージインストーラ]を選び、「CRAN(バイナリパッケージ)」となっているドロップダウンを「このコンピュータ上のパッケージディレクトリ」を選択*1
「インストール…」ボタンをクリックして、 .tgz ファイルを解凍して出来たフォルダを選択。
→インストール完了。特に問題なし。

動かしてみました。

> library(RMeCab)
> res <- RMeCabC("すもももももももものうち")
> res
[[1]]
    名詞 
"すもも" 

[[2]]
助詞 
"も" 

[[3]]
  名詞 
"もも" 

[[4]]
助詞 
"も" 

[[5]]
  名詞 
"もも" 

[[6]]
助詞 
"の" 

[[7]]
  名詞 
"うち" 

> res2 <- unlist(res)
> res2
    名詞     助詞     名詞     助詞     名詞     助詞     名詞 
"すもも"     "も"   "もも"     "も"   "もも"     "の"   "うち" 
> res2[names(res2) == "名詞"]
    名詞     名詞     名詞     名詞 
"すもも"   "もも"   "もも"   "うち" 
>

うん、簡単ですね。

RCaBoCha*2 のインストール

こちらも同じ人が同じサイトで公開しているのですが、実はちょとバージョンが古く、そのままではうまく動きませんでした。例によって成功例を先に載せて後で Point としてどこで引っかかって何をしたか記します。

基本は RCaBoCha のサイトからバイナリパッケージをDLして、同ページ内「Macintosh 版バイナリのインストール方法」に従ってインストール。詳細は以下。

DLしたアーカイブ（RCaBoCha_0.29.tgz）を適当なフォルダに解凍
R64(GUI)のメニューから [パッケージとデータ]→[パッケージインストーラ]を選び、「CRAN(バイナリパッケージ)」となっているドロップダウンを「このコンピュータ上のパッケージディレクトリ」を選択
「インストール…」ボタンをクリックして、 .tgz ファイルを解凍して出来たフォルダを選択。
インストール完了後、ターミナルを開いて以下を実施：
- /Library/Frameworks/R.framework/Versions/ に移動 (cd /Library/Frameworks/R.framework/Versions/)
- ディレクトリ 2.13 が存在しないことを確認し、2.15（現在利用している R のバージョンディレクトリ）のシンボリックリンクとして 2.13 を作成 (sudo ln -s 2.15 2.13)

$ cd /Library/Frameworks/R.framework/Versions/
$ ls -la
total 8
drwxrwxr-x  5 root  admin  170  9  1 17:00 .
drwxrwxr-x  8 root  admin  272  9  1 16:59 ..
drwxrwxr-x  3 root  admin  102  4  5 21:39 2.14
drwxrwxr-x  7 root  admin  238  9  1 22:08 2.15
lrwxr-xr-x  1 root  admin    4  9  1 16:59 Current -> 2.15
$ sudo ln -s 2.15 2.13
$ ls -la
total 16
drwxrwxr-x  6 root  admin  204  9  1 22:08 .
drwxrwxr-x  8 root  admin  272  9  1 16:59 ..
lrwxr-xr-x  1 root  admin    4  9  1 22:08 2.13 -> 2.15
drwxrwxr-x  3 root  admin  102  4  5 21:39 2.14
drwxrwxr-x  7 root  admin  238  9  1 22:08 2.15
lrwxr-xr-x  1 root  admin    4  9  1 16:59 Current -> 2.15

Point:

用意されている RCaBoCha が R ver.2.13 での利用を想定したもの*3で、そのまま実行しようとすると、/Library/Frameworks/R.framework/Versions/2.13/Resources/lib/libR.dylib を読みに行こうとしていた（R2.13は利用しておらずそのバージョンディレクトリがないので、「libR.dylib が読み込めない」旨のエラーメッセージが出た）
→ 同名のファイルが …/2.15/Resources/lib/ 内に出来ていたので、こっちを見に行ってくれるように、バージョンディレクトリを偽装。

※もし R 2.13 も過去に利用していて、同じエラーメッセージが出たら、/Library/Frameworks/R.framework/Versions/2.13/Resources/lib/ ディレクトリ内で ln -s /Library/Frameworks/R.framework/Versions/2.15/Resources/lib/libR.dylib すれば良いんじゃないかと思います。

さてこれでインストールできましたが、RMeCab をロードしたままの状態で RCaBoCha をロードしようとすると、

> library(RCaBoCha)

 次のパッケージを付け加えます: '‘RCaBoCha’' 

The following object(s) are masked from ‘package:RMeCab’:

    entropy, globalEntropy, globalIDF, globalIDF2, globalIDF3,
    localBin, localLogTF, localTF, mynorm, removeInfo

>

こんな感じに警告が出ます。つまり関数名がかぶってるんですね。さすが同じ作者。
ということで、実際に使用するときは、RMeCab と RCaBoCha は同時にロードしないようにするか、両方必要なときは RMeCab::entropy() とか RCaBoCha::mynorm() とかパッケージ名修飾するですね。

さてさて、色々実験。

> RCaBoCha("それは面白い本であった。")
   Term1  Term2    POS D1 D2
1   それ   それ   名詞  0  0
2     は     は   助詞  0  0
3 面白い 面白い 形容詞  0  0
4     本     本   名詞  0  0
5     で     だ 助動詞  0  0
6   あっ   ある 助動詞  0  0
7     た     た 助動詞  0  0
8     。     。   記号  0  0
FROMAT_TREE =
       それは---D
        面白い-D
    本であった。
EOS
> RCaBoChaFreq("それは面白い本であった。しかし、この本に比べると面白くはない。")
string 2 = "ない" setted: length = 0 
          Term           Pos Freq
1      、+ない   記号+形容詞    1
2           。          記号    2
3         ある        助動詞    1
4    この+ない 連体詞+形容詞    1
5  しかし+ない 接続詞+形容詞    1
6         それ          名詞    1
7           た        助動詞    1
8           だ        助動詞    1
9      と+ない   助詞+形容詞    1
10     に+ない   助詞+形容詞    1
11          は          助詞    1
12     は+ない   助詞+形容詞    1
13          本          名詞    1
14     本+ない   名詞+形容詞    1
15 比べる+ない   動詞+形容詞    1
16      面白い        形容詞    1
17 面白い+ない 形容詞+形容詞    1
>

基本的な使い方は問題ないようです。が…。

問題発生

同サイトで用意されているサンプルデータで動作を試してみて、いくつか動作しないものが。例えば。

> dat <- read.csv("H18koe.csv")
> res <- RCaBoChaDF(dat[,"opinion"])
no terms larger than  minFreq = 1
 以下にエラー RCaBoChaDF(dat[, "opinion"]) : give less number to minFreq!
>

つまり読み解くと、「出現頻度が（指定された）最小値（＝1）より大きな言葉がありません」というエラー。おそらく RCaBoCha が出しているエラーだと思われます。
でも元データは作者が用意しているデータで、明らかに2回以上出てくる言葉がたくさんちりばめられています。
これは、ウチの CaboCha に何か問題がある？それとも RCaBoCha の不具合？
後者だとしたら、1年もほっとかないでなんとかしてほしいなー。

*1:RMeCab サイト上の説明では、「このコンピュータ上のバイナリパッケージ」を選択して、DLした .tgz ファイルを指定すればインストールできるとあるのですが、うまくいかなかったので解凍してそのディレクトリを指定する方法にしました。

*2:ところでなんで RCaboCha じゃなくて RCaBoCha なんでしょう。"B" の大文字小文字を何度打ち間違えたことか…。

*3:てか1年ちょっと前で更新が止まってる…。

2012-08-13

ナイーブベイズフィルターの実装と考察

機械学習 Ruby MeCab 勉強

やっぱ実家はいろいろ捗る。洗濯とか料理とか自分では何もやらなくていいから。
というわけで、前回の勉強の続き、予告通りナイーブベイズフィルターを実装してみました。

実装

と言っても、実装に関しては手抜きです。
gihyo.jpの連載機械学習はじめようの、第3回「ベイジアンフィルタを実装してみよう」のサンプルコードをほぼそのまま流用しました。
ただし。以下のような手を加えています。

せっかくなので別言語(=Ruby)に翻訳することで、一応自分でコードを写経して覚えるという行動を取る。
せっかくなので先日インストールした MeCab を利用するようにコードを修正する。

どちらもすんなりうまくいったのですが、後者でちょっと躓きました。それは追々触れていきます。

まずは翻訳（写経）*1。

あ、もちろん morphological.rb の動作には MeCab と MeCab Ruby のインストールが前提です。

さて、これをそのまま動かしてみたのですが…。

$ ruby naivebayes.rb
ヴァンロッサム氏によって開発されました. => 推定カテゴリ: Ruby
豊富なドキュメントや豊富なライブラリがあります. => 推定カテゴリ: Python
純粋なオブジェクト指向言語です. => 推定カテゴリ: Ruby
Rubyはまつもとゆきひろ氏(通称Matz)により開発されました. => 推定カテゴリ: Ruby
「機械学習 はじめよう」が始まりました. => 推定カテゴリ: 機械学習
検索エンジンや画像認識に利用されています. => 推定カテゴリ: 機械学習

最初の文章が、期待通り「Python」に分類されず「Ruby」と分類されてしまっています。
他の「『機械学習はじめよう』のサンプルを試してみた」サイトをいくつか見てみても、きちんと Python と分類されているところしか見付からないのに、なぜだろう。ロジックは間違っていないし。
ということで、色々見てみました。

MeCab の形態素解析の精度、あるいは、クセ？

色々調べた結果、以下のことが分かりました。

まず、Python→Ruby に翻訳したコードには問題なし（ロジック的に）。
主な原因は全て、MeCab 側にある、と。
その問題点は、いくつか出てきたのですが大きなのは以下の2つ。

人名「グイド・ヴァンロッサム」が、この塊で一般名詞として認識されている。
- 特に「ヴァンロッサム」が単語として認識されていないのでカウントに引っかからない。
日本語的には助動詞の「れ(る)」「られ(る)」が、動詞（接尾）として認識されている。
- 動詞は検索対象なので、これが含まれる量がカテゴリ推定に大きく左右されてしまっている。

つまり簡単に言うと、

「ヴァンロッサム」という単語が出てきたがこれは学習データに存在しない → 推定材料にならない。
「開発された」の「れ（動詞）」が、Ruby カテゴリの学習データに一番多く存在している → 推定材料。

ということで結果として Ruby に推定分類されてしまった、というわけ。

MeCab 側の解決方法模索

まず、MeCab の解析精度にも問題があることは確かです。
考えてみれば先日インストールしてから全然触ってないですし、辞書も mecab-ipadic をデフォルトのままいじってないですから。
あと morphological.rb での品詞分類も MeCab デフォルトの品詞分類に基づいたモノになっています。

このあたりの、設定を弄るなり、学習データを用意して学習させるなり、mecab-dic-overdrive等のツールを使って辞書へのパッチ適用・ノーマライズ、ユーザ辞書の拡充、などを行っていくべきでしょう。行く行くは。

ま、それはそれとして、課題として残しておいて、今回は別のアプローチを考えます。

学習

MeCab にも「学習課題」が出てきたわけですが、そっちのせいにしないで、自分自身にもっと何かできることを探す、というのがやっぱ筋でしょうね。
じゃ、今回実装したフィルターに、「学習」させれば良い。
そもそも現在の学習データが各カテゴリ1件ずつというのは少なすぎますからね。
連載記事の最後の方にもきちんと書いてあります。

訓練データが増えることによって，より正確な分類ができるようになるので興味のある方はご自身で試してみてください

ということで、訓練データを追加してみましょう。
というか、naivebayes.rb には追加訓練データがコメントアウトしてあります。

  #nb.train("ヴァンロッサム氏によって開発されました。", "Python")

そう、正しく分類されなかったら、その正しい分類先を教えてやる。これが機械学習の「学習」の基本。

ということでこのコメントアウトを外して実行してみました。

$ ruby naivebayes.rb
ヴァンロッサム氏によって開発されました. => 推定カテゴリ: Python
豊富なドキュメントや豊富なライブラリがあります. => 推定カテゴリ: Python
純粋なオブジェクト指向言語です. => 推定カテゴリ: Ruby
Rubyはまつもとゆきひろ氏(通称Matz)により開発されました. => 推定カテゴリ: Ruby
「機械学習 はじめよう」が始まりました. => 推定カテゴリ: 機械学習
検索エンジンや画像認識に利用されています. => 推定カテゴリ: 機械学習

無事、同文章が Python に分類されました(^-^)

でも、これで「めでたしめでたし」じゃないんですけどね。「れ(る)」「られ(る)」問題はまだ残ってます。
例えば「必要に迫られて開発された発明品。」という文章*2を分類してみると、この時点での結果は「Python」となります。
同じ文章を、先ほどコメントアウトを外した追加訓練データをもう一度コメントアウトして実行すると、今度は「Ruby」になります。

ここまで来ると、やっぱり MeCab 側の調整か、プログラム側の調整が必要になってくるでしょうね。
元々連載記事の方にも、推定フェーズの出現頻度のところで「ゼロ頻度問題」の解決策とその精度問題のことが取り上げられています。今回は同連載のコードをそのまま翻訳したので、決して精度の高くない「加算スムージング」という補正方法を用いています。これをもっと精度の良いモノに置き換えることで解決する問題とかもきっと出てくるでしょう。

まとめ

勉強と言っても今回はほぼ写経で終わりましたが、これで大体「どういうことをすれば単純分類器が作れるか」ということは掴めたと思います。
もっと実用的な or 個人的ニーズのある文例の分類を試しながら、色々調整していこうと思います。

次回

次…どうしよう。
連載機械学習はじめようの流れに沿って、次の記事、理論編の分布の話を掘り下げてみるか。
それとも、ここでもうちょっと統計よりに方向転換を図るか…。

*1:gist に上げてありますので、ご自由に Fork してご使用ください。もっとこうしたら良いよとかあればそれも大歓迎。

*2:この例文には意味はありません。主に「れ(る)」「られ(る)」によって分類先の推定が左右されることを見るためだけの例文です。

2012-08-12

Four Fours by @miruka_bot

整数ゲーム数学ガール

8月限定で、Twitter上でちょっと面白いゲームが開催されています。


      【夏休み企画】4つの4 (Four Fours)で表そう。ハッシュタグ #FourFours を付けて、私 @miruka_bot 宛に。重複しない正しい式はリツイートします。詳しくは URL へ。今日(8月12日)は 12 です。
2012-08-12 07:12:21 via ミルカさん bot
      使えるのはすべて半角で『4つの4』と『括弧』と『四則演算(+-*/)』、並べる(44)、小数点(.4)、余剰(%)、冪乗(**)、階乗(4!,(式)!)。平方根、循環小数などは不可です。締め切りは23時まで。#NotRT がある時はリツイートしません。 #FourFours
2012-08-12 07:12:54 via ミルカさん bot

2年前にも同じ企画を行っていて、その時は全日参加した（全日リツイートしてもらえたかは記憶が定かでは無い）のですが、今年は出遅れました(^-^; というか2年前に出た答えを知っているのでそれが出てなかったらチャチャ入れる程度の参加にしようかな、と。

取り敢えず今日（答えが12になる式）2つだけ投げておきました。
今日は答えがたくさんある日だし日曜日だしお盆だし、ご興味と暇があればみなさんも腕試しに参加してみてください(^-^)

少しだけ考察

この問題、まとめると、以下のようになります。

数字の「4」を4つ使う（10進数）。他の数字は一切使用しない。
以下の演算及び記号を用いて、与えられた数（今日なら 12）を算出する式を作成する。
使える演算記号その他は以下の通り：
- 括弧（「(」、「)」）
- 四則演算（「+」、「-」、「*」、「/」）
- 剰余（「%」）
- 小数点（「.」、ただし「.4」と書いてもOK（0.4という数とみなされる））
- 並べて2桁以上にする（「44」、「444」、「4444」のほか、「4.4」「.44」などもOK）
- 冪乗（「**」、4**4 なら 256）
- 階乗（「!」、4! なら 1*2*3*4=24）

これ、もっと制限を厳しくしたタイプも考えられます。つまり使える演算を少なくする、ということ。例えば、

- 剰余（「%」）

これ意外と自由度が高く、凄くバラエティに富んだ計算式が生まれます。
逆にこれを条件から外しても、問題としては成立します。少なくとも1〜31の整数は剰余無しですべて計算式を作ることが可能です（ @antimon2 調べ）。

さらに四則演算だけに絞ったり、並べるのを不可（「4」または「.4」のみ許可）にするなど、の制限もアリだと思います。

逆に制限をもっと緩くすることも可能例えば、

- 平方根（「√」、√(4) なら 2）
- log（自然対数、log(4) は 1.386294…）

を使えるようにするとか。
「これ何の役に立つの？」て思われるかも知れませんが、実はこの2つが使えちゃうと、どんな整数も作り出すことが出来ちゃいます(^-^;

log(log(4)/log(4))/log(log(√(4))/log(4)) = 0
log(log(√(4))/log(4))/log(log(√(4))/log(4)) = 1
log(log(√(√(4)))/log(4))/log(log(√(4))/log(4)) = 2
 :
log(log(√(√(…√(4)…)))/log(4))/log(log(√(4))/log(4)) = n（√のネストがn回）

参考：数学パズル「４つの４」の進展：の下の方

まーこれでは問題としては面白くないので、

- logは使用不可
- 平方根（「√」）はOK。
- その代わり『並べる』か『階乗（「!」）』のどちらかを使用不可

あたりにすれば、ゲームバランスが大きくは崩れず面白くなるんじゃないかなー、と個人的には思います。

2012-08-11

フィボナッチ数(2) ver.末尾再帰

フィボナッチ数 Ruby 整数

暑中残暑お見舞い申し上げます。
この記事書き始めた頃は暑中だったのに色々あって立秋過ぎちゃった。

さておき、前回予告した、「フィボナッチ数を求めるための末尾再帰」について。

おさらい …フィボナッチ数

詳細は、前回の記事とか、適当にググってみたりとかしてください。

簡単に言うと、

$F_0 = 0, F_1 = 1, F_n = F_{n-1} + F_{n-2}(n\geq2)$

で定義される数列で、例えば

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0  # for negative n
  a, b = 0, 1
  n.times { a = b + b = a }
  a
end

こんな関数（メソッド）を書くと一般の n 番目のフィボナッチ数を求められます。
ちなみに最初の方を列挙すると以下の通り。

(0, )1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, …

末尾再帰について少しだけ

一般に関数の「再帰呼び出し」というのは、関数の実装の中で自分自身を呼び出す手法のこと。
「末尾再帰」というのは、その再帰呼び出しを関数の最後のステップで（のみ）行う手法のこと。

例えば、再帰呼び出しを利用した階乗の定義を考えると分かりやすいですね。
↓

def factorial n
  return 1 if n < 2
  n * factorial(n-1)
end

最後に自分自身を呼び出しています。逆にそれ以外の部分では余分な再帰呼び出しは行っていません。
これも末尾再帰と考えることも出来ますが、さらに推し進めて「最後の行は本当に再帰呼び出ししか行わない」ようにしてみましょう。この場合↓こうなります。

def factorial n, r=1
  return r if n < 2
  r *= n
  n -= 1
  factorial n, r
end

本当に最後の行が、元の関数の定義式とほぼ同じになるように、少し冗長に書きました。とにかく最後の行は、パラメータを少し変えて自分自身を再帰呼び出ししている、ただそれだけのことをしているのが分かると思います。
これが厳密な「末尾再帰」のカタチ。

この形になっていると、再帰呼び出しをジャンプに変換することで実行効率を上げる、所謂「末尾最適化」が可能になる*1、という特徴もあるのですが、今回はその話はパス。本題はそっちじゃないので。

フィボナッチ数を末尾再帰で(1) O(n) バージョン

まず前回考えた「関数指向の例」を考えてみましょう。
コードを再掲します。

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0  # for negative n
  return n if n < 2
  fib(n - 2) + fib(n - 1)
end

この例も、再帰呼び出しになっています。しかし最後のステップは、自分自身を2回呼び出しています。
これは末尾再帰のカタチではありません。

これをなんとか末尾再帰のカタチに持って行きたいと思います。どうすれば良いでしょう？

ここで、前節の階乗の例で、2つの実装がどう違っていたのかを思い出してみましょう。

最初の例の最後のステップは、詳しく見るとこうなっています。

再帰呼出しで、nが一つ少ない場合の階乗を実際に計算している。
それを計算した上で、階乗の定義に合わせて今回のnの場合の答えを算出している。

この仕組み、先ほどのフィボナッチ数の関数と同じ構造になっていますよね。

再帰呼出しで、nが1少ない場合と2少ない場合のフィボナッチ数を実際に計算している。
それを計算した上で、フィボナッチ数の定義に合わせて今回のnの場合の答えを算出している。

それに対して、階乗の2つめの例は、最後のステップがこんな構造になっています。

最後のステップに行く前に、カウンタとなる変数と、実際の値を計算するパラメータを変更している。
新しいパラメータで再帰呼出している。

この「実際の値を計算するパラメータ」が重要で、階乗の例の場合、実際に階乗をステップごとに計算しています。
ただし、先の例は「1に、1〜nをこの順に掛けている」形になっているのに対して、後の例は「1に、n〜1をこの順に掛けている」形になっています。
でもこれは大きな問題ではありません。この場合、結果はどちらも同じになりますから。

それを踏まえて。
フィボナッチ数の場合は、どうすれば良いのかというと。
まず、カウンタ以外のパラメータは少なくとも2つ必要です。なぜなら、フィボナッチ数の定義に則った計算を行う以上、少なくとも2つのフィボナッチ数からしか算出できない構造になっているのですから。
そして、その2つの値を使って計算して、その結果をまたパラメータとして次の再帰呼出に与える、ということを繰り返していくことを考えればOK、というわけです。

その2つをこのようにすれば、実はうまくいきます。

def fib n, a=0, b=1
  return n.even? ? -fib(-n) : fib(-n) if n < 0  # for negative n
  return a if n < 1
  #return b if n < 2
  fib n-1, b, a+b
end

つまり。
最初はa=0(=F₀),b=1(=F₁)とし、再帰呼出のたびに、aにb、bにa+bを代入することでF_n-2とF_n-1をF_n-1とF_nに変換している、というわけ。
これを所定の回数繰り返せば、aにF_nの計算結果が入っているはずだから、それを返して終わり。

もちろんすぐに分かるように、この関数の場合の計算量は O(n) です。単純にn回繰り返しているだけですからね。

改良(?)バージョン

少し考えると、この方がもっと効率よいんじゃないか？という以下の実装も思いつきます。

def fib n, a=0, b=1
  return n.even? ? -fib(-n) : fib(-n) if n < 0  # for negative n
  return a if n < 1
  return b if n < 2
  fib n-2, a+b, a+b+b
end

つまり1回の再帰呼出で、カウンタを1ずつ減らす代わりに2ずつ減らして、与えるパラメータも「F_n-2とF_n-1をF_nとF_n+1に変換」する、というわけ。
計算量はやはり O(n) ですが、実際は、というと。

実は実行速度の面ではむしろ遅くなります。1回繰り返すときに、先の例では足し算1回だけだったのに対して、今回は1回繰り返すとき（先の2回分に相当）に足し算3回。つまり1.5倍になってしまうからです。
ただし再帰呼出の回数は減るので、関数呼び出しのオーバーヘッドとか、スタックオーバーフローの問題は軽減できます。ただこの話も今回の本題ではないので言葉だけ。

フィボナッチ数を末尾再帰で(2) O(logn) バージョン

次は、倍数公式をうまく活用して計算量が O(logn) になるようにできないか考えてみましょう。
「倍数公式」とは、以下のようなモノでした。

$F_{2n}=(2F_{n-1}+F_n)F_n\\F_{2n+1}=F_n^2+F_{n+1}^2$

実はこれをもっと一般化した、以下のような「和の公式」が存在します。

$F_{n+m}=F_nF_m+F_{n-1}F_m+F_nF_{m-1}\\F_{n+m-1}=F_nF_m+F_{n-1}F_{m-1}$

これを利用するとさらにこんな風に書けます：

$F_{2^k+m}=F_{2^k}F_m+F_{2^k-1}F_m+F_{2^k}F_{m-1}$

ここで0≦m＜2^kと考えれば、nを2進数表記して、それを利用してまずF_nをF_2^kとF_m（およびその周辺のフィボナッチ数）で表して、F_mをまた同様に分解して…という繰り返しで計算していくことが出来そうです。

実際にやってみたのが、こんな例になります：

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0
  fibFastIter = lambda do |a, b, p, q, c|
    return a if c < 1
    if c.even?
      a1, b1, p1, q1, c1 = a, b, p**2 + q**2, (2*p+q)*q, c.div(2)
    else
      a1, b1, p1, q1, c1 = a*p + b*q, a*q + b*q + b*p, p, q, c-1
    end
    fibFastIter[a1, b1, p1, q1, c1]
  end
  fibFastIter[0, 1, 0, 1, n]
end

ちょっと先ほどまでと書き方が変わってますが、自分自身を再帰関数としていたか、別の再帰関数を（ラムダ式で）用意したかの違いだけなのであまり気にしないでください。
p, q の方が、倍数公式を利用した、F_2^kを表すフィボナッチ数となっており、
a, b の方が、和の公式を利用した、F_2^k+mを表すフィボナッチ数となっています。
で、繰り返していって最後に、aに求めるフィボナッチ数が入っている、という仕組み。

これ、結構速いです。前回の記事の『メモ化＋倍数公式』に肉薄するくらい。どちらも O(logn) ですしね。
ベンチマークを取ってみました。

                      user     system      total        real
fib(50)x1000      0.010000   0.000000   0.010000 (  0.008081)
fib(100)x1000     0.010000   0.000000   0.010000 (  0.010353)
fib(500)x1000     0.020000   0.000000   0.020000 (  0.025588)
fib(1000)x1000    0.040000   0.000000   0.040000 (  0.033170)
fib(5000)x1000    0.080000   0.000000   0.080000 (  0.082407)
fib(10000)x1000   0.190000   0.000000   0.190000 (  0.192046)

与える n が十分に小さいときは、メモ化＋倍数公式の時よりも速くなっているのが分かります。
でも n=5000 の時に同じ数値となり、n=10000 のときは、逆に少し時間かかってますね。
それは、この方式にまだ少し無駄があるから。

改良(!)バージョン

実は先のコードは、ググって見付けたHaskelか何かの関数型言語の実装例をRubyに翻訳したモノです。
でも、明らかな無駄があるの、分かりますよね？
カウンタ変数 c を、どんどん2で割れば良いのに、奇数の時は1引いて再帰呼び出しして、偶数の時に2で割っています。これをまとめちゃえば、効率よくなるはず。
ということでまとめちゃいました。

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0
  return 0 if n < 1
  fibFastIter = lambda do |a, b, p, q, c|
    return a*p + b*q if c < 2
    a, b = a*p + b*q, a*q + b*q + b*p if c.odd?
    fibFastIter[a, b, p**2 + q**2, (2*p+q)*q, c.div(2)]
  end
  fibFastIter[0, 1, 0, 1, n]
end

ベンチマーク。

                      user     system      total        real
fib(50)x1000      0.000000   0.000000   0.000000 (  0.005049)
fib(100)x1000     0.010000   0.000000   0.010000 (  0.006286)
fib(500)x1000     0.020000   0.000000   0.020000 (  0.016950)
fib(1000)x1000    0.030000   0.000000   0.030000 (  0.029278)
fib(5000)x1000    0.060000   0.000000   0.060000 (  0.062313)
fib(10000)x1000   0.130000   0.000000   0.130000 (  0.132980)

うん。やっぱり先のモノより速くなりました(^_^) n=10000 の時も『メモ化＋倍数公式』の時と同じ数値を出しています。
ていうか、これ以上大きな n のことを考えると、やっぱりメモ化を利用した場合の方が速くなる傾向にあるってことですよね、これ。
メモ化利用の方が再帰呼出の回数は多いけれど、演算（足し算や掛け算）の総回数は下回っていくから、なんでしょうね。

*1:実際には厳密バージョンではなくても末尾最適化は適用可能です。

2012-07-25

最適化の手前の数学

整数フィボナッチ数

というタイトルで、発表してきました。
一昨日のイベント数学ナイト2 in 名古屋GeekBar にて。

その時の発表資料を↓に埋め込んでおきます。

最適化の手前の数学 from Shunsuke GOTOH

内容のメインは、「フィボナッチ数を高速に求める方法」。
このブログの過去記事を再構成したモノです。→ フィボナッチ数 - 名古屋で数学するプログラマ(仮)

ただ。
この2か月前の記事似ない内容も盛り込まれています。
それは、「末尾再帰」について。

近いうちに、このスライドで紹介しているフィボナッチ数を求めるための末尾再帰について、それと「末尾再帰」そのものについて深く掘り下げた内容について、記事にしたいと思っています。

2012-07-15

ベイズの展開公式

確率勉強

前回から2か月以上間が空いてしまいました。
勉強の続きです。というかまともな連休がないと勉強が進まない自分をなんとかしたい…。

事前確率・事後確率

事前確率・事後確率というのは、ベイズの定理やその応用を勉強していて初めて触れた言葉です。

事前確率: 先に知られている（あるいは"事前に与える"）確率。
事後確率: 条件を加えて算出された条件付き確率。

ベイズの公式：

$\displaystyle P(X|Y) = \frac{P(Y|X)P(X)}{P(Y)}$

で言うと、

事前確率…P(X)
事後関数…P(X|Y)

となります（確率変数 X に着目した場合）。

確率変数の独立性

確率変数の独立性: 確率変数 X と Y について、P(X|Y) = P(X) が成立するとき、X と Y は「独立」である、と言う。

※P(X|Y) = P(X) なら、乗法定理より P(Y|X) = P(Y) も言えます。

この定義から、以下が成り立ちます。

確率変数 X と Y が独立ならば、 P(X, Y) = P(X)P(Y) 。

証明は乗法定理からあっという間なので省略。
この式を利用して、独立ではない（もしくは独立性が確認されていない）2つの確率変数の同時確率を P(X)P(Y) で計算することもあるそうです。多少の誤差が生じるけれど「パラメータが多すぎて計算できない」よりはマシ、ということで。

原因の確率、結果の確率

ベイズの公式に意味付けします。
今まで確率変数として、無機質な X, Y を使ってきましたが、ここで一時的に H, D に置き換えます。その意味は以下の通り：

H: 原因や仮定（Hypothesis）
D: H のもとで得られる結果やデータ（Data）

ベイズの公式は以下のようになります。

$\displaystyle P(H|D) = \frac{P(D|H)P(H)}{P(D)}$

この式に現れる確率のウチ、

P(D|H)…結果の確率（原因Hが与えられたときに得られたデータDの確率なので）
P(H|D)…原因の確率（データDが得られたときの原因がHである確率なので）

と呼びます。
つまりベイズの公式は、結果の確率から原因の確率を導出する公式、になっているというわけですね。

ベイズの展開公式

今、原因が複数考えられて、それらを H₁, H₂, … , H_n とします。
結果Dは、原因 H₁〜H_n のいずれかから起きるものであることが分かっていて、漏れはない（つまり全射）とします。

さらに。
H₁〜H_n が互いに独立と仮定します。

すると以下が成り立ちます。

$P(D) = P(D, H_1) + P(D, H_2) + \cdots + P(D, H_n)$

さらに乗法定理より

$P(D) = P(D|H_1)P(H_1) + P(D|H_2)P(H_2) + \cdots + P(D|H_n)P(H_n)$

これを使ってベイズの公式を変形します。
確率変数 H_i(i=1, 2, … , n) について、

$\displaystyle P(H_i|D) = \frac{P(D|H_i)P(H_i)}{P(D|H_1)P(H_1) + P(D|H_2)P(H_2) + \dots + P(D|H_n)P(H_n)}$

ここでもう一つの言葉を定義。

尤度: 事前確率から事後確率を算出する際の「もっともらしさ」を表す確率。

ベイズの展開公式（H_iに着目）で言うと、

事前確率…P(H_i)
尤度…P(D|H_i)
事後確率…P(H_i|D)

となります。
先ほど定義した別の言葉がありましたが、それらとの関係は

結果の確率＝尤度
原因の確率＝事後確率

となります…なんか字面だけ見るとちぐはぐしているように見えますね。
文脈で呼ばれ方が変わってくるワケなのですが、詳細はまた機会があれば。

応用例

こんな問題を考えます*1。

[問題]
赤玉と白玉合わせて3個入った壺が3つある。壺1には赤玉が1つ、壺2には赤玉が2つ、壺3には赤玉が3つ入っている。今どれか1つの壺から1個取り出してみたところ、赤玉であった。取り出した壺が壺3である確率を求めよ。

これをベイズの展開公式にあてはめると、以下のようになります。

P(壺3|赤) = P(赤|壺3)P(壺3) / (P(赤|壺1)P(壺1)+P(赤|壺2)P(壺2)+P(赤|壺3)P(壺3))

つまり

事前確率…P(壺3)
尤度…P(赤|壺3)
事後確率…P(壺3|赤)

であり、この問題は「事前確率や尤度から、事後確率を求める問題」になっていることはよく分かります。

まず、尤度（＝結果の確率）から見てみましょう。
これはすぐに分かります。

P(赤|壺1)…壺1から取り出した玉が赤玉である確率＝1/3。
P(赤|壺2)…壺2から取り出した玉が赤玉である確率＝2/3。
P(赤|壺3)…壺3から取り出した玉が赤玉である確率＝3/3（=1）。

次に事前確率、なんですが…。
あれ。どの壺が選ばれる確率、というのは、条件としてどこにも出てきませんね。
そこで、ここでは各壺が選ばれる確率を勝手に仮定します。例えば、

P(壺1)…壺1が選ばれる確率＝1/3。
P(壺2)…壺2が選ばれる確率＝1/3。
P(壺3)…壺3が選ばれる確率＝1/3。

うん、これ、一番ストレートで直感的で、合ってるような気がしますよね。ということでこれで行きましょう*2。

さぁ、これで全て揃いました。事後確率を求めます。

P(壺3|赤) = 1×1/3 / (1/3×1/3 + 2/3*1/3 + 3/3*1/3) = 1/2

ということで答えは1/2。

この例題は、ベイズの公式（ベイズの展開公式）の使用例であるとともに、ベイズの公式の意味付けの理解のための例示にもなっています。

H_i…どの壺を選んだかという「原因」
D…どの色の玉が取り出されたかという「結果」

例示は理解の試金石、ですねっ（テトラちゃん風に）

他の例（取り敢えずトピック上げだけ）

今のは本当に分かりやすい（取っつきやすい）例ですが、他に

ベイズ更新：計算した事後確率を事前確率に設定してベイズ計算を繰り返す
モデルの設定：先ほど「事前確率は同等」と仮定したところを、事前確率をある分布に従うモノと仮定してベイズ計算

などなどありますが、今日はここまで。

次回

ナイーブベイズフィルターを実装したい。

*1:書籍「史上最強図解これならわかる！ベイズ統計学」(ISBN978-4-8163-5181-5)第3章§10の例題より。

*2:このように、「何も情報が無ければ確率は同等と仮定する」ことを「理由不十分の原則」と呼ぶそうです。