可視化で理解しない「負の二項分布」

匿名知的集団であるホクソエムの親分が可視化で理解する「負の二項分布」 - ほくそ笑むで言及している負の二項分布ですが、シンクロニシティというか、同じ匿名知的集団ホクソエムなので当たり前と言えば当たり前なのですが、私も負の二項分布、特に負の二項分布のパラメータを推定するための関数であるglm.nb関数の挙動を調べる必要があったのでまとめる。数式変形も大体理解しているが、参考テキストや資料によって表記のギリシア文字のぶれがあるので、その辺どれを使うのがいいのか決めてからまとめたい。

ダミーデータの作成

glm.nb関数に食わせるためのダミーデータをrnegbin関数(負の二項分布に従う乱数を生成)作成する。もちろんシードは71だ。
ここで

  • 負の二項分布の平均はexp(3.7*x+2)、形状パラメータは10.7

と適当に設定している。このあたりのパラメータの呼称が混乱の原因になるように思うが、ここでは

でいう`shape parameter`を用いて書いた形で議論する(そのほうがおいおいGamma–Poisson mixtureとしての負の二項分布としての見通しがよくなる)。

set.seed(71)
size <- 10^4
x <- -(1:size)/size
shape <- 10.7
mu <- exp(3.7*x+2)
y1 <- sapply(mu, function(m)rnegbin(1, mu=m, theta=shape)) 

モデルの推定

負の二項分布に従うモデルのフィッティングを行うためにglm.nb関数で推定を実施する。表記法は(g)lm関数などと同様formula式でOK。

> library("MASS")
> model_glmnb <- glm.nb(y1 ~ x)
> summary(model_glmnb)

Call:
glm.nb(formula = y1 ~ x, init.theta = 9.775593004, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.1997  -0.9142  -0.4375   0.5292   3.6430  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.02512    0.01309   154.7   <2e-16 ***
x            3.76445    0.03668   102.6   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(9.7756) family taken to be 1)

    Null deviance: 24389  on 9999  degrees of freedom
Residual deviance: 10194  on 9998  degrees of freedom
AIC: 29351

Number of Fisher Scoring iterations: 1


              Theta:  9.776 
          Std. Err.:  0.746 

 2 x log-likelihood:  -29345.396 

ここまでの使い方は適当にググるとでてくるのだが、問題はこの結果の解釈で、それについてはあまり言及されていない気がするので詳細に書くと。。。

  • 切片(Intercept)のはおおよそ: 2.02512
  • 特徴量(説明変数)xの係数: 3.76445
  • Theta(shape変数の値): 9.776

と読む。この係数(2.02512、3.76445、9.776)は元のダミーデータを設定した際の値がほぼ再現されており、良い推定を与えていることがわかる。


また、ここで推定されているもの(predict, fitted関数でとれるもの)は、(lm関数などの結果から類推できるが、)上でいう負の二項分布の平均値(mu)の推定値になる。実際、モデルオブジェクトに入っている値は以下のようにexp(係数×データ)で計算されるもので、上で書いたmuの定義exp(3.7*x+2)そのものになっている。以下では同じ答え(推定値)を出す三通りの方法+厳密な値(私が指定)を算出している。推定値と厳密な値も当然近しいものになる。

> head(fitted(model_glmnb, x))
       1        2        3        4        5        6 
7.574199 7.571348 7.568499 7.565650 7.562803 7.559956 
> head(model_glmnb$fitted.values)
       1        2        3        4        5        6 
7.574199 7.571348 7.568499 7.565650 7.562803 7.559956 
> head(as.numeric(exp(cbind(1, x) %*% beta)))
[1] 7.574199 7.571348 7.568499 7.565650 7.562803 7.559956
> head(mu)
[1] 7.386323 7.383590 7.380859 7.378128 7.375399 7.372671

元のデータの復元

ここから先は上述の親分のブログにあるように、元のデータと同じようにサンプリングされたデータを

  • ガンマ分布からサンプリングする
  • ポアソン分布からサンプリングする

