Wolfeyes Bioinformatics beta

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

内容：一様分布を元にして様々な確率分布に従う乱数を生成する

統計解析ソフトであるRには，様々な分布に対応した組み込み関数が用意されている．本章では，そういった組み込み関数を使わずに，一様分布から生成される乱数を逆関数で変換することで，他の確率分布の乱数を表現する．

本書で扱う乱数とは，完全なランダム性を持つ乱数ではなく擬似乱数である．擬似乱数はset.seed()関数で設定した値を種として乱数を生成するため，どのような環境においてもset.seed()関数で同じ値を使うことで乱数を再現することができる．

今回は，一様乱数を元にして指数乱数を生成し，指数乱数からガンマ分布やベータ分布の乱数へと変換していく．これらの確率分布は全てRの関数で用意されているものなので，一様乱数から生成した乱数とRの関数から生成した乱数を比較することによって，乱数の生成がうまくいっているかどうかを判断する．以下のコードで描写したヒストグラムはすべて，左側が一様分布から生成した乱数，右側がRの関数を用いて生成した乱数(N=10⁴ )．赤い曲線はどちらもRの関数を用いて分布を示したもの．

例 2.1 一様乱数から指数乱数を作成する

これはテキストにある通りのコードと作図．mcsmパッケージのdemo(Chapter.2)にも同様のコードがある．

以下のコードでは $-\log{(1-u)}$ をそのままコードに落とし込んでいる．テキストにもある通り， $U \sim \mathcal{U}_{[0,1]}$ ならば0から1の間で一様に等しい確率分布なのだがら， $U$ も $1-U$ も一様分布になる．

par(mfrow=c(1,2))
Nsim <- 10^4
U <- runif(Nsim)
X <- -log(1-U)
Y <- rexp(Nsim)
par(mfrow=c(1,2))
hist(X, freq=F, main="Exp from Uniform", nclass=50)
curve(dexp(x), add=T, col="red", lwd=2)
hist(Y, freq=F, main="Exp from R", nclass=50)
curve(dexp(x), add=T, col="red", lwd=2)

練習問題 2.2 逆変換法を用いてロジスティック分布とコーシー分布の乱数を生成する

逆変換法を用いてロジスティック分布とコーシー分布の乱数を生成する．変換で求まる関数はいわゆる逆累積分布関数というもの．

a.ロジスティック分布

$u = \frac{1}{1+e^{-\frac{x-\mu}{\beta}}}$ を変形して， $x = -\beta \log{\frac{1-u}{u}} + \mu$

以下のコードは $\mu=0$ ， $\beta=1$ の場合．

Nsim <- 10^4
U <- runif(Nsim)
par(mfrow=c(1,2))
X <- -log((1-U)/U)
Y <- rlogis(Nsim, 0, 1)
hist(X,freq=F,main="Logis from Uniform")
curve(dlogis(x), add=T, col="red", lwd=2)
hist(Y, freq=F, main="Logis from R")
curve(dlogis(x), add=T, col="red", lwd=2)

b.コーシー分布

$u = \frac{1}{2} + \frac{1}{\pi} \arctan{\frac{x-\mu}{\sigma}}$ を変形して， $x = \mu + \sigma \tan{\pi(u-\frac{1}{2})}$

以下のコードは $\mu=0$ ， $\sigma=1$ の場合．

Nsim <- 10^4
U <- runif(Nsim)
x <- tan(pi*(U-0.5))
Y <- rcauchy(Nsim)
hist(X, freq=F, main="Cauchy from Uniform", xlim=c(-10,10))
curve(dcauchy(x), add=T, col="red", lwd=2)
hist(Y[Y>=-10 & Y<=10], freq=F, main="Cauchy from R", xlim=c(-10,10))
curve(dcauchy(x), add=T, col="red", lwd=2)
# cauchyは値が両端に飛んでヒストグラムが綺麗に書けないので[-10,10]で整えている

練習問題 2.12 指数分布からガンマ分布とベータ分布の乱数を生成する

a.ガンマ乱数

以下の作図では， $\beta$ の値を固定して $a$ の値を1,2,5,9と変化させたときの分布の変化を見ている．

