Haskell模写修行（2）

問題: Codeforces Beta Round #84 (Div. 2 Only) B. Lucky String

自分が書いたコード:

import Control.Monad

main = do n <- liftM read getLine
          putStrLn . take n . cycle $ "abcd"

お手本コード（bjin氏参考）:

main = interact$(`take`cycle"abcd").read

interact使って書くとごちゃごちゃしそうだったので、doブロックで書いた。Control.MonadのliftMは便利。適用する関数の引数の数に応じて、liftM2, liftM3, liftM4, liftM5 などもある。

bjinさんはゴルフ好きな感じである。takeをinfix operationにしてセクション（二項演算子の部分適用をするための手法）を実現してるのは面白いテクニックだなーと。

Haskell模写修行（1）

"New Zealand Autumn" by Abaconda

　

　モナドを実装したり、メモ化したり、コンビネータを書いたりとかしたけど、Haskellらしいコードが全然書けない。一応は動くものは書けるけど、どうも命令型で関数型っぽくない。ということで、Haskellらしいコードが書けるように修行する。

　"らしい"コードを書けるようになるためには、素晴らしいプログラマが書いた"らしい"コードをたくさん読んで吸収していくのが一番じゃないかなと思う。ということで、しばらく以下をやっていこうと思います。

1. codeforceの問題をHaskellで解く
2. ratingの高いHaskellerが書いたコードと自分のコードを比べる
3. いいところを吸収する

問題: Codeforces Beta Round #89 (Div. 2) A. String Task

自分が書いたコード:

import Char

isVowel 'a' = True
isVowel 'i' = True
isVowel 'u' = True
isVowel 'e' = True
isVowel 'o' = True
isVowel 'y' = True
isVowel 'A' = True
isVowel 'I' = True
isVowel 'U' = True
isVowel 'E' = True
isVowel 'O' = True
isVowel 'Y' = True
isVowel x = False

addDot   [] = []
addDot (x:xs) = '.' : x : addDot xs

main = do 
         input <- getLine
         putStrLn . map toLower . addDot . filter (not . isVowel ) $ input
         

お手本コード（bjin氏参考）:

import Char
g x | x < '0' || elem x "aoyeui" = ""
    | 1 > 0           = '.':[x]
main = interact $ concatMap (g.toLower)

ポイント:

1. elem
2. concatMap
3. interact

　”Haskellらしい”うんぬんより、まず便利な関数を覚えていくところからかな。

UNIX環境でサービスのポート番号を調べる方法

教えてもらったので忘れないうちにメモ。

例えば、MySQLのポートを調べたいときは、以下のようにします。

$ sudo lsof -P | grep mysql | grep TCP

　lsofはプロセスによって開かれているファイルの情報を出力するらしいです。ここでいう"ファイル"とは、通常のファイルやディレクトリ、さらにはライブラリ、ストリーム、ネットワークファイルまで含まれます。

-Pオプションは、ポート番号をポートの名前（サービス名）に変換するのを抑制するという意味です。-Pオプションを付けない場合は

TCP *:mysql (LISTEN)

と表示されるのが、-Pをつけると、

TCP *:1234 (LISTEN)

のようにポートの番号で表示されます。

ちょっと調べてみたところ、-iオプションというものもあって、これを使うとネットワーク関連の情報だけを取り出せます。上のコマンドは以下のようにしてもいいです。

$ sudo lsof -i tcp | grep mysql

　lsofの他にも、socklistというコマンドがあるみたい。このコマンドはオプションなしで

$ sudo socklist

とすると、lsofと同様に通信プロトコル、ポート番号、プロセスID、サービス名などが表示されます。こっちの方が楽ですねー。

MCMCで二次元正規分布に従う乱数を生成する

"The Borneo office" by angela7dreams

　マルコフ連鎖モンテカルロシミュレーション（MCMC）を用いて、事象(x, y)の発生確率が二次元正規分布

に従う事象をサンプリングしてみます。この分布はグラフに書くと下のようになります。

さて、この確率分布に従う乱数を発生させたいと思います。今回は一番簡単なメトロポリス法を用いて推移確率を決定します。

Haskellで書きました。自分の実力足らずでいろいろとやばいです。IOモナドのところが何かぎこちないです。実行速度も遅いです。Haskellもっと勉強しないとっ。。

import System.Random
import Text.Printf

select :: (Double, Double) -> IO (Double, Double)
select (x, y) = do
                    dx <- getStdRandom $ randomR (-5,5)
                    dy <- getStdRandom $ randomR (-5,5)
                    return (x+dx, y+dy)