という手順で再現する。そして、その頻度分布が近しいものですねと、ただしくサンプリングされていますねということをチェックする。

まずは、実際に設定したパラメータを用いてサンプリングする。ここで、scale=m/shapeと設定しているが、これは数式変形をがんばると見えてくるところだ。

lambda <- sapply(mu, function(m){rgamma(1, shape=shape, scale=m/shape)})
y2 <- sapply(lambda, function(l)rpois(1, l))
plot(density(y1), col="red")
lines(density(y2), col="blue")


となり、ほぼ一致する。

同様の作業を推定したモデルで実行すると

beta  <- model_glmnb$coefficients
theta <- model_glmnb$theta
lambda <- sapply(model_glmnb$fitted.values, function(m){rgamma(1, shape=theta, scale=m/theta)})
y3 <- sapply(lambda, function(l)rpois(1, l))
hist(y1, breaks=0:max(y1, y3), col="red")
hist(y3, breaks=0:max(y1, y3), col="blue", add=TRUE)


となり、こちらもほぼ一致する。

層別サンプリングした結果を綺麗に持ちたい

層別サンプリングをするには、samplingパッケージのstrata関数使えばいいんだけど、出力とインターフェイスがイケてないのでラップする。


下記の`stratified_sampling`関数を使えば、各層(クラス)の割合を一定に保ったままサンプリングしてくれる。

library("sampling")
stratified_sampling <- function(df, name, size)
{
  sizes <- df %>% 
    count_(name) %>% 
    mutate(size_sampling=round(n/nrow(df)*size)) %>% 
    select(size_sampling) %>% unlist %>% as.numeric
  getdata(df, strata(df, name, size=sizes, method="srswor"))[, colnames(df)]
}

実行結果

> stratified_sampling(iris, "Species", 9)
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
5            5.0         3.6          1.4         0.2     setosa
14           4.3         3.0          1.1         0.1     setosa
25           4.8         3.4          1.9         0.2     setosa
65           5.6         2.9          3.6         1.3 versicolor
82           5.5         2.4          3.7         1.0 versicolor
93           5.8         2.6          4.0         1.2 versicolor
102          5.8         2.7          5.1         1.9  virginica
109          6.7         2.5          5.8         1.8  virginica
147          6.3         2.5          5.0         1.9  virginica

これでCross Validationするのに使うには、どのデータがどこから来たのか覚えるようにしないとだめだなぁ。。。

(g)lmのNA省略(デフォルトの挙動)でうっかりうっかり

データにうっかり欠損(NA)が入ってると、デフォルトでその行がモデル構築&予測結果からごっそり削られるうっかりミスを防止するための備忘録うっかりうっかり。

> library(dplyr)
> #デフォルト(自分で設定したオプション値だけど)のNAなし時の挙動
> getOption("na.action")
[1] "na.omit"

深い意味はない変形+欠損(NA)をirisに加える

> # 適当に追加のラベル(これを予測する)
> df <- iris %>% mutate(label=c(rep(1, 75), rep(0, 75)))
> df <- rbind(df, tail(df, 5) %>% mutate(Species=NA))
> #適当に欠損入れる
> df[1:4, 2] <- NA
> #データをチェック
> head(df)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species label
1          5.1          NA          1.4         0.2  setosa     1
2          4.9          NA          1.4         0.2  setosa     1
3          4.7          NA          1.3         0.2  setosa     1
4          4.6          NA          1.5         0.2  setosa     1
5          5.0         3.6          1.4         0.2  setosa     1
6          5.4         3.9          1.7         0.4  setosa     1
> tail(df)
    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species label
150          5.9         3.0          5.1         1.8 virginica     0
151          6.7         3.0          5.2         2.3      <NA>     0
152          6.3         2.5          5.0         1.9      <NA>     0
153          6.5         3.0          5.2         2.0      <NA>     0
154          6.2         3.4          5.4         2.3      <NA>     0
155          5.9         3.0          5.1         1.8      <NA>     0
> #総行数(155)
> nrow(df)
[1] 155