Nsim <- 10^4
beta <- 2.0
par(mfrow=c(4,2))
for(a in c(1,2,5,9)){
  U <- runif(a*Nsim)
  m <- matrix(U, nrow=a)
  X <- beta * apply(-log(m), 2, sum)
  Y <- rgamma(Nsim, shape=a, scale=beta)
  hist(X, freq=F, main=paste("Gamma from Uniform:","a=",a," beta=",beta), nclass=50)
  curve(dgamma(x, shape=a, scale=beta), add=T, col="red", lwd=2)
  hist(Y, freq=F, main=paste("Gamma from Uniform:","a=",a," beta=",beta), nclass=50)
  curve(dgamma(x, shape=a, scale=beta), add=T, col="red", lwd=2)
}

a.ベータ乱数

この手法では $(a \in \mathbb{N^*})$ ， $\mathbb{N^*} = 1,2, \dots$ という制約があるため， $a=0.1$ ， $b=0.1$ のようなベータ分布は作ることが出来ない．以下はa=1,b=1，a=2,b=3，a=8,b=4の場合．

# a=1,b=1の場合
par(mfrow=c(3,2))
Nsim <- 10^4
a <- 1
b <- 1
U <- runif((a+b)*Nsim)
m <- matrix(U,nrow=(a+b))
X <- -log(m[a,]) / apply(-log(m),2,sum) # a>=2の場合はapply(-log(m[1:a,]),2,sum) / apply(-log(m),2,sum)
Y <- rbeta(Nsim,a,b)
hist(X, freq=F, main=paste("Beta from Uniform:","a=",a,"b=",b), nclass=50, ylim=c(0,3))
curve(dbeta(x,a,b), add=T, col="red", lwd=2)
hist(Y, freq=F, main=paste("Beta from R:","a=",a,"b=",b), nclass=50, ylim=c(0,3))
curve(dbeta(x,a,b), add=T, col="red", lwd=2)

b. 一様分布から指数乱数を作る(逆変換)

$U = e^{-\lambda x}$ より $x = - \frac{\log{U}}{\lambda}$

c. 一様分布からロジスティック乱数を作る(逆変換)

$U = \frac{e^x}{1+e^x}$ より $x = \log{\frac{U}{1-U}}$

「Rによるモンテカルロ法入門」

本書は，モンテカルロ法の実践的な解説書であり，統計解析ソフトのRを用いた豊富な実例と練習問題が組まれている．モンテカルロ法とは乱数を用いて数値計算を行う手法の総称であり，本書で扱う内容は乱数の発生からモンテカルロ積分，そしてマルコフ連鎖モンテカルロ法（MCMC）の各種アルゴリズムに至るまで非常に幅広い．たいていの解説には理論に実践演習が付随した形となっており，数学的な理論を軸にして実際にRを用いたコード例が示される．

練習問題を解きつつ読書ノートをまとめてみる

そんなこんなで，久保本と並行する形で「Rによるモンテカルロ法入門」を読んでいる．一応MCMCの部分だけひと通り目を通したのだが，最終的にMCMCの実装までひと通りやるにしても一連の流れを簡単にでも追っておかなければと思って，最初の乱数の部分からじっくり読み進めている．これがなかなか難しくて，手も足も出ないところをなんとかRのコードを動かして理解を補っているわけなのだが……取り敢えず練習問題を解ける範囲で解きつつRのコードを書いていこうというこでとはセクションごとに小分けてして，じっくり読書ノートを付けて勉強することにしてみた．

以下に，読書ノートのリンクを書き加えていく

読書ノートまとめ

