Wolfeyes Bioinformatics beta

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

例 6.1 メトロポリス・ヘイスティング・アルゴリズムでベータ分布の乱数を生成する

例2.7では受理・棄却法を用いてベータ分布の乱数を生成したが，本問題ではメトロポリス・ヘイスティング・アルゴリズムを使ってベータ分布からサンプリングを行うことを目的とする．

a <- 2.7
b <- 6.3
c <- 2.669
Nsim <- 5000
X <- rep(runif(1), Nsim)

for(i in 2:Nsim){
  Y <- runif(1)
  rho <- dbeta(Y, a, b)/dbeta(X[i-1], a, b)
  X[i] <- X[i-1] + (Y - X[i-1])*(runif(1)<rho)
}

plot(1:Nsim, X, type="l", xlab="Iterations", main="Sequence X^t")
plot(4500:4800, X[4500:4800], type="l", xlab="Iterations", main="Sequence X^t for t = 4500 - 4800")

まず，以下の図は系列 ${X^{(t)}}$ をプロットしたものである．

そして，その中からシミュレーション回数が4500〜4800回までの系列を抜き出したのが以下の図となっている．このように一部分を拡大して見てみると，系列を表す線が変化せずにx軸に並行に遷移している部分が見られる．これは $\rho(x^{(t)},Y_t)$ の確率で棄却されて値が変化していないことを表している．これは対応する $Y_t$ に依存し，マルコフ連鎖の性質から次の状態にも影響するため，値が変化しない状態が一定期間続くことがあるからである．

また，得られた系列がベータ分布の乱数となっているかを確かめるため，左側に今回のメトロポリス・ヘイスティング・アルゴリズムで作成した乱数，右側にRの組込み関数rbeta()で生成した乱数をヒストグラムで図示し，今回実験で設定したベータ分布 $\mathcal{Be}(2.7,6.3)$ を赤の曲線で表した．

par(mfrow=c(1,2))
hist(X, nclass=100, freq=F, xlim=c(0,1.0), main="Metropolis-Hastings")
curve(dbeta(x,a,b), add=T, col="red", lwd=2)
hist(rbeta(Nsim,a,b), freq=F, nclass=100, xlim=c(0,1.0), main="Direct Generation")
curve(dbeta(x,a,b), add=T, col="red", lwd=2)

この図を見ると，メトロポリス・ヘイスティング・アルゴリズムを使った場合でも，受理・棄却法の時と同じく乱数の生成に成功していることがわかる．

例 6.2 メトロポリス・ヘイスティング・アルゴリズムでコーシー分布の乱数を生成する

例6.1のベータ分布の乱数生成における提案分布は $\mathcal{U}_{0,1}$ だったが，コーシー分布の乱数生成における提案分布には，

標準正規分布： $\mathcal{N}(0,1)$
自由度0.5のスチューデントt分布： $\mathcal{T}_{1/2}$

を使うことができる．このような提案分布は独立提案分布であり，マルコフ連鎖ではない．

なお，本問題ではこれらの分布を表す単語として候補分布と提案分布という2つの言葉が出てくるが，今回は表記として提案分布で統一する．

# Standard Normal Distributionp
Nsim <- 10^4
X <- c(rt(1,1))
for(t in 2:Nsim){
  Y <- rnorm(1)
  rho <- min(1, dt(Y, 1)*dnorm(X[t-1]) / (dt(X[t-1],1)*dnorm(Y)))
  X[t] <- X[t-1] + (Y - X[t-1])*(runif(1) < rho)
}

# t-distribution (0.5 degree of freedom)
X2 <- c(rt(1,1))
for(t in 2:Nsim){
  Y <- rt(1, 0.5)
  rho <- min(1, dt(Y, 1)*dt(X2[t-1], 0.5) / (dt(X2[t-1],1)*dt(Y, 0.5)))
  X2[t] <- X2[t-1] + (Y - X2[t-1])*(runif(1) < rho)
}

まずはそれぞれの系列をプロットして特徴を見ていく．左側に提案分布を標準正規分布とした場合，右側に提案分布をスチューデントt分布とした場合で，上から順に系列のプロット，系列のヒストグラム，系列の自己相関プロットを図示したのが以下の図になっている．