このまま一般化回帰(回帰でも同じ)しちゃうと・・・

> model <- glm(factor(label) ~ Sepal.Length + Sepal.Width + Species, data = df, family=binomial)
> summary(model)

Call:
glm(formula = factor(label) ~ Sepal.Length + Sepal.Width + Species, 
    family = binomial, data = df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.46047  -0.00005  -0.00003   0.00005   1.43731  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)
(Intercept)         18.5831  2608.2132   0.007    0.994
Sepal.Length         0.7637     0.6688   1.142    0.254
Sepal.Width         -0.5336     1.0893  -0.490    0.624
Speciesversicolor  -21.6373  2608.2112  -0.008    0.993
Speciesvirginica   -42.6643  3595.1381  -0.012    0.991

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 202.289  on 145  degrees of freedom
Residual deviance:  67.956  on 141  degrees of freedom
  (9 observations deleted due to missingness)
AIC: 77.956

Number of Fisher Scoring iterations: 19
> #総行数が155ではなくなっている(欠損分削除されている)
> length(predict(model, type = "response"))
[1] 146
> length(predict(model, newdata=df, type = "response"))
[1] 155

データ削られてるので注意!(一番最後の実行結果はデータ数はいいが、中にNAが詰まってる)

可視化で理解する中心極限定理(ベルヌーイ分布編)

こういう話がある。

非常に素晴らしいので、指数分布じゃなくて、ベルヌーイ分布版をアニメーションにしてみた。
下図は

  • ベルヌーイ分布(コイン投げでいうところの表が出る確率p=0.2)からサンプリングしたデータのヒストグラム(灰色)
  • 平均値の推定値(理論値、青い密度描画、中心極限定理から平均0.2, 標準偏差sqrt(0.2*0.8/データ数)の正規分布に従う)
  • データから推定した平均値、およびその1σ信頼区間(赤色縦棒、1σ信頼区間正規分布の1σにそのまま対応)

となっている。

データが溜まってくると平均値の推定値の精度があがる(赤線の間隔 or 正規分布標準偏差が縮まる)ことが可視化されてわかりやすい。

コードは以下で、ほぼ@hoxo_m氏の上記の記事のパクリ。ggplot2の勉強になりました。

library(animation)
library(ggplot2)
saveGIF( {
  prob <- 0.2
  mean_p <- prob
  sd_p <- sqrt(prob*(1-prob))
  x <- c()
  for(i in 1:50) {
    x <- c(x, sample(c(1,0), 10, prob=c(prob, 1-prob), replace=TRUE))
    mean_x <- mean(x)
    sd_x <- sd(x)/sqrt(length(x))
    p <- ggplot(data.frame(x=x), aes(x=x)) + 
      geom_histogram(aes(y=..density..), binwidth=0.1, alpha=0.5) + 
      stat_function(fun=dnorm, geom="density", color="blue", fill="blue", alpha=0.5, arg=list(mean=mean_p, sd=sd_p/sqrt(length(x)))) + 
      geom_vline(xintercept=mean_x,      color="red", size=2) +
      geom_vline(xintercept=mean_x+sd_x, color="red", size=1) +
      geom_vline(xintercept=mean_x-sd_x, color="red", size=1) +
      scale_y_continuous(labels = function(x) sprintf("%0.2f", x)) + 
      xlim(0, 1.1) + xlab("ベルヌーイ分布の平均値の分布")      
    print(p)
  }
}, cmd.fun = system, interval = 0.4, ani.width=600, ani.height=400)

推定した平均値のバラつき(正規分布の1σ or データから計算した1σ信頼区間)も、理論・データからの推計値がほぼ一致。

> sd_x
[1] 0.01823362
> sd_p/sqrt(length(x))
[1] 0.01788854

ガウシアンを一様乱数でモンテカルロ積分した際の一様乱数の幅依存性

ガウシアンを、一様乱数でモンテカルロ積分した際の結果の、一様乱数の範囲依存性がみたい。

こんな感じか。