「Rによるモンテカルロ法入門」読書ノート：2章逆変換・一般変換法
- 例 2.1，練習問題 2.2，練習問題 2.12
「Rによるモンテカルロ法入門」読書ノート：2章受理・棄却法
- 例 2.7，例 2.8，練習問題 2.7，練習問題 2.8
「Rによるモンテカルロ法入門」読書ノート：2章練習問題
- 練習問題 2.18，練習問題 2.19，練習問題 2.22
「Rによるモンテカルロ法入門」読書ノート：3章モンテカルロ積分
- 例 3.3，例 3.4
「Rによるモンテカルロ法入門」読書ノート：3章モンテカルロ積分その2
- 練習問題 3.1，練習問題 3.9
「Rによるモンテカルロ法入門」読書ノート：3章重点サンプリング
- 例 3.5，練習問題 3.5
「Rによるモンテカルロ法入門」読書ノート：5章モンテカルロ最適化数値的最適化法
- 例 5.1，例 5.2
「Rによるモンテカルロ法入門」読書ノート：5章モンテカルロ最適化確率的探索
- 例 5.3，練習問題 5.3
「Rによるモンテカルロ法入門」読書ノート：5章モンテカルロ最適化確率的探索その2
- 例 5.4，例 5.6
「Rによるモンテカルロ法入門」読書ノート：5章モンテカルロ最適化確率的勾配法
- 例 5.7
「Rによるモンテカルロ法入門」読書ノート：6章マルコフ連鎖とランダム・ウォーク
- 練習問題 6.1，練習問題 6.2
「Rによるモンテカルロ法入門」読書ノート：6章メトロポリス・ヘイスティング・アルゴリズムその1
- 例 6.1，例 6.2
「Rによるモンテカルロ法入門」読書ノート：6章メトロポリス・ヘイスティング・アルゴリズムその2
- 例 6.3，練習問題 6.11
「Rによるモンテカルロ法入門」読書ノート：6章ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズム
- 例 6.4，例 6.10，練習問題 6.14
「Rによるモンテカルロ法入門」読書ノート：6章メトロポリス・ヘイスティング・アルゴリズムその3
- 練習問題 6.13
「Rによるモンテカルロ法入門」読書ノート：6章メトロポリス・ヘイスティング・アルゴリズムを用いたモデル選択
- 例 6.8，練習問題 6.5

私的正誤表

全体的な修正点
- Rのoptimize()関数で最大値を求める際には「max=T」→「maximum=T」を使う．
P.143 例題5.1内のコード8行目
- 「maximum=max=T」→「maximum=T」
P.155 例5.7内のコード4行目
- 「diff > 10^-5 」→「dif > 10^-5 」
- ただしこれは原書でも同様の表記をしている
- mcsmパッケージのデモスクリプトでは直っている
P.171 式(5.12)右辺
- 「 $\max_{\theta}Q(\theta|\hat{\theta}_{(j-1)}\bf{x})$ 」→「 $\max_{\theta}Q(\theta|\hat{\theta}_{(j-1)},\bf{x})$ 」
P.199 下から10行目のplot()の中
- 「nsim」→「Nsim」
P.200 下から4行目
- 「 $\varepsilon_{ij} \sim N(0,\sigma^2)$ 」→「 $\varepsilon_{ij} \sim \mathcal{N}(0,\sigma^2)$ 」
- これは原書でも同様の表記をしている
- ここでの $\mathcal{N}$ は正規分布， $N$ はサンプルの総数

実行環境

R version 2.15.1 (2012-06-22)

お決まりごと

この読書ノート（上記リンク先を含む）は個人的なメモであって，解答およびコードの正確さや厳密さを保証するものではありません．間違いが多分に含まれている可能性があるので，参考にされる際には十分ご注意下さい．

参考：サンプルコード

CRANのmcsmパッケージに，本書の一部コードのスクリプトがデモとして入っている．

http://cran.r-project.org/web/packages/mcsm/index.html

> install.packages("mcsm")
> library(mcsm)
> demo(Chapter.1) # Chapter.8まである

参考：基本的な確率分布の略称

一様分布: $\mathcal{U}$
指数分布: $\mathcal{Exp}$
カイ二乗分布: $\mathcal{\chi^2}$
ガンマ分布: $\mathcal{G}$
ベータ分布: $\mathcal{Be}$
正規分布: $\mathcal{N}$
二項分布: $\mathcal{Bin}$
負の二項分布: $\mathcal{Neg}$
ポアソン分布(※): $\mathcal{P}$
コーシー分布: $\mathcal{C}$
スチューデントのT分布: $\mathcal{T}$
二重指数分布（ラプラス分布）: $\mathcal{L}$
逆ガンマ分布 : $\mathcal{IG}$