まずは標準正規分布を使った場合からみていく．系列のプロットは例6.1と同じような形を示しているが，イテレーションが8500を過ぎたあたりで値が一定期間変化しない部分がある．これは乱数から生成する正規候補 $Y_t$ が大きい（または小さい）値を取ってしまい，連鎖で長い間引きずることになったからである．この現象は，系列のヒストグラムの右端において値が少し大きくなっていることからも確認することができる．

一方，スチューデントt分布を使った場合においては，系列のプロットは0付近で集中しており，たまに正規候補 $Y_t$ が非常に大きい値を取ることはあるものの，それが継続して続くわけではないことがわかる．系列のヒストグラムを見ても，標準正規分布の場合と比べて分布がなだらかでコーシー分布に当てはまっていることがわかる．ただし，実際には両端の外れ値を除外して系列のヒストグラムをプロットしているので，このヒストグラムの描写範囲から外れている値があることには注意が必要である．

次に，2つの提案分布におけるメトロポリス・ヘイスティング・アルゴリズムによって生成したコーシー乱数の累積収束プロット，すなわち乱数の平均値がどのように収束していくかを図示する．

黒線が提案分布に標準正規分布を使った場合，黄色の線が提案分布にスチューデントt分布を使った場合となっている．今回シミュレートしているコーシー分布から求めた正確な値は0.896なので，スチューデントt分布はイテレーション回数を重ねるごとに真の値に収束していくことがわかる．一方で標準正規分布においては，今回のイテレーション回数(N=5000)では真の値に収束していないようにみえる．しかし，平均値の傾向を見ると値が徐々に上昇しては一気に下がり…といった傾向が見られ，理論的にも標準正規分布を提案分布として用いた場合でも真の値に収束することがわかっている．なので，今回の場合はイテレーション回数が少なすぎるために，収束の様子が見られない．

Sublime Text 2でOctopressの基本コマンドを実行できるパッケージを使ってみた．

http://blueplanet.github.com/blog/2012/10/16/sublime-text-2-octopress/

1. sublime-text-2-octopressをインストールする

Package Controlでsublime-text-2-octopressをインストールする．これは定石通り．

2. Octopress.sublime-settingsを編集する

sublime-text-2-octopressの環境設定を行う．Sublime Text 2のメニューバーから”Sublime Text 2”→”Preferences”→”Package Settings”→”Octopress”→”Settings-User”でOctopress.sublime-settingsを開く．

今回設定するのは”octopress_path”と”octopress_cmd_before_rake”の2箇所．”octopress_path”はOctopressのディレクトリを指定する．そして”octopress_cmd_before_rake”にはrubyやrakeのパスを通す設定を読み込むコマンドを指定するのだが，下手に.zshrcなどを読み込むと不要な設定が多すぎてエラーになってしまうので注意が必要になる．自分の環境の場合，RubyはRVMで管理しているため，”octopress_cmd_before_rake”にはRVMの設定ファイルである“$HOME/rvm/scripts/rvm”を読み込むように指定する．

{
  // path to your octopress
  "octopress_path": "/Users/yag_ays/Dropbox/octopress",  
  // command to run before calling rake, eg source ~/.bash_profile to set up your local environment inc paths to ruby, rake etc.
  "octopress_cmd_before_rake" : "source /Users/yag_ays/.rvm/scripts/rvm",
  // set to generate, deploy or generate_and_deploy if you wish to have your changes generated into the /public folder and/or deployed upon file save
  "octopress_onsave_action": "",
  // true or false
  "use_bundle": false
}

3. OctopressのNew Postを実行してみる

それでは試しに新しい記事を書いてみる．Command Paletteで「octopress」と入力すると，sublime-text-2-octopressでサポートされているコマンド一覧が表示される．

「Octopress: New Post」を選択すると，画面下にタイトルを入力する欄が表示される．

タイトルを入力してEnterを押すと，Octopressのテンプレートが表示され，New Postが実行される．

この他にも，既存のPostやPageを編集したり，GenerateやDeployができる仕様となっている．

参考