prob :: (Double, Double) -> Double
prob (x, y) = 1/(2*pi) * exp(-(x*x+y*y)/2)

nextStatus :: (Double, Double) -> IO (Double, Double)
nextStatus (x, y) = do
                 (xt, yt) <- select(x, y)
                 let p = min 1  (prob(xt, yt) / prob(x, y))
                 q <- getStdRandom $ randomR (0, 1)
                 if p >= q then return (xt, yt) else return (x, y)

markovChain :: (Double, Double) -> [IO (Double, Double)]
markovChain (x, y) = iterate (>>= nextStatus) (return (x, y))

printPoints :: (Double, Double) -> IO()
printPoints (x, y) = printf "%f %f\n" x y

main = mapM_ (>>= printPoints) . take 3000 . drop 1000 . markovChain $ (0, 0)

それっぽい形になっています。E(x) = -0.00368、E(y) = -0.0121、V(x) = 1.03、V(y) = 1.02でした。試行回数を増やせば、それぞれ0, 0, 1, 1に近づいていくでしょう。

MCMC入門まとめ

マルコフ連鎖モンテカルロシミュレーション（MCMC）の入門資料をいくつか読んだので、分かったことと考えたことをメモしておく。

［分かったこと］

■アルゴリズムの流れ

1. 不定分布が分かっているものから、マルコフ性を仮定して遷移確率を推定する。

2. 適当な初期条件を選んで、1.の推移確率を用いて推移させる。

3. 十分な回数ループしたら終了。そうでなければ、2.へ戻る。

　1. では、ある特定の不定分布におちつく遷移確率行列は複数あることに注意。推移確率の選定には｢詳細釣り合い｣の条件（状態遷移図を書いたときに、確率のフローのようなものが釣り合っていると考える）が一般的に使われる。詳細釣り合いだけでは、自由度が残り値がfixしないため以下いずれかの方法で値を決める。

• メトロポリス法
• ギブスサンプラー
• メトロポリスヘイスティング（M-Hアルゴリズム）

　2. では、いろいろな初期解に対するある時刻における各事象の発生確率（パターン数を無限にする）と、ある初期解について無限回推移させた場合の各事象の発生確率が等しくなる（エルゴード性）という特性が活かされている。

■用途

1. 任意の分布に従うデータをサンプルできる
2. モンテカルロ積分

1. は、最初よく分からなかったが、良く考えると任意の確率分布にしたがう乱数をつくるのは結構むずかしい（[参考文献3] 参照）。普段よく利用するのは一様分布にしたがう乱数だが、正規分布N(0, 1)にしたがう乱数とか、もっとわけの分からない適当なp(x)にしたがう乱数発生器を作れて言われると難しい。MCMCを用いると、簡単に任意の乱数発生器を作成できる。

2.は、1.の応用のように思える。例えば、事象xが起こる確率をp(x)としてf(x)の期待値を求めよ。という問題を考える。∫ f(x) p(x) d(x)を計算することになるが、xの次元が大きいと計算が困難なためMCMCを利用したモンテカルロ積分がしばしば利用される。

［考えたこと］

　MCMCはマルコフ性のある事象に対するシミュレーション手法だと勘違いしていたが、どちらかというと任意の確率分布にしたがうサンプリングデータを生成できるということに利用価値がありそうだ。

　MCMCでマルコフ性のある事象をシミュレーションする場合は、確率遷移行列の選び方がとても重要だと思う。パラメータの取り方に自由度があるメトロポリスヘイスティングが一番使えそうだが、提案密度の選び方にはどうしても主観が入ってしまいそう。MCMCはベイズ識別にも応用されているらしいので、そのあたりを詳しく調べる必要がありそうだ。（というかそこがそもそも知りたかったところなので。）

[参考文献]

1. マルコフ連鎖モンテカルロ法入門-1
2. マルコフ連鎖モンテカルロ法入門-2
3. Markov chain Monte Carlo

確率遷移行列の形と定常状態

マルコフ性（単純マルコフ過程）をもつ事象がある確率推移のもと、どのような定常状態に落ち着くのか実験してみた。

　Haskellで書きました。こういうのはHaskellが綺麗に書けますねー。

genNextPrb a b c d = \(x, y) -> (a*x+b*y, c*x+d*y)
main = mapM_ print . take 200 . iterate nextPrb $ (1, 0)
       where nextPrb = genNextPrb 0.2 0.2 0.8 0.8