f <- function(x){exp(-0.5*x^2)}
N <- 10^3
L <- c(1:10, 20, 50, 100)
y <- numeric(length(L))
for(i in 1:length(L)){
  y[i] <- L[i]/N*sum(f(runif(N, min=-0.5*L[i], max=0.5*L[i])))
}
plot(L, y, type="b", col=2)
abline(h=sqrt(2*pi))

横軸が一様乱数の幅で、縦軸が評価値。黒線が答えで赤線がモンテカルロ積分で評価した値。3σもあれば十分かなと思っていたが、6σくらいないと結果が収束してないように見える。また、幅がでかすぎてもいらないところを評価しまくることになるので、推定精度は下がる。

ありがたい事に添削してもらった。RMSE(=sqrt(M)×エラーバー)も算出してもらっているので、こっちのほうがいいね&可視化ggplot2苦手勢にはありがたい。

時系列データに対するブートストラップ法(ブロック・ブートストラップ法)について

あたまだし

検定やクロスバリデーション等への応用を企図した、サンプル数を水増しするための手法としてブートストラップ法がある。これをRで実行するにはsample関数を使って自分でリサンプリングするコードを実装するか、あるいはbootパッケージのboot関数を用いればいい。

ただ、通常このようなリサンプリングにおいては、例えば、データのレコードの行番号を一様にリサンプリングするなど、"データの順序"を考慮したものとはなっておらず、これはデータの順序に意味があるデータ、特に(時)系列データに対して問題となってくるので、通常のブートストラップ法を適用することはできない。

時系列データに対するブートストラップ法に関しては、まず、大枠としてのブロック・リサンプリング法があり、その構成要素としてブロック・ジャックナイフ法、ブロック・ブートストラップ法が研究されてきた。ブロック・ジャックナイフ法はさほどメジャーではないということなので、ここではブロック・ブートストラップ法のみを取り上げる。また、時系列データには定常性を仮定する。

詳しくは

を見るとよい。論文へのリファレンスも多くオススメ。本記事の内容は当該書籍のRのコード部分の要約に近い。

ブロック・ブートストラップ法

ブロック・ブートストラップ法のアイディアは非常に単純で、

  • データの順序に意味があるのならば、データを1つだけではなく、塊(ブロック)としてリサンプリングすればよかろう

というものである。そうすることで、大局的なデータの相関構造は壊れるが、ブロック内での局所的なデータの相関構造は担保されるだろうというものである。

ブロック・ブートストラップ法の中でも、ここでは

  • 移動ブロック・ブートストラップ法(Moving Block Bootstrap, 以下MBB法)
  • 円形ブロック・ブートストラップ法(Circular Block Bootstrap, 以下CBB法)
  • 定常ブートストラップ法(Stationary Bootstrap, 以下SB法)

の3つを取り上げる。上述の書籍ではこれに加え「非重複ブロック・ブートストラップ法」も取り喘げているが、私的にあまりイケてない感があったので、こいつは端折る。

以下ではbootパッケージのtsboot関数を使うので、これはインストールして読み込んでおく。

install.packages("boot")
library("boot")

移動ブロック・ブートストラップ法(Moving Block Bootstrap, MBB法)

MBB法は例えば、

> x <- 1:5
> x
[1] 1 2 3 4 5

というデータがあった時に、ブロックサイズを3とすると

c(1:3)
c(2:4)
c(3:5)

という3つのブロックにデータ分割し、このブロック自身をリサンプリングするものである。なお、ブロックは元のデータと同じ数だけデータがリサンプリングされるよう、リサンプリングされるその個数が決まる。今回の場合は2個ブロックがあればデータが6個できるので、ブロックは2個リサンプリングされ、はみ出した1個のデータは捨てられる。

tsboot関数を使ってMBB法を行うためには、引数として

  • sim="fixed"
  • endcorr=FALSE

と設定する。その他の引数については、第一引数(tseries)にデータを、第二引数(statistic)にブートストラップ法で計算したい統計量を算出したい関数を指定する。