最近ローカルの開発環境をEmacsからSublime Text 2に乗り換えたので，長い間.emacsに溜め込んだ環境をSublime Text 2でも再現するために，ちまちまと移行作業をしている．Sublime Text 2の強みというのはPythonで拡張パッケージが書ける上にパッケージ管理が非常に簡単というところにあると思っていて，とりあえずデフォルトの環境をなんとかしたいと思えばパッケージを探せば何とかなる．

ということで，今回はEmacsで構築していたClojureの開発環境を，SublimeREPLパッケージを使ってSublime Text 2で再現してみる．Emacsではswank-clojureとslime-replで書いたコードを逐次実行できる環境を作っていたが，今回もSublime Text 2内でClojureを動かして，エディタ内でコードの編集と実行を両立できるようにする．なお，今回の環境構築ではClojureやLeiningenのインストールはひと通り終えているものとする．

ちなみに，私の過去のEmacs＆Clojure環境は以下のblogにまとめてある．

Clojure + Emacs環境のセッティング（MacOSX Lion 10.7） - やぐブロ

1. Package Controlをインストールする

これについては他のサイトでもよく解説されていることなので，簡単に紹介とインストール方法だけ．Sublime Text 2ではSublime Package Controlというパッケージ管理が定番となっていて，Sublime Text 2内のコマンドから検索・インストール・最新バージョンへの更新などを簡単に行うことができる．

まず，Sublime Text 2のコンソールに以下のコマンドを入力（コピペ）する．

import urllib2,os;pf='Package Control.sublime-package';ipp=sublime.installed_packages_path();os.makedirs(ipp) if not os.path.exists(ipp) else None;open(os.path.join(ipp,pf),'wb').write(urllib2.urlopen('http://sublime.wbond.net/'+pf.replace(' ','%20')).read())

インストールが終了したら，Sublime Text 2を再起動する．メニューバーの”Tools”からCommand Paletteを開いて「package」と入力して，以下のようにPackage Controlのコマンドが表示されればOK．これでパッケージのインストールが簡単に行えるようになる．

2. SublimeREPLをインストールする

次に，Clojureをエディタ内で動かすためにSublimeREPLというパッケージをインストールする．ちなみにこのパッケージはClojure専用というわけではないので，RubyやPython，Rなどの他の言語のREPLの機能も兼ね備えている．

まず，先ほどの確認のときのようにCommand Paletteを開いて「install」と入力し，「Package Control: Install Package」を選択する．

すると，登録されているパッケージ一覧が表示されるので，「SublimeREPL」と入力してパッケージを選択してインストールを行う．

これで，SublimeREPLパッケージのインストールは完了．メニューバーのToolsにSublimeREPLという項目が表示されればOK．

3. SublimeREPLの環境設定にLeiningenのPATHを追記する

最後に，SublimeREPLの環境設定を行う．具体的にはLeiningenがインストールされているディレクトリにPATHを通すのだが，これをしないと「OSError(2, ‘No such file or directory’)」といったエラーが表示されるので注意．

今回の場合は，パッケージ内の設定ファイルに直接追記することにする．まず，SublimeREPLのパッケージのディレクトリ内にあるClojureのMain.sublime-menuというファイルを開く．これには，Sublime Text 2の”Preference”→”Browse Packages…”でFinderを表示してディレクトリを辿っても良いが，今回はTerminal.appで以下のコマンドを入力する．

$ open $HOME/Library/Application\ Support/Sublime\ Text\ 2/Packages/SublimeREPL/config/Clojure/

そして，Finderで表示された中にあるMain.sublime-menuをエディタで読み込み，29行目のextend_envの行にLeiningenがあるディレクトリのPATHを入力する．今回はMacPortsで入れたLeiningenを使用するために，/opt/local/binを指定する．

"extend_env": {"INSIDE_EMACS": "1" }

と書かれている行の末尾に「”PATH”: “{PATH}:/opt/local/bin/”」を追記する．

"extend_env": {"INSIDE_EMACS": "1", "PATH": "{PATH}:/opt/local/bin/"}

4. SublimeREPLの起動

メニューバーの”Tools”から”SublimeREPL”→”Clojure”→”Clojure”を選択して，以下のようにREPLが起動すれば，無事にSublimeREPLの環境設定は終了．

