2011-12-11

abcのABC (R Advent Calendar 2011 11日目)

近似ベイズ計算 (ABC: Approximate Bayesian Computation) に基づいたパラメータ推定のできるパッケージabcを紹介します．

ABCとは

ABCは，尤度関数の計算が不可能・困難なときに，尤度関数を求めることなく事後分布を推定できる手法です．比較的新しい手法で，集団遺伝学・生態学，疫学などさまざまな分野で応用と開発が進んでいます．
ABCに関しては，kos59125 さんの以下のエントリなどがとても参考になります．

近似ベイズ計算によるベイズ推定 (slideshare)
近似ベイズ計算によるベイズ推定 (捨てられたブログ)

パッケージabc

パッケージabcは，

観察データから計算した要約統計量 (引数 target)
事前分布からサンプリングしたパラメータ (引数 param)
それらをもとにシミュレートした要約統計量 (引数 sumstat)

を与えると，観察データと十分に近いシミュレーションセットを許容して，それらのパラメータを事後分布として返してくれます．要約統計量やパラメータが複数ある場合には，それぞれベクターやマトリックスの型にします (行が各シミュレーション，列が各要約統計量・パラメータ)．
許容範囲の大きさは，引数tolで指定します．
引数methodでは3種のアルゴリズムを選択できます．内訳は以下のとおり．

"rejection" 単純にもっとも近いn%のセットを許容
"loclinear" 局所線形回帰によりパラメータを重み付け
"neuralnet" 非線形回帰によりパラメータを重み付け

引数はほかにもありますが，とりあえず重要なのはこれくらいでしょうか．

作者自身の手による論文 [1] もありますので，ヘルプドキュメントとあわせてご参照ください．

デモ

では，実際にデモをしてみましょう！
5つの値からなる観察データ $\{D_{obs} | 2.2, 3.5, 2.9, 3.6, 2.8\}$ があったとします．簡便のため，正規分布にしたがうことがわかっており，分散 $\sigma ^2$ は既知で~~0.025~~ 0.25 (訂正 111213) とします．近似ベイズ計算のフレームワークで，平均値を算出してみます．真の平均値 $\hat{\theta}$ は3.0としました．

まずこの場合，データを生じるモデルは $f(D|\theta) \sim ND(\theta, 0.25)$ ということになります．
モデルのパラメータ $\theta$ には，事前分布として，平均値0で分散100の正規分布でも与えておきましょう．
データをあらわす要約統計量 $S$ には，5つの値の平均値を用います．
シミュレーションデータと観察データの要約統計量は，それぞれ $S_i$ ， $S_{obs}$ とします．
それらの距離を $\rho(S_i, S_{obs})$ とおきましょう．

試行回数は50万回で，以下のような手順になります．

事前分布から50万個のパラメータ $\theta _i$ をサンプリング
それらの $\theta_i$ をもとに，モデル $f(D|\theta)$ を利用して50万セットのシミュレーションデータ $D_i$ を生成
$D_i$ の要約統計量 $S_i$ を計算
観察データ $D_{obs}$ の要約統計量 $S_{obs}$ に十分近いシミュレーションを採用
$\rho$ が十分小さければ，採用されたシミュレーションのパラメータ群が事後分布に近似できる

abc は，このうち4, 5を計算・評価してくれるパッケージになっています．

まず手順1-3をやってみます．

# Set observed data, iteration number, and others.
data.obs <- c(2.2, 3.5, 2.9, 3.6, 2.8)
iteration <- 500000
sigma.hat <- 0.5
num <- 5
library(abc)

# Set prior distribution for theta.
SetPrior <- function(n){
  rnorm(n, 0, 10)
}

# Simulate data from sampled theta.
SetNormData <- function(theta){
  rnorm(num, theta, sigma.hat)
}

# Model
Model <- function(iteration){
  theta <- SetPrior(iteration)
  data <- sapply(theta, SetNormData)
  ss.sim <- apply(data, 2, mean)
  return(list(theta = theta, ss = ss.sim))
}

# Calculate summary statistics for observed data.
ss.obs <- mean(data.obs)

# Simulate theta and  summary statistics.
sim <- Model(iteration)

この段階で $\theta$ の分布を見ると，当たり前ですが，きちんと事前分布 $ND(0, 10^2)$ にしたがっています．

hist(sim$theta, breaks = 50, xlab = expression(theta),
  main = "All sampled theta")

次に，パッケージabcを利用して手順4-5をやってみます．許容範囲は1%としておきましょう．

# Rejection algorithm (tolerance = 1%).
abc.rej <- abc(target = ss.obs, param = matrix(sim$theta, dimnames = list(NULL, "theta")),
  sumstat = sim$ss, tol = 0.01, method = "rejection")

# Local linear regression algorithm (tolerance = 1%).
abc.lin <- abc(target = ss.obs, matrix(sim$theta, dimnames = list(NULL, "theta")),
  sumstat = sim$ss, tol = 0.01,method = "loclinear", transf = "none")