（※ $\mathcal{P}$ は場合によって変わる）

参考

(pdf) http://www.ceremade.dauphine.fr/~xian/shortManual.pdf

公式ページには，番号が奇数の練習問題だけ解答がある．

(pdf) http://www.stat.ufl.edu/~casella/ShortCourse/MCMC-UseR.pdf

UseR!の発表スライドっぽい資料．

http://waxworksmath.com/Authors/N_Z/Robert_N_Casella/robert_n_casella.html

1章から3章途中までの解説とRのコードがある．

概要

Jekyll/Octopressで数式を使った記事の作成を簡単にします．

Google Chart Toolsには，Tex形式で書かれた数式をURLに含めて投げると変換した画像を返してくれるAPIがあります．しかし，URLで数式を投げるには特定の文字のエスケープ(変換)をする必要があり，とても面倒です．ということで，Jekyll/Octopress内に書かれた数式を自動でGoogle Chart APIのURL形式に変換するプラグインを作りました．これを使えば，記事の中で数式を書くときに，面倒なHTMLのタグやURLの変換を気にする必要がありません．

使い方

以下のようにgctexというタグとTex形式の数式を埋め込むと，Google Chart APIの画像のURLに変換されます．

{% gctex \frac{1}{\sqrt{2\pi}\sigma} \exp(\frac{-(x-\mu)^2}{2\sigma^2}) %}

また，文字の途中に埋め込むことも可能です．

これは{% gctex \mathcal{N}(x|\mu,\sigma) %}です．

これはです．

注意

ここで使用しているLiquid Tagプラグインは，中に改行を入れてしまうとコンパイルエラーを起こして正しく変換されません．なので，

{% gctex x \\ y \\ z %} 

TeX内で改行を表すことはできますが，

{% gctex x \\
         y \\
         z% } 

のように，実際にテキストの中で改行を含めるような書き方はできません．

インストールとコード

以下のgooglechart_tex.rbをダウンロードして，プラグインのコードがあるディレクトリに入れるだけです（Jekyllの場合は_plugins，Octopressの場合はplugins）．

# Convert tex formula into google chart API format (Jekyll/Octopress plugin) # # Usage: {% gctex tex_formula %} # # Author: yag_ays (yanagi.ayase [at] gmail.com) # blog http://yagays.github.com/ # # Ref : http://d.hatena.ne.jp/repose/20110506/1304612242 require "cgi" module Jekyll class GoogleChartTexTag < Liquid::Tag def initialize(tag_name, text, tokens) super @text = text end def render(context) "<img src=http://chart.apis.google.com/chart?cht=tx&chl=#{CGI.escape(@text)}>" end end end Liquid::Template.register_tag('gctex', Jekyll::GoogleChartTexTag)

参考

ウェブデザイン・インフォグラフィクスを含めたデザイン全般に使えるリファレンス的な1冊

一般の社会人や学生が「デザイン」を意識するのは，プレゼンテーションのスライドを作成するときだろう．エンジニアや研究者などはそれに加えて，ポスターやウェブサイトの作成，そして最近よく言われるようになったUX/UIなどの要素が加わる．そういう人に私が問答無用でお薦めするのが「ノンデザイナーズ・デザインブック [フルカラー新装増補版]」なのだが，本書「デザインする技術 ~よりよいデザインのための基礎知識」もなかなかの良書だったので紹介してみる．

本書は，デザイナーの仕事としての企画・発案から制作に至るすべての工程におけるデザインの技法を包括的に取り扱った技術書である．それらデザイナーの技法を，5つの要素「考」「図」「文字」「面」「色」に分けて，各要素が含む細かく技法をそれぞれ体系的に解説している．例えば「考」では，アイデアの出し方から組み立て方，そして最終的なカタチに仕上げていく過程に必要な技法が説明される．個人的な創作活動では境界が曖昧になりがちなプロトタイプ・ラフ・カンプにおいて，それぞれの段階で何をすべきでどのような箇所に気をつけなければいけないかが述べられ，それぞれの役割が明確になる．その他にも多人数で知恵を出しあってモノを作っていく方法にも触れられており，ブレインストーミングにおいては相手の1次的な意見に自分の意見を重ねることが大切だといった実践的なテクニックも入っている．そういった形で，本書をひと通り読めばある程度デザイナーの仕事をなぞることができるレベルの内容が詰まっている．