ここで第二引数にはMBB法を使って計算したい統計量を算出する関数を渡すが、統計量ではなく、リサンプルされたサンプルそのものが欲しい場合には、引数をそのまま返す関数を指定すればよい。Rには引数をそのまま返す関数、force関数があるので、これを指定している。force関数の実装は要するに

force <- function(x){x}

ということだ。

実際にMBB法を使って見ると

> # ブロックの長さを3としたMBB法によりテストデータと同じ長さのブートストラップ標本を4回発生させる
> tsboot(x, force, R=4, l=3, sim="fixed", endcorr=FALSE)$t
[,1] [,2] [,3] [,4] [,5]
[1,]    2    3    4    2    3
[2,]    2    3    4    1    2
[3,]    2    3    4    1    2
[4,]    1    2    3    1    2

となる。使用した引数の意味を改めて書いておくと

  • tseries: (時系列)データ
  • stat: ブートストラップ法により計算したい統計量
  • R: リサンプリング回数
  • l: ブロックの長さ
  • sim: リサンプルの発生方法
  • endcorr: ブロックの端点補正有無

となる。
tsboot関数の返り値自体はbootクラスのオブジェクトになっており、リサンプリングの際の条件なども持っているが、ここでは特段いらないので、その返却値に対して$tと付けることで、リサンプリングの結果だけを取得している。

円形ブロック・ブートストラップ法(Circular Block Bootstrap, CBB法)

これは物理でいう”周期境界条件”をつけたものと考えると良い。具体的には、先ほど同様のデータ

> x
[1] 1 2 3 4 5

に対して、ブロックサイズを3とすると

c(1:3)
c(2:4)
c(3:5)
c(4,5,1)
c(5,1,2)

という5つのブロックにデータ分割し、このブロック自身をリサンプリングするものである。こうすることで先ほどのデータでいうと1や5という数値が2,3に比べて出にくいというMBB法に存在した、"リサンプリングデータに対する偏り"がなくなる。

RでCBB法を実行する上で、MBB法と違う点はendcorr=TRUEとするだけである。また、このendcorr引数はデフォルトでTRUEとなっているので、実際には書かなくていい。

# ブロックの長さを3としたCBB法によりテストデータと同じ長さのブートストラップ標本を4回発生させる
> tsboot(x, force, R=4, l=3, sim="fixed")$t
[,1] [,2] [,3] [,4] [,5]
[1,]    3    4    5    3    4
[2,]    3    4    5    2    3
[3,]    3    4    5    4    5
[4,]    3    4    5    1    2
定常ブートストラップ法(Stationary Bootstrap, SB法)

理論的な背景を調べきれていないのだが、元の時系列が定常であっても、MBB・CBBでリサンプリングした場合には定常性が崩れる場合がある(らしい)。そういう状況を避けるためには、SB法を使うとよい。RでSB法を実行する上で、CBB法と違う点はsim="geom"とするだけである。こうすることで、平均l(ブロックサイズ)になる幾何分布からブロックサイズをサンプリングするようになる。なので、SB法の場合、ブロック数は一定とならない。

# ブロックの長さを3としたSB法によりテストデータと同じ長さのブートストラップ標本を4回発生させる
> tsboot(x, force, R=4, l=3, sim="geom")$t
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    4    5    5    1
[2,]    5    1    1    1    2
[3,]    1    2    3    4    5
[4,]    1    2    3    4    5

ブロックサイズの選び方、

上述の参考書籍に「サブサンプリング法に基づくブロック長さの選択法」なるものが載っていた。疲れたから省略。

まとめ

結局定常ブートストラップ法使っておけばいいのか?なんかタラタラ書いてたら堅い日本語になっちゃった・・・

リサンプリング(復元抽出)で積分値を評価する

の応用かつ、「これで間違ってねぇよな?」という自分への備忘録。インポータンスサンプリングや測度変換とは微妙に分母の形式が違ってくるので間違いやすいんだ。

数式の整理

以下の積分I
 I = \int_1^{L} w(x) f(x) dx