# Non-linear regression algorithm (tolerance = 1%).
abc.net <- abc(target = ss.obs, matrix(sim$theta, dimnames = list(NULL, "theta")),
  sumstat = sim$ss, 0.01, method = "neuralnet", transf = "none")

abcの適用結果は，print, summary, hist, plotを用いて表すことができます．
hist()でそれぞれ許容されたパラメータの分布を見てみると，以下のような感じになります．
許容範囲が十分小さければ，これらは，事後分布からのサンプルと見ることができます．

range.x <- seq(-10, 10, x <- 0.05)
oldpar <- par(mfrow = c(3, 1))
hist(abc.rej, breaks = range.x)
points(range.x, sapply(range.x, function(y){dnorm(x = y, sd = 10)}) * iteration * 0.01 * x,
  type = "l", lty = "dotted")
hist(abc.lin, breaks = range.x)
points(range.x, sapply(range.x, function(y){dnorm(x = y, sd = 10)}) * iteration * 0.01 * x,
  type = "l", lty = "dotted")
hist(abc.net, breaks = range.x)
points(range.x, sapply(range.x, function(y){dnorm(x = y, sd = 10)}) * iteration * 0.01 * x,
  type = "l", lty = "dotted")

(※ 事前分布の式が間違っていたので訂正しました．誤: dnorm()，正: dnorm(sd = 10)．図は未訂正です．20111211)

上から順に，"rejection"，"loclinear"，"neuralnet"の結果で，点線は事前分布です．

summary()で，以下のように95%信用区間や中央値を見ることができます．

summary(abc.rej)
Call: 
abc(target = ss.obs, param = matrix(sim$theta, dimnames = list(NULL, 
    "theta")), sumstat = sim$ss, tol = 0.01, method = "rejection")
Data:
 abc.out$unadj.values (5000 posterior samples)

              theta
Min.:        2.1308
2.5% Perc.:  2.5357
Median:      2.9951
Mean:        2.9951
Mode:        2.9780
97.5% Perc.: 3.4587
Max.:        3.8066

summary(abc.lin)
Call: 
abc(target = ss.obs, param = matrix(sim$theta, dimnames = list(NULL, 
    "theta")), sumstat = sim$ss, tol = 0.01, method = "loclinear", 
    transf = "none")
Data:
 abc.out$adj.values (5000 posterior samples)
Weights:
 abc.out$weights

                        theta
Min.:                  2.1389
Weighted 2.5 % Perc.:  2.5542
Weighted Median:       2.9940
Weighted Mean:         2.9936
Weighted Mode:         2.9915
Weighted 97.5 % Perc.: 3.4140
Max.:                  3.7821

summary(abc.net)
Call: 
abc(target = ss.obs, param = matrix(sim$theta, dimnames = list(NULL, 
    "theta")), sumstat = sim$ss, tol = 0.01, method = "neuralnet", 
    transf = "none")
Data:
 abc.out$adj.values (5000 posterior samples)
Weights:
 abc.out$weights

                        theta
Min.:                  2.1388
Weighted 2.5 % Perc.:  2.5542
Weighted Median:       2.9941
Weighted Mean:         2.9936
Weighted Mode:         2.9920
Weighted 97.5 % Perc.: 3.4150
Max.:                  3.7763

重みをつけると，結果が多少変わっているのがわかります．
全体的に， $\theta$ の事後分布として，3.0のまわりを最頻値に2.5-3.5の範囲が推定できており， $\hat{\theta} = 3.0$ と整合的な結果です :D

パラメータの重み付け

abcの"loclinear"と"neuralnet"アルゴリズムが何をしているのか，簡単に説明しておきます．

これらは，許容範囲内のパラメータを許容した後， $S_i$ と $S_{obs}$ が完全には一致していないセットに対して，以下のような補正を行ないます．
$\theta_i = m(S_i) + \epsilon _i$
$m$ は回帰式で， $\epsilon _i$ は共通の分散から生成される誤差項です．
さらに，許容されたシミュレーションには， $S_{obs}$ との一致に応じて重みが与えられます．重みは，カーネル $K$ を用いて， $K(\rho(S_i , S_{obs}))$ で与えられます．
”loclinear”は $m$ に線形回帰を，"neuralnet"は非線形回帰を仮定します．"neuralnet"を利用すれば要約統計量セットの次元を減らすことができます．

さらに，引数hcorrにTRUEを指定すると (デフォルト)，以下の式を用いて $S_i$ と $S_{obs}$ の不等分散を補正してくれます．
$\theta_i = m(S_i) + \frac{\sigma(S_{obs})}{\sigma(S_i)} \epsilon _i$

最適な許容範囲の検討

ABCでは，許容範囲の設定が肝になります．大きすぎると事後分布を正確に近似しませんし，小さすぎるとサンプリング効率が低下します．この問題に対し，パッケージabcのcv4abc関数を利用すれば，最適な許容範囲の検討が可能です．

cv4abc関数は，任意の許容範囲について，leave-one-out cross-validationによって，パラメータ推定の正確さとロバストネスを評価します．アルゴリズムは以下のような感じ．