そして本書の最大の特徴は，見開き1ページ単位で各技法が紹介されているところだ．何かデザインしようと思ったら，ペラペラ本をめくりながら自分がやりたいことを探して，そのページの内容をじっくり読むなり他の文献への足がかりにすれば良い．そういったリファレンスとしての使い方ができるというのは，無限の可能性があるデザインの分野において作りたい作品のアタリを付けるために非常に有効である．

しかし，じゃあそれだけで満足のいくデザインが作れるようになるかといえば，私は必ずしもそれだけでは足りないと考える．そこには，デザインという分野における技法がある程度体系化しながらも，「やってはいけないこと」という人間の感性に基づいた禁則があるからだ．デザインは最終的に人間が見て感じて判断するものだから，いくら技法やテクニックが優れていようとも，人間側が受け付けなければ意味がない．色に関して言えば，例えば同系色の背景と文字が重なって読みにくかったりだとか，原色が強すぎて印象が強烈になるだとかいったことは，初心者が陥りがちな部分だ．そういった「やってはいけないこと」という禁則が，デザインには少なからず存在する．そういった部分において，本書はデザインの成功例は豊富に載っているが，失敗例は非常に少ない．

その点，冒頭で私が薦めた「ノンデザイナーズ・デザインブック」は，「やってはいけないこと」に関する情報が充実している．デザインを体系的に学んでいくなかで，ある手法に関してひと通り説明されたあとには必ず，その手法を使った作品の失敗例と，それを改良した成功例が対になって出てくる．そのため，読者は反例を学びながら手法自体の理解を深めていくことができる．そして，自分がいざデザインするというときには，そういった反例と見比べながら自分のデザインを修正することができる．このような真似や反面教師という関係は，初学者にとってはとっつきやすい部分だと思う．

といったように，これ1冊ですべてを網羅できるといったわけではないものの，本書はデザインにおいて網羅的でまとまりのあるリファレンス的な1冊だといえる．プレゼンテーションなどを作る際には，手元に置いておきたい本だ．

「データ解析のための統計モデリング入門」をひと通り読んだ．本書はGLMからMCMCによる分布推定までの一連の統計モデリングの流れを，生態学における研究の問題に即したテストデータやRを使った解析例とともに解説した本である．本書を書かれた久保さんの講義資料は前々から拝見していたのだが，今回はそれが全体を通して非常によくまとまっている印象を受けた．やはり実例に沿った例題があって，それを解決するためのストーリーが組まれていると，何が問題で何をすべきなのか，そしてその評価方法を含めてハッキリとしていて読みやすい．Rのコードに関しても，コマンドの実行方法からその解釈の仕方まで丁寧に解説が組まれており，数式とのつながりもわかりやすい．個人的には，後半のMCMCの実験に関しては大部分をWinBUGSにお任せで，シミュレーションの過程が少し不明瞭だった感じもするのだが，限られた紙面でMCMCの細かい実装方法を説明するよりかは既存のツールを使ったほうが全体の流れとしても良かったのだろう．

さて，本書を読み始める前に必ず目に通すべきなのが，公式ページにも貼られている各章のつながりや線形モデルの発展の図だ．全体の流れを俯瞰できるほか，実際に中を読みながら今ある問題と図とを照らし合わても，自分が今どういった位置づけの中で何を解いているのかということを把握することができる．

これが非常に参考になったので，個人的にも自分の読書ノートに書きなぐったまとめ図を清書して作ってみた．ほとんど久保さんの図を真似ただけだが，2つの図を組み合わせ感じでストーリーに沿った形で線形モデルの発展をなぞっている．基本的に章の間をつなぐような自問自答の形を取っており，各章には推定するパラメータと手法の代表的なものを付け加えている．