を計算することを考える。積分範囲は以下の関数例にあうように適当に設定してしまった。

これはモンテカルロ法により以下のように近似される。UNIFORMは一様分布のこと。
 I = (L-1)\int_1^{L} \frac{1}{L-1} w(x) f(x) dx \sim \frac{L-1}{N} \sum_{i=1}^{N} w(x_i)f(x_i), \,\,\,\,\, x_i \sim UNIFORM(1, L)

ここで、積分Iの式を、以下のように式変形して考える。
 I = \frac{L-1}{N} sumw \sum_{i=1}^{N} \frac{w(x_i)}{sumw} f(x_i), \,\,\,\,\, sumw=\sum_{i=1}^{N} w(x_i)

ここで定義した
 p(x_i)=\frac{w(x_i)}{sumw}
は各サンプルiに対する重み(確率)とみなすことができるので、この値に比例するようにリサンプリングすると
 I = \frac{L-1}{N} sumw \frac{1}{N}\sum_{i=1}^{N} f(x_i)=\frac{L-1}{N} weight \sum_{i=1}^{N} f(x_i), \,\,\,\,\, x_i \sim p(x_i), \,\,\, weight=\frac{sumw}{N}
として計算することができる*1


できる・・・はずだと信じてやってみる。

やってみる

#適当な関数w, f
w <- function(x){log(x)}
f <- function(x){exp(-x^2)}
#Lを4にしてどんなもんになるか描画


L <- 4
x <- seq(1, L, length.out = 100)
plot(x, w(x)*f(x))
#ふつーのモンテカルロ法(1万乱数)
N <- 10^4
x <- runif(N, min=1, max=L)
(L-1)*sum(w(x)*f(x))/N
#際サンプリングして計算
index <- sample(N, replace=TRUE, prob=w(x))
weight <- sum(w(x))/N
(L-1)*sum(weight*f(x[index]))/N
#実際に積分させてみた答え
stats::integrate(function(x){w(x)*f(x)}, 1, L)

実行結果の箇所だけ示すと

> (L-1)*sum(w(x)*f(x))/N
[1] 0.03627271
> (L-1)*sum(weight*f(x[index]))/N
[1] 0.03558152
> stats::integrate(function(x){w(x)*f(x)}, 1, L)
0.03588273 with absolute error < 4.1e-12

というわけで、ちゃんと一致してそうだ。

補記

ここで定義した
 p(x_i)=\frac{w(x_i)}{sumw}
は各サンプルiに対する重み(確率)とみなすことができるので、この値に比例するようにリサンプリングすると

・・・の流れは
 I = \int_1^{L} w(x) f(x) dx = (\int_1^{L} w(x)dx) \int_1^{L} \frac{w(x)}{\int_1^{L} w(x)dx} f(x) dx
=(\int_1^{L} w(x)dx) \frac{1}{N}\sum_i^N f(x_i), \,\,\,\, x_i \sim \frac{w(x)}{\int_1^{L} w(x)dx}
\sim \left( \frac{L-1}{N} \sum_i^N w(z_i)\right) \frac{1}{N}\sum_i^N f(x_i), \,\,\,\, x_i \sim \frac{w(x)}{\left( \frac{L-1}{N} \sum_i^N w(z_i)\right)}, z_i \sim UNIFORM(1, L)
\sim \frac{L-1}{N} sumw \frac{1}{N}\sum_i^N f(x_i), \,\,\,\, x_i \sim \frac{w(x)}{\frac{L-1}{N} sumw}, z_i \sim UNIFORM(1, L), sumw=\sum_i^N w(z_i)
\sim \frac{L-1}{N} weight \sum_i^N f(x_i), \,\,\,\, x_i \sim \frac{w(x)}{\frac{L-1}{N} sumw}, z_i \sim UNIFORM(1, L), sumw =\sum_i^N w(z_i), weight = \frac{sumw}{N}
とした方が見通しがよさげ

*1:ここで、"二回Nで割る"操作が入っているのを一回だと勘違いしてひどい目にあった。そのためのメモ