j番目のシミュレーションがランダムに選ばれ，その $S_j$ が $S_{obs}$ のダミーとして使われる
j番目をのぞいたすべてのシミュレーションから，abcでパラメータ推定
これをn回繰り返して，正確に推定できたか評価

評価には，以下の式を用います (summaryで返されるのがこの値)．
$E_{pred} = \frac{\Sigma (\theta_j - \hat{\theta_j})^2}{Var(\hat{\theta_j})}$
$\hat{\theta_j}$ はj番目のシミュレーションの真のパラメータで， $\theta_j$ は推定されたパラメータです．

全シミュレーションの回数分繰り返すのが理想ですが，時間がかかりすぎるので，ランダムにいくつか選んで評価します．回数は，引数nvalで与えます．

rejection samplingで，許容範囲を1%と50%にとって，実際にやってみるとこんな感じになります．

# Rejection sampling (tolerance = 1%).
cv.rej.1 <- cv4abc(param = matrix(sim$theta, dimnames = list(NULL, "theta")),
   sumstat = sim$ss, tols = 0.01, nval = 10, method = "rejection")
summary(cv.rej.1)
Prediction error based on a cross-validation sample of 10

           theta
0.01 0.000571858

# Rejection sampling (tolerance = 50%).
cv.rej.50 <- cv4abc(param = matrix(sim$theta, dimnames = list(NULL, "theta")),
   sumstat = sim$ss, tols = 0.5, nval = 10, method = "rejection")
summary(cv.rej.50)
Prediction error based on a cross-validation sample of 10

        theta
0.5 0.2077764

真の値と推定値の適合をみると以下のようになります (左は許容範囲1%，右は許容範囲50%)．

oldpar <- par(mfrow = c(1, 2))
plot(cv.rej.1)
plot(cv.rej.50)

summaryの結果と図から，許容範囲1%だとよく適合しているが，50%だと推定誤差が大きくなっていることがわかりますね．

考察

ABCの最大の難点のひとつは，計算が非常に非効率なところです．許容範囲を十分に小さくするほど，また要約統計量の次元が増えるほど，より大部分のシミュレーションが棄却されることになり，計算時間が必然的に増大します．
この問題を解決するため，MCMCやSMC (Sequential Monte Carlo) を利用した効率的なサンプリングが考案されていますが [2-5]，パッケージabcではそれらを利用できません．
実際にABCを実装しようとすると，パッケージabcで採用されているrejection samplingベースのアルゴリズムは効率が悪すぎて使いものにならないことが多いように思いますが，ABCのフレームワークで試しにちょっと解析してみたい，というようなときには手軽にできて良いのではないかと思います．

参考

[1] Csillery et al. 2010. abc: an R package for Approximate Bayesian Computation (ABC). arXiv:1106.2793.

[2] Marjoram et al. 2003. Markov chain Monte Carlo without likelihoods. PNAS 100:15324-15328.

[3] Sisson et al. 2007. Sequential Monte Carlo without likelihoods. PNAS 104:1760-1765.

[4] Beaumont et al. 2009. Adaptive approximate Bayesian computation. Biometrika 96:983-990.

[5] Toni et al. 2009. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. J. R. Soc. Interface 6:187-202.

2011-12-09

カテゴリカル変数に対するロジスティック回帰分析

カテゴリカル変数を目的変数にとったロジスティック回帰分析についてまとめてみます．
主に『カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)』藤井良宜 2010年共立出版の第6章を参考にしました．

ロジスティック回帰分析とは

カテゴリカル変数として表された目的変数を，その他の変数で説明するモデルを立てて分析する手法のひとつです．ロジスティック回帰分析では，説明変数 x と目的変数 y の関係が直線的ではなく，シグモイド型に変化していくモデル (たとえば，x が大きくなれば y も大きくなるが，ある程度のところを超えると増加の度合いが鈍っていくというようなかたち) を仮定します．目的変数には，2値のカテゴリー，3つ以上のカテゴリー，順序のあるカテゴリーをとることができます．

目的変数が2値の場合

まず，目的変数が2値で，説明変数が連続変数の場合を説明します．
書籍にならい，パッケージ vcd にあるスペースシャトルのデータ SpaceShuttle を例にとります．スペースシャトルの打ち上げテスト時に，O-リングがひとつでも故障したかどうかを調べたデータです．

library(vcd)
head(SpaceShuttle)
  FlightNumber Temperature Pressure Fail nFailures Damage
1            1          66       50   no         0      0
2            2          70       50  yes         1      4
3            3          69       50   no         0      0
4            4          80       50 <NA>        NA     NA
5            5          68       50   no         0      0
6            6          67       50   no         0      0