環境

OS
- Mac OS X 10.8.2
Sublime Text 2
- Version 2.0.1, Build 2217
Leiningen
- Leiningen 1.6.1.1

参考

相変わらずデザイナーの実験場のような雑誌だ．ページレイアウトを重視するあまり規定ページ数から溢れた記事の一部を巻末に押しやるという普通なら考えられないことをするのは，この雑誌だけなんじゃないかと思う．そんな雑誌の作り方をしていることもあってか，ページをペラペラめくりながら目に留まった記事を片っ端から読んでいくという読み方も楽しいし，やはり1冊の雑誌として興味あるもの無いものを混ぜこぜにして読めるのは，紙媒体だからこそできる体験だろう．個人的には，なんだかんだ言って雑誌は紙の方が好きだ．机の上に置いてあるだけで何度も開いて読みたくなる．

本号の特集は「教育」ということで，オンラインで大学の授業を受けられるUdacity・Courseraなどのプログラムや，TEDやMITメディアラボなどの取り組みも紹介される．この特集の他には，Twitter創業者の一人のジャック・ドーシーのインタビューやスタートアップの紹介，宇宙やSF，面白ガジェットなどなど，普段WIRED.jpで見かけるようなタイプの記事を，多彩な写真・画像とボリュームある記事で読むことができる．個人的には伊藤計劃などのSFまわりの海外展開やマルチメディア展開が興味深かった．こういった裏で着々と進行していくプロジェクトなんかは，ニュース記事などで断片がポツポツと出てくるくらいで忘れた頃にまた続報がといったことになることが多いので，こうやってひとつの記事として出てきてくれると記憶にも残りやすい．

そういえば，最近定期購読している雑誌はWIREDくらいだ．去年に新装してから買い続けていることになる．前々からよく読んでいるWEB+DB PRESSとSoftware Designは気になる記事があれば買う程度なので，不定期ながらも情報源やネタ帳として軽い感じで楽しめるWIREDはある意味貴重な存在になっている．

前回の「佐藤可士和の超整理術」を読んでる途中にコンビニで整理術の雑誌を見つけたので読んでみた．相変わらずAssocieは読者ターゲットがブレない記事構成でいつも感心する．自分の手の届く範囲でなんとかやりくりしないといけない新人ビジネスマンも，なかなか大変なご様子．

内容としては，やはり机で働くビジネスマン向けの実践的なテクニックが多め．机の上を片付けたり電子化したり名刺を整理したりと王道は抑えつつ，タクシー運転手やカメラマン，料理人などの少し変わった整理術も紹介される．整理といえばやはりバインダーや手帳などの文房具に自然と話が進むわけで，定番文房具の使い方から最新文房具の紹介までひと通りカバーされているのもいつも通り．

そういった内容がいろんな分野の中堅社員の具体例とともに語られるのだが，前回の佐藤可士和の本を読んだ後に眺めると，それぞれのテクニックは一体何がしたいのかということがクリアに理解できて面白かった．佐藤可士和のように何もモノを置かないというのは現実的に厳しいだろうし，本の中ではあまり語られなかった具体的な方法をこういった雑誌で穴埋めできると思う．ポイントとしては，どうプライオリティをつけて取捨選択するのか，そしていかに整理のための規格を統一するかの2点に集約されるというあたりで，あとは自分に合った方法で実践するのみ．

ちなみに，個人的には針無しホッチキスに非常に興味を惹かれたが，わざわざ買うものだろうかと考えた結果購入は見送りになった．針無しだとシュレッダーにそのまま入れられるメリットはあるが，キチンと書類を留めてくれる方が大事ということで，今のところは物欲を抑えている．ちなみに，ホッチキスの針を簡単に取ることができる「サンスター文具はりトルプロ 4765-800」は本当に便利なのでオススメ．ハサミを握る程度の軽い力で簡単に針を抜くことができるので，量が多くても全然疲れない．これのお陰で，針を取る手間をそれほど煩わしく感じていないというのもある．一方，「リムーバーブルー」は，針を抜くときに力をいれなければいけないのでちょっと大変．それに分厚い紙を束ねている針はほとんど取れない．