分かったこと。確率遷移行列をP, 事象の分布の初期値をx0, 収束時の事象の分布をxnとおきます。

1. Pが単位行列の場合は、xn = x0となる。
2. すべての行iにおいてPi0 = Pi1 = Pi2 =...の場合（列の値が同じ）、xn = [Pi0, Pi1, Pi2, ...]となる。
3. Pが対象行列の場合、xnはすべての要素の値は等しい。（一様分布になる）
4. 1.以外の場合は、xnはx0の値によらない。（エルゴード性）

　少しだけ補足。

　1.はベクトルに単位行列かけるだけなので確率分布がかわらないのは当たり前。

　2.は、次に特定の事象がおこる確率は、現在の状況に依存しないことを表している。例えば、サイコロをふる場合など。これはマルコフ性を持つ事象とは呼ばないのかもしれない。

　3.はおもしろい。x ∈ R^2のときは感覚的に分かりやすい。

例えば、天気の推移にマルコフ性があると仮定して

• 晴れの次が晴れの確率を0.9
• 晴れの次が雨の確率を0.1
• 雨の次が晴れの確率0.1
• 雨の次が雨の確率0.9

とする。このとき確率推移行列は対称行列である。一旦晴れになるとなるべくその状態を保とうとするし、一旦雨になるとなるべくその状態を保とうとする。よって仮にx0 = [1, 0]、つまり晴れの確率100%などとしても、次の日には雨が10%になって、その10%は晴れ側にほとんど戻ってこないのでそれが蓄積されていって最終的にxn = [0.5, 0.5]で落ち着く。これは感覚的に分かりやすい。

　では一般にx ∈ R^nの場合はどうか？これは対称行列を写像として見たときの意味を考えるか、対角化して極限を取るみたいなことをして数学的に証明できるのかもしれないが、大変そうなので直感的に理解する方法を考えてみる。オートマトンみたいな状態遷移図を書けば分かりやすそう。。3次元の場合も、以下の図をみると直感的には定常状態では均等に落ち着きそうです。

たくさん持っているところからはたくさん出て行き、少ししかないところからは少ししかでていかない。同じ色の線は同じ遷移確率を持るのがポイントです。（電子回路のオペアンプをバッファとして使う場合に定常状態がどうなるかを考えるのに似ているかも。）

4. は、確率遷移行列が単位行列以外の場合は、定常状態が初期状態によらないことを意味している。初期状態が何でもいいというのは、シミュレーションをする場合にとてもありがたい特性である。また、初期状態によらないということは、マルコフ性を持つ確率モデルは遷移行列によってのみ特徴づけられることを意味している。固有値とか固有ベクトルとかがここでも使えそうである。

sedの最短一致

最近知ったけど、sedはstream editorの略らしい。

それはどうでもよくて、sedの正規表現のマッチングの小技的なものをメモとして書いておく。

HTMLをパースする場合を考える。

　

　例えば、

<tr><td>id</td><td>data</td></tr><tr><td>1</td><td>AKB48</td></tr><tr><td>2</td><td>perfume</td></tr><tr><td>3</td><td>bump.y</td></tr>
…

のようなデータからHTMLのタグだけを消して、意味のある値だけを取り出したいとする。つまり、

id data

1 AKB48

2 perfume

3 bump.y

…

のような出力をえたい。

　このときに、

sed s/\<.*\>//g test.html

みたいな書き方をすると悲しいことになる。（行全体にマッチしてしまい、すべて空行になってしまう。）

これは、複数のマッチングが考えられる場合なるべく長いものにマッチする（最長一致）からである。この性質はsedに限ったことではなくて正規表現全般に言えることで、デフォルトでは最長一致である。

　最短一致を利用するためのオプションを持つ言語もあるが、sedには無いので何らかの工夫が必要となる。今回の場合は、閉括弧が現れた時点で一旦マッチングをやめるようにすればいいので、以下のように書くと思惑のデータがえられる。

sed s/\<[^\>]*\>//g test.html

　データの間にスペースを入れたい場合は、

sed 's/<[^>]*>/ /g' test.dat

とする。（スペースを用いるので''で第一引数を囲む。''を使うとエスケープ文字の扱いが変わるので<>はエスケープしない。※注）

※注)

　シングルクオーテーションで囲んだ場合、シングルクオーテーション以外の文字は普通の文字と同じ扱い。ダブルクオーテーションで囲んだ場合は、"$`\以外の文字は普通の文字と同じ扱いになるそうです。

http://genki01.cc.hokudai.ac.jp/reo/pukiwiki/index.php?shellscripts