ここでは，Temperature と Pressure (どちらも連続的な説明変数) によって，Fail (故障ありかなしという2値の目的変数) を説明するモデルを立ててみます．一般的なかたちは
$log(\frac p{1 - p}) = a + bx$
となります (a, b は定数)．二項分布 (または多項分布) にしたがうデータに対して，連結関数にロジット変換を用いた一般化線形モデル (GLM) ※1 をあてはめて解析する手法と言うこともできます．そうすると，[tex:frac p{1-p} は，ふたつのカテゴリーの確率の比と考えることができます．

解析では，GLMの枠組みでモデルを構築しますので，glm 関数を使用します．

shuttle.tp <- glm(Fail ~ Temperature + Pressure, data = SpaceShuttle, family = binomial)
summary(shuttle.tp)

Call:
glm(formula = Fail ~ Temperature + Pressure, family = binomial, 
    data = SpaceShuttle)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2130  -0.6089  -0.3870   0.3472   2.0928  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept) 14.359775   7.442899   1.929   0.0537 .
Temperature -0.241540   0.109722  -2.201   0.0277 *
Pressure     0.009534   0.008703   1.095   0.2733  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 28.267  on 22  degrees of freedom
Residual deviance: 18.972  on 20  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 24.972

Number of Fisher Scoring iterations: 5

結果を見ると，Pressure は Fail の説明にそれほど効いていないようです．
ですので，Temperature だけで Fail を説明するモデルを立ててみます．

shuttle.t <- glm(Fail ~ Temperature, data = SpaceShuttle, family = binomial)
summary(shuttle.t)

Call:
glm(formula = Fail ~ Temperature, family = binomial, data = SpaceShuttle)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0611  -0.7613  -0.3783   0.4524   2.2175  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  15.0429     7.3786   2.039   0.0415 *
Temperature  -0.2322     0.1082  -2.145   0.0320 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 28.267  on 22  degrees of freedom
Residual deviance: 20.315  on 21  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 24.315

Number of Fisher Scoring iterations: 5

AICを比較してみると，前者 (24.972) に比べて後者 (24.315) のモデルのほうが若干小さい値となっているので，Temperature のみで Fail を説明するモデルのほうがリーズナブルなようです．

モデルで予測される値を出力してみるとこんな感じ．

library(ggplot2)
fail.t.pred <- predict(shuttle.t, type = "response")
fail.obs <- as.integer(SpaceShuttle$Fail[!is.na(SpaceShuttle$Fail)]) - 1
shuttle.df <- data.frame(
	Temperature = rep(SpaceShuttle$Temperature[!is.na(SpaceShuttle$Fail)], 2),
	Fail = c(fail.obs, fail.t.pred),
	Class = c(rep("Observed", length(fail.obs)), rep("Predicted", length(fail.t.pred))))
gg <- ggplot(shuttle.df, aes(Temperature, Fail))
gg + geom_point(aes(colour = Class))

でも，あまりフィットしてるように見えませんね…

目的変数が3つ以上でカテゴリー間に順序がないとき

ここでも書籍にならって，iris のデータを使います．フィッシャーが計測したアヤメのデータですね．

head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

図示するとこうなります．

library(ggplot2)
gg <- ggplot(iris, aes(Petal.Length, Petal.Width))
gg + geom_point(aes(color = Species))

setas, versicolor, virginica の3種を目的変数に，雄しべの計測値を説明変数にとり，アヤメの種を推定するモデルを立ててみます．
今度は，目的変数のカテゴリー数が2より多いので，多項分布モデルを適用して，各カテゴリーの確率を考える必要があります．ひとつのカテゴリーを基準として，それと他のカテゴリーの確率の比を考えます．
以下が，このときのモデルになります．
$log(\frac {P(versicolor)}{P(setosa)}) = a_1 + b_1 Petal.length + c_1 Petal.Width$
$log(\frac {P(virginica)}{P(setosa)}) = a_2 + b_2 Petal.length + c_2 Petal.Width$

解析には，パッケージ nnet の multinom 関数を使用します．これは，多項分布モデル専用の glm 関数のような感じです (使い方は glm とほぼ同じだが，多項分布専用なので，family の指定はない)．

install.packages("nnet")
library(net)
iris.multi <- multinom(Species ~ Petal.Length + Petal.Width, data = iris, maxit = 1e4)
# weights:  12 (6 variable)
initial  value 164.791843 
iter  10 value 12.657828
iter  20 value 10.374056
iter  30 value 10.330881
iter  40 value 10.306926
iter  50 value 10.300057
iter  60 value 10.296452
iter  70 value 10.294046
iter  80 value 10.292029
iter  90 value 10.291154
iter 100 value 10.289505
iter 110 value 10.289258
iter 120 value 10.287417
iter 130 value 10.286160
final  value 10.286136 
converged

summary(iris.multi)
Call:
multinom(formula = Species ~ Petal.Length + Petal.Width, data = iris, 
    maxit = 10000)

Coefficients:
           (Intercept) Petal.Length Petal.Width
versicolor   -24.48672      7.26604    8.626438
virginica    -69.57214     12.99170   19.045286

Std. Errors:
           (Intercept) Petal.Length Petal.Width
versicolor    59.26616     50.07007    110.0120
virginica     60.79077     50.12256    110.0756

Residual Deviance: 20.57227 
AIC: 32.57227

モデルが複雑なので，内部でニューラルネットワークを利用した計算を動かしているそうですが，勉強不足なのでここでは触れません (汗)．
しかし何にせよ，setas に比べて versicolor と virsinica は Length も Width も大きく (Coefficients が正のため)，それらが大きくなるほど virsinica である確率が高くなる (virsinica の Coefficients のほうが大きいため)，という解釈ができるようです ※2．

目的変数が3つ以上でカテゴリー間に順序があるとき

パッケージ MASS に含まれる polr 関数で比例オッズモデルを利用するようですが，理解しきれていないのでまた今度．

注

※1 一般化線形モデルについては，生態学のデータ解析をご参照ください．
※2しかしこの場合は，判別分析を採用したほうが良いような気もしますが…

2011-12-08

ループのスピードアップ

Rのループを (しかたなく) 使うときに．
このあたりのエントリを参考にしました．

Another aspect of speeding up loops in R (me nugget)
How to speed up loops in R (Revolutions)

例1のように，まわすたびにオブジェクトのサイズが変わっていくようなループでは，そのたびごとに新たなメモリを確保する必要があるため，パフォーマンスが極端に低下します．
その一方，例2のように，最初に大きなオブジェクトを指定して，その中に値を代入していくようなループなら，極端な速度低下を防ぐことができます．
代入していくオブジェクトのサイズがわからなければ，例3のようにlistに入れていけば良いとか．

…基本中の基本ですが，パフォーマンスが100-150倍くらい向上したのには驚きました．
例4では，同じことをsapply関数を使ってやっていますが，やはりこちらも早いですね．

nrow <- 100
ncol <- 5000

例1

x.bind <- rnorm(nrow)
system.time(
		for(i in 2:ncol){x.bind <- cbind(x.bind, rnorm(nrow))}
)

   ユーザ   システム       経過  
    28.474     11.258     40.090

例2

x.array <- array(0, dim = c(nrow, ncol))
system.time(
		for(i in 1:ncol){x.array[ , i] <- rnorm(nrow)}
)

   ユーザ   システム       経過  
     0.306      0.003      0.308

例3

x.list <- list(rnorm(nrow))
system.time(
		for(i in 2:ncol){x.list[[i]] <- rnorm(nrow)}
)

   ユーザ   システム       経過  
     0.475      0.017      0.492

例4

x.apply <- rep(nrow, ncol)
system.time(
		x.apply <- sapply(x.apply, rnorm)
)

   ユーザ   システム       経過  
     0.290      0.006      0.294

2011-10-24

行列の各要素ごとに異なる引数を与えて，同じ関数を適用

applyを使って，行列の行や列ごとに，異なる引数を与えて同じ関数を適用したいときの方法です．

まずこんな行列を作成

> (xx <- t(array(1:12, dim = c(4, 3))))
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12
>

列ごとに異なる引数を与えて，各行に同じ関数を適用したいとき．
例えば，第1列なら1，第2列なら2，第3列なら3，第4列なら4を引きたいとき．

> t(apply(xx, 1, function(x){x - 1:4}))
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    4    4    4    4
[3,]    8    8    8    8
>

行ごとに異なる引数を与えて，各列に同じ関数を適用したいとき．
例えば，第1行なら0，第2行なら4，第3行なら8を引きたいとき．

> apply(xx, 2, function(x){x - c(0, 4, 8)})
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
[3,]    1    2    3    4
>

それだけです．

2011-09-19

3次元のデータをグラフにする

3次元のベクトルで表されるデータをグラフ化する方法について，いくつかまとめてみました．

使用するデータ

volcano

Rに組み込まれているデータセットで，Maunga Whau山の標高データだそうです．10m * 10mのグリッドで，87rowと61columnからなります．rowは東から西へ，columnは南から北へのラインと対応しており，最小値は94，最大値は195です．

まずはデータセットの準備．

# あらわしたいデータのラベルを作成
EastWest <- 1:nrow(volcano) * 10  # rowの数は87
SouthNorth <- 1:ncol(volcano) * 10  # columnの数は61
# 頂上の位置を変数に入れておく
mountaintop <- which(volcano == max(volcano), arr.ind = TRUE) * 10

グラフを描く

image
値に応じた色を格子上に示す関数です．

主要な引数は以下のとおり．
x, y: プロットされるzの位置を表す座標．昇順でなければならない．
z: プロットされる値のmatrix．NA値も許される．
col: プロットされる値につける色．

備考
xはnrow(z)+1かnrow(z)と等しくなる必要がある．前者ならセルの境界に，後者ならセルのmidpointに，値がプロットされる．
NA値は透明になる．
zのmatrixはf(x[i], y[i])のtableとして解釈される．よって，xはrowの数に，yはcolumnの数に，対応している必要がある．つまり，データセットとグラフは90度回転した関係にあるということですね．

実例

image(EastWest, SouthNorth, volcano, col = heat.colors(100))
# col の値には他にこんなのもあります．
# rainbow(7), heat.colors(10), topo.colors(100), terrain.colors(100), gray((0:100)/100), gray((100:0)/100)
# つづけて頂上の位置にプロットを置きます
points(mountaintop[1], mountaintop[2], pch =  20, col = "blue")
# 等高線を描くcontour関数も合成してみましょう．今回は標高175mの位置に．
contour(EastWest, SouthNorth, volcano, levels = 175, add = TRUE)
# levels = seq(90, 200, by = 5) などもできます．

以下のようなグラフが描けます．

persp

x-y平面に遠近図を描く関数です．

主要な引数
x, y: プロットされるzの位置を表す座標．昇順でなければならない．
z: プロットされる値のmatrix．NA値も許される．
col: プロットされる値につける色．
theta, phi: 視座の角度．前者は横向きの回転，後者は縦向きの回転．
r: 視野とプロットの距離．
d: 遠近感の強さ．1以上なら遠近感が弱まり，1以下なら強調される．
scale: FALSEなら，全軸共通のスケーリングがなされる．
border: 表面に描かれる線の色．NAなら線を描かない．
ltheta, lphi: 指定した値の視座から光が投影されたような陰がつく．
shade: 陰がつく．1に近いほど濃い陰になる．
ticktype: "detailed"なら軸にスケールがプロットされる．デフォルトでは"simple"．
nticks: スケールの数を指定．

実例
ここでは特に，2次元のグラフに使われる関数を3次元のグラフに描くtrans3d関数を使って，点や線も追加してみます．

mtMaugnaWau <- persp(EastWest, SouthNorth, volcano, theta = 25, phi = 30, scale = FALSE,
  col = "green", border = NA, ltheta = 120, shade = 0.7, ticktype = "detailed", cex.axis = 0.8)
# 頂上の追加
points(trans3d(mountaintop[1], mountaintop[2], max(volcano), pmat = mtMaugnaWau), col = "red", pch = 16)
# 頂上から標高150mラインにある一点に直線を引く
level150 <- which(volcano == 150, arr.ind = TRUE) * 10
x <- 6  # 150mラインのうちのどの点に直線を引くか．適当に決めて6にしました．
lines(trans3d(c(mountaintop[1], level150[x, 1]), c(mountaintop[2], level150[x, 2]),
  c(max(volcano), 150), pmat = mtMaugnaWau), col = "blue")

以下のようなグラフが描けます．

scatterplot3d

値を3次元にプロットする関数です．CRANからパッケージのダウンロードが必要です．
scatterplot3d: 3D Scatter Plot

主要な引数
x, y, x: それぞれの座標．
scale.y: y軸のスケール (奥行き) の調節．
angle: x-y軸の角度
highlight.3d: TRUEなら，y軸の値に応じて色の濃さが変わる．

実例

# データセットの準備．
EastWest3d <- rep(EastWest, length(SouthNorth))
SouthNorth3d <- vector(length = 87 * 61)
for(i in 1:61){
  SouthNorth3d[((i - 1) * 87 + 1):(i * 87)] <- SouthNorth[i]
}
# プロットの実行．
library(scatterplot3d)  # パッケージの呼び出し．
mtMaugnaWau3d <- scatterplot3d(EastWest3d, SouthNorth3d, volcano, scale.y = 0.5, highlight.3d = TRUE)
# plane3dメソッドやpoints3dメソッドで平面や点の追加が可能です．
# 標高150mラインに平面の追加．
mtMaugnaWau3d$plane3d(150, x.coef = 0, y.coef = 0, lty.box = "dashed") 
# 頂上の追加．
mtMaugnaWau3d$points3d(mountaintop[1], mountaintop[2], max(volcano), col = "blue", pch = 19)

以下のようなグラフが描けます．

どんなときに使うか

まず，連続的な値をべたーっと面で描きたいときには，imageやperspが適しているでしょう．imageは直感的にイメージがつかみにくいですが，全体が見渡せるので，ひとつひとつの値をある程度きちんと検証できます．一方，perspは直感的にイメージがつかめますが，山や谷で見えなくなる部分があり，データをすべて見せきれないかもしれません．
次に，線・点など，面に比べるととびとびになる値を3次元にプロットしたいときは，scatterplot3d関数が適しているでしょう．引数も充実しており，ラベルやスケールなどもきちんと描くことができます．

※ 他にもこんな関数がある，などご存じでしたら，ぜひトラックバックなどでお教えいただければ幸いです！

2011-09-15

確率密度関数の生成

正規分布を例に，Rで確率密度関数を生成してみます．

まずは準備

> N <- 100				# 各系列の長さ
> mean <- 0; sd <- 0.5		# 分布のパラメータ
> x <- ppoints(N)			# 0から1までの100分位点

確率密度関数の作成とグラフ化

> y <- qnorm(x, mean, sd)	# 分布を x : (1-x) に分割する値を返す
> z <- dnorm(y, mean, sd)	# yにおける密度関数の値を返す
> plot(y, z)

プロットの結果は以下のようになります．

確率の計算

> calcZ <- function(x){ dnorm(x, mean, sd)}		# 密度関数の定義
> integrate(calcZ, - sd, sd)					# ±1SD の範囲に入る確率
0.6826895 with absolute error < 7.6e-15
> integrate(calcZ, -2 * sd, 2 * sd)				# ±2SD の範囲に入る確率
0.9544997 with absolute error < 1.8e-11

こんな感じに，かの有名な0.683とか0.955とかが返されてますね．

- - - -

参考
Rにおける確率分布 (RjpWiki) http://www.okada.jp.org/RWiki/?R%A4%CB%A4%AA%A4%B1%A4%EB%B3%CE%CE%A8%CA%AC%C9%DB
Rで学ぶデータサイエンスベイズ統計データ解析, 姜興起, 2010, 共立出版社

2011-09-06

ベイズ統計入門

statistics

涌井良幸. 2009, 『道具としてのベイズ統計』, 日本実業出版社, 東京.
を読んでの覚え書きです．私の勝手な解釈・誤植・間違いなどあるかと思いますので，参照される際は原典にあたることを強くおすすめします．
ちなみに，私の読んだベイズ統計の入門書のなかでは，この本が一番，わかりやすさと専門性のバランスが程良いと感じました．実例で用いているのがExcelなのが，唯一の難点かもしれませんが…

ベイズ理論とは

ダイナミックな統計手法
・事前情報を取り込み，確率 (確率分布) を更新していく柔軟さ
・階層モデルで多数の分布を融合
・「真の母数」自体を確率変数として考える

例えばコイントスをしたとして，表裏の出る確率が不明であっても，事前知識によって値を設定して計算できる．また，計算をくり返すなかで事前知識を更新していき，よりもっともらしい確率モデルに近づけていける．
言い換えると，母集団における「真の値」を仮定せず，得られたデータを素直に理論に取り込んでいくということ．

ベイジアンネットワーク

因果関係をノードとリンクで表現した確率モデル．それぞれのノードで事象の確率を与えておくことで，ある変数の値が決まったとき，他の未観測の変数の確率分布を求めることができる．

例えば，
$T$ : 友達からメールが来た
$M$ : 迷惑メールが来た
$F$ : 迷惑フォルダにメールが振り分けられた

について，以下のような確率を仮定する．

$P_t(T= 0) = 0.95$ 友だちからメールが来ない確率は0.95
$P_t(T = 1)= 0.05$ 友だちからメールが来る確率は0.05

$P_m(M = 0) = 0.90$ 迷惑メールが来ない確率は0.90
$P_m(M = 1) = 0.10$ 迷惑メールが来る確率は0.10

$P_f(F = 0 | T = 0 \cap M = 0) = 0.99$ メールが来ていない
$P_f(F = 1 | T = 0 \cap M = 0) = 0.01$ 迷惑メールが来たことになっている！
$P_f(F = 0 | T = 1 \cap M = 0) = 0.95$ 友達からのメールを正しく受信
$P_f(F = 1 | T = 1 \cap M = 0) = 0.10$ 友達からのメールを誤って迷惑フォルダに
$P_f(F = 0 | T = 0 \cap M = 1) = 0.09$ 迷惑メールを誤って受信フォルダに
$P_f(F = 1 | T = 0 \cap M = 1) = 0.91$ 迷惑メールを正しく排除
$P_f(F = 0 | T = 1 \cap M = 1) = 0.45$ 友達が送ってきた迷惑メールをそのまま受信
$P_f(F = 1 | T = 1 \cap M = 1) = 0.55$ 友達が送ってきた迷惑メールを正しく排除

これを用いることで，例えば $P(F|T)$ などを計算できる．
ネットワークを重層化していくことで，連鎖する現象の一部について推定が可能になる．
ノードの部分においた仮定を，実験・観測・経験などで更新していくことも可能．

条件付き確率とベイズの定理

まず

$A$ と $B$ の同時確率 $P(A \cap B)$
$A$ のもとで $B$ が起こる条件付き確率 $P(A|B)$
$P(B|A) = P(B \cap A) / P(A)$
$P(A \cap B) = P(B|A)P(A)$ …乗法定理

乗法定理より，
$P(A \cap B) = P(B|A)P(A) = P(A|B)P(B)$
$P(A \cap B)$ を消去して
$P(A|B) = P(B|A)P(A) / P(B)$ …ベイズの定理

$A$ , $B$ をそれぞれ $\theta$ (仮説を代表する母数) と $D$ (data) と置き，離散的な点確率ではなく連続的な確率分布を仮定して，
$\pi(A|B) = f(B|A) \pi( \theta) / P(D)$
これは， $D$ が与えられたときに $\theta$ が成り立つ確率 (事後確率分布) を， $\theta$ のもとで $D$ が成り立つ確率を与える関数 (尤度関数) と $\theta$ についての事前知識 (事前確率分布) から逆算している，ということになる．

さらに

ベイズ統計では $D$ が与えられた後のことを考えるので， $P(D)$ は定数となる．
よって，
$\pi(A|B) \propto f(B|A) \pi( \theta)$
「事後分布は，尤度と事前分布の積に比例する」ということ．

ところが，尤度と事前分布の積の計算は複雑になることが多い．そこで，以下2つの方法が取られる．
a. 自然な共役分布を利用 (尤度関数と相性の良い事前分布を用いて，公式的に事後分布を算出)
b. MCMC (複雑なモデルをそのまま用いて，多数回サンプリングによって近似的に分布を抽出)

自然な共役分布

尤度関数と相性のよい事前分布を用いて，計算を簡単にする．
一覧

事前・事後分布	尤度	表記
ベータ分布 [1, 2]	二項分布[3, 4]	$Be(p, q)$
正規分布 [5]	正規分布	$D( \mu, \sigma^2)$
逆ガンマ分布 [6]	正規分布 (分散未知のとき)	$IG( \alpha, \lambda)$
ガンマ分布	ポアソン分布	$Ga( \alpha, \lambda)$

ベータ分布: 0-1の値をとる一様分布に従う確率変数 $(p + q - 1)$ 個を大きさの順に並べかえたとき，小さい方から $p$ 番目の確率変数 $x$ の分布のこと．ただしベイズ統計では，この確率的性質というよりは，関数の形 (計算のしやすさ) が利用されている．
正規分布に従う尤度関数: 分散が未知のときは， $\mu$ は正規分布に， $\sigma^2$ は逆ガンマ分布に従うと仮定して計算する．
ガンマ分布: 統計的性質はあまり調べられていない．ベイズ統計で重宝されるのは，関数の形ゆえ．

[1] ベータ分布 (Wikipedia) http://ja.wikipedia.org/wiki/%E3%83%99%E3%83%BC%E3%82%BF%E5%88%86%E5%B8%83
[2] ベータ分布 (NtRand) http://www.ntrand.com/jp/beta-distribution/
[3] 二項分布 (Wikipedia) http://ja.wikipedia.org/wiki/%E4%BA%8C%E9%A0%85%E5%88%86%E5%B8%83
[4] 二項分布 (青木繁伸) http://aoki2.si.gunma-u.ac.jp/lecture/Bunpu/nikou.html
[5] 正規分布 (Wikipedia) http://ja.wikipedia.org/wiki/%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83
[6] ガンマ分布 (Wikipedia) http://ja.wikipedia.org/wiki/%E3%82%AC%E3%83%B3%E3%83%9E%E5%88%86%E5%B8%83

MCMC

複雑なモデルをそのまま用いて，多数回サンプリングによって近似的に分布を抽出する．
マルコフ連鎖: ひとつ手前の状態の性質のみから，次の状態が決まる確率過程．適度なランダムさを有する (完全にランダムだと，ひとつ手前の状態すら関与しないことになる)．
モンテカルロ法: 与えられた関数を再現するように点列を抽出して，近似的に関数の積分値を求める方法． $f(x)$ の関数値に比例した密度の有限個の点をサンプリング．

MCMCを実行する方法には以下のようなものがある．

ギブス法

各変数に対し順番に，それ以外の変数を固定したときの条件付き確率分布から点をサンプリングする試行をくり返す方法．

2変数からなる確率分布を例にとると，
1. 初期値 $y_0$ を指定
2. $P(x, y_0)$ から $x_1$ をサンプリング
3. $P(x_1, y)$ から $y_1$ をサンプリング
4. $P(x, y_1)$ から $x_2$ をサンプリング
5. 以下繰り返し
となる．

事後分布が複雑なものはサンプリングが困難なため，自然な共役分布と相性が良い．

メトロポリス法

確率が高い方へはすみやかに移動し，低い方へはなかなか移動しないような動きでサンプリングする方法．
1. 現在の値 $x_t$ からステップ間隔 \epsilonを決め， $x'( = x_t + \epsilon)$ を調べる
2. 以下を実行
$P(x') \geq P(x_t) \to x_{t+1} = x'$
$P(x') < P(x_t) \to$ 確率 $r$ で $x_{t+1} = x$ '，確率 $1-r$ で $x_{t+1} = x_t$
3. $r = P(x') / P(x_t)$

多変数の場合は， $x$ や $\epsilon$ がベクトルになる．
ステップ間隔の設定が重要で，小さすぎると局所最大値にとらわれ，大きすぎると計算がたくさん必要になる．

メトロポリス・へイスティング法では， $r$ の計算式に提案密度が付加されより一般的になる．

MCMCに必要な性質

1. 詳細釣り合い条件
任意の $t$ に対し，
$P(w_t)S(w_t \to w_{t+1}) = P(w_{t+1})S(w_{t+1} \to w_t)$
が成り立つこと．
これより，必然的に，
$P(w_t) > P(w_{t+1})$ なら
$S(w_t \to w_{t+1}) < S(w_{t+1} \to w_t)$
となる．逆も成り立つ．
つまり，確率の高いほうに動きやすく，低いほうには動きにくい，という性質．

2. エルゴート性
任意の2つの状態 $w$ と $w'$ のあいだの遷移確率が，有限個の0でない $S$ の積で表される，という性質．
有限個のステップで，任意の2つの状態を行き来できるということ．

階層ベイズ法