餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

RでべたべたとBHHH法

Berndt-Hall-Hall-Hausmanアルゴリズムは、化石的な最適化アルゴリズムで、実証分析を行っている経済学徒でお世話になったことがある人は、そこそこいると思います。TSPと呼ばれる統計解析パッケージ*1では標準的に使えたのですが、RではmaxLikパッケージ*2などを入れる必要があります。利用困難と言うわけではないですが、人気ではなさそうです。マルコフ連鎖モンテカルロ法が当たり前の時代、古の技法になりつつあると言えるでしょう。しかし、ちょっとした都合でBHHH法の計算手順を確認したいと思います。

1. BHHH法の概要

目的関数として対数尤度関数 \mathcal{l}(y; \beta)が与えられたときのことを考えましょう。yは被説明変数のベクトル、\betaはパラメーターです*3。もっとも基本的なものはニュートン・ラフソン法になります。真の係数を \beta、推定量 b^{(m)} mは反復回数で指数ではない)とすると、

 b^{(m)} = b^{(m-1)} - \alpha^{(m-1)} \frac{\partial \mathcal{l}^{(m-1)}}{\partial{b^{(m-1)}}\partial{b^{(m-1)}}} \frac{\partial{\mathcal{l}^{(m-1)}}}{\partial{b^{(m-1)}}}

になります。ここで、二階微分の計算があると上手く計算されないことが多々あるので、一階微分しかないフィッシャー情報行列の符号を反転させたもので代用しようと言うのが、BHHHアルゴリズムです。原理上、最尤法の計算ぐらいにしか使えません。

2. データセット

サンプルサイズ10のデータセットです。

X <- t(matrix(c(1, 5.75095, 7.13628, 1, 8.12624, 3.43466, 1, 0.470068, 8.32101, 1, 6.97381, 8.56067, 1, 3.01102, 1.86405, 1, 2.73131, 9.86068, 1, 0.639832, 7.30122, 1, 2.02848, 3.03172, 1, 2.27631, 0.179666, 1, 2.59997, 7.71923), 3))

y <- t(matrix(c(-8.67053, 6.59171, -21.7837, -10.9077, 1.78061, -23.3735, -21.2741, -3.17336, 4.87379, -16.6442), 1))

3. 対数尤度関数とグラディエント

正規分布を仮定した線形回帰を最尤法で求めます。対数尤度関数を一階微分したグラディエントは、フィッシャー情報行列の計算のために、説明変数と従属変数を一つづつ取れるようにします。

# 対数尤度関数
objf <- function(p){
  theta <- p[1]
  len <- length(p)
  b <- matrix(p[2:len], len-1, 1)
  s <- sum((y - X %*% b)^2)
  n <- length(y)
  -n*log(theta^2)/2 - s/theta^2/2
}

# グラディエント
objg <- function(p, y=get("y"), X=get("X")){
  n <- length(y)
  theta <- p[1]
  len <- length(p)
  b <- matrix(p[2:len], len - 1, 1)
  res <- y - X %*% b
  dtheta <- -1*n/theta + sum(res^2)/theta^3
  dbeta <- numeric(length(b))
  for(i in 1:length(b)){
    dbeta[i] <- sum(X[, i]*2*res)/theta^2/2
  }
  c(dtheta, dbeta)
}

4. フィッシャー情報行列の計算

フィッシャー情報行列の計算にヘッシアンが要りそうな感じがする*4のですが、ij列が \frac{1}{N} \sum_{n=1}^N \frac{ \partial \mathcal{l}(y_n; \beta) }{ \partial \beta_i } \frac{ \partial \mathcal{l}(y_n; \beta) }{ \partial \beta_j } となる分散共分散行列なので素朴に計算できます。これを観測数で割ったもの*5に-1を乗じて、ヘッシアンの代わりにします。

FIM <- function(p){
  m <- matrix(0, length(y), length(p))
  for(i in 1:length(y)){
    m[i, ] <- objg(p, y[i], matrix(X[i, ], 1))
  }
  (t(m) %*% m)/nrow(m)
}

なお、サンプルサイズを増やすと各セルの値が小さくなりすぎて「システムは数値的に特異です: 条件数の逆数 = ・・・」とエラーが出たりするので、実用的に使うには、もっと改良が要ります。

5. ステップ幅の計算

ステップ幅 \alphaの具体的な計算方法については原論文やアルゴリズム紹介には言及が無いので、「Rで黄金分割探索」で定義した関数を流用します。BFGS法の実装でも、ステップ幅の計算はまちまちのようで、プログラマに任されているのでしょう。

6. BHHH法による推定

ニュートン・ラフソン法とほとんど一緒です。

init <- c(1, 1, 1, 1) # 初期値
param <- init # 更新していくパラメーターの変数
pv <- objf(param) # 初期値の目的関数の値を保存
for(i in 1:300){ # 収束が遅いので繰り返し回数上限は長めに
  g <- matrix(objg(param, y, X), length(param))
  R <- -1*( FIM(param) / nrow(X) )
#
# 以下に置き換えるとニュートン・ラフソン法になる
#  library("numDeriv")
#  R <- jacobian(function(p){ objg(p, y, X) }, param)
#
  # ステップ幅を計算
  a <- gss(function(a){
    # solve(R) %*% g = solve(R, g)
    objf(param - a * solve(R, g))
  }, min=0, max=1, maximize=TRUE) # データセットを変えてみて、うまく収束しなくなった場合はmin=-1にすると・・・

  # パラメーターを更新
  param <- param - a * solve(R, g)

  # 新たなパラメーターで目的関数の値を計算
  v <- objf(param)

  # 改善幅が小さければ終了
  if(abs(pv - v) < 1e-12){
    print(sprintf("pv - v = %f", pv - v))
    break;
  }

  # 計算の途中経過をばらばらと表示
  print(sprintf("i:%d a=%f", i, a))
  print(round(c(param), 6))

  # 現在のパラメーターの目的関数の値を保存
  pv <- v
}

# 推定結果をリストにまとめる
r_bhhh <- list(
  param = param,
  g = matrix(objg(param, y, X), length(param)),
  R = -1*( FIM(param)/ nrow(X) ) ,
  vcov = FIM(param)
)

7. 推定結果の比較

係数はほぼ同じモノになります。

r_nlm <- nlm(function(p){ -1*objf(p) }, init, hessian=TRUE)
r_ols <- lm(y ~ 0 + X)
r <- matrix(c(r_bhhh$param, r_nlm$estimate, NA, coef(r_ols)), ncol(X) + 1)
colnames(r) <- c("BHHH", "nlm/DSM", "lm/OLS")
r

f:id:uncorrelated:20200629151144p:plain
分散共分散行列は・・・もっとサンプルサイズが要るのでしょう。

*1: ライブラリと言う意味ではなく、箱詰めで売られていると言う意味。

*2: 商用ソフトと言う意味ではなく、ライブラリと言う意味。

*3: パラメーターには分散などが含まれるので、\thetaの方が良かったかもです。

*4:教科書に数式の展開で出てくるわけですが、直接計算する練習をさせるものは(たまたまかも知れませんが)見た記憶が無いです。

*5:誤植みたいな部分もあって自信がないのですが、原論文ではそうなっているように読めました。

Rで黄金分割探索

ある化石的な最適化アルゴリズムを書くための準備*1なのですが、黄金分割探索をRで実装してみましょう。計算するだけならばoptimize関数を関数を用いればよいのですが、建築や美術など見た目以外でも黄金比が実用的に使えるのを実感できます。

黄金分割探索は、微分法を使わずに一変数の最適化問題を解くことができるもので、引数の有限探索区間を4分割し、内側の2箇所の区間境界における目的関数の値に応じて、左右どちらかの区間1つを探索区間から除く作業を、収束条件を満たすまで繰り返すアルゴリズムで、探索区間の分割を黄金比にするところに名前の由来があります。数値演算の教科書にはよく載っていますし、Wikipediaあたりには実装に必要な量の説明がされているので、図つきの詳しい説明はそちらを見てください。

#
# 黄金分割探索(目的関数を極大化する引数を0から1までの範囲で探す)
#
gss <- function(vf, min=0, max=1, maxit=30, steptol=1e-6, maximize=FALSE){

# a:=p2-p0, b:=p1-p2, c:=p3-p2と置いて、b/a, c/aが黄金比になるように、p0 < p2 < p4 < p1の初期値を定める。
  gr <- (1 + sqrt(5))/2 # 黄金比
  p0 <- min
  p1 <- max
  l <- (p1 - p0)
  p2 <- p0 + l*1/(1 + gr)
  p3 <- p2 + l*1/(1 + gr)/(1 + gr)

# 初期値における目的関数の値を計算
  v2 <- vf(p2)
  v3 <- vf(p3)

  it <- 0 # 繰り返し回数を初期化

# 終了条件は目的関数の引数が変化しなくなったときと、繰り返し回数が上限に達したとき
  while(steptol < p1 - p0 && it < maxit){

    # 目的関数の大小関係を示す差分を計算
    # 最大化のときは符号を逆にする
    v_diff <- (v2 - v3)*(1 - 2*maximize)

    # 探索区間を縮小し、p1,p2,p3,p4,v2,v3の値を更新する
    if(0 > v_diff){
      p1 <- p3
      p3 <- p2
      v3 <- v2
      l <- (p1 - p0)
      p2 <- p0 + l*1/(1 + gr)
      v2 <- vf(p2)
    } else {
      p0 <- p2
      p2 <- p3
      v2 <- v3
      l <- (p1 - p0)
      p3 <- p2 + l*1/(1 + gr)/(1 + gr)
      v3 <- vf(p3)
    }

    it <- it + 1
  }

# 収束時はp1=p2=p3=p4, 収束しないことは考えない
  p1
}

試しに二次関数の最小化/最大化問題を解いてみましょう。

gss(function(x){
  (x - 0.23)^2
})

gss(function(x){
  -1*(x - 0.61)^2
}, min=0.5, max=1, maxit=30, maximize=TRUE)

そこそこループするので速度的には微妙な気もしますが、0.2300004と0.6100005が計算できます。

*1:準ニュートン法と呼ばれる多変数の最適化アルゴリズム群において、一変数の最適化問題を求めたいところが出てきます。目的関数 \mathcal{l}(y; \beta)が与えられたときのことを考えましょう。もっとも基本的なものはニュートン・ラフソン法になります。真の係数を \beta、推定量 b^{(m)} mは反復回数で指数ではない)とすると、  b^{(m)} = b^{(m-1)} - \alpha^{(m-1)} \frac{\partial \mathcal{l}^{(m-1)}}{\partial{b^{(m-1)}}\partial{b^{(m-1)}}} \frac{\partial{\mathcal{l}^{(m-1)}}}{\partial{b^{(m-1)}}} になります。 \alphaはステップ幅です。このヘッセ行列の部分があるニュートン・ラフソン法では \alpha^{(m-1)}=1でよいのですが、二階部分の計算があると上手く計算されないことが多々あり、一階微分しかないフィッシャー情報行列の符号を反転させたもので代用したり、一回計算したら簡素な更新を続けて使いまわしたりします。代用する場合、移動量が十分に調整されないわけで、ステップ幅を毎回調整しないと上手く計算できないようです。 ステップ幅自体は一定の条件を満たせば何でもよいのですが、説明用にあまり雑なコードを書いているとあれこれ言われてしまいます。ステップ幅の計算のために、黄金分割探索のソースコードを準備します。

信頼区間がマイナスに入らない二項分布の区間推定

感染症の罹患経験など、無作為抽出のyes/noの観測値は二項分布に従うことが知られていますが、その場合も正規分布を用いて区間推定を行うことが多いです*1中心極限定理正規分布に漸近することを利用しているのですが、稀にしか観測されない事象で、全体の生起数が10にも満たない場合は不適切な推定になりがちです。区間の片方の端がマイナスになります。

例えば厚生労働省新型コロナウイルスへの抗体保有率の調査の場合、東京都で観測数が1971名、抗体保有者が2名となり、正規分布から区間推定を行うと0.0010(95%信頼区間-0.00039~0.0024)と言う結果になります。2項分布の生起確率でマイナスはあり得ないので、これはちょっとまずいです。プラスの区間になるように、計算精度を上げましょう。

ベイズの定理を使って、事前分布と尤度関数の積を1に正規化して事後確率を求める方法で、信頼区間を計算するための確率分布をつくります*2。生起確率を求めるときの2項分布の尤度関数は、観測数と生起数を固定した2項分布の確率密度関数となることに注意すると、計算機にやらせる分にはそんなに難しくない作業になります。

n <- 1971 # 観測数
m <- 2 # 生起数

p <- m/n # 生起確率
v <- n*p*(1-p) # 観測数の分散

# 正規化するための尤度関数の積分
s <- integrate(function(p){
    # pを求めるときの2項分布の尤度関数は、観測数と生起数を固定した2項分布の確率密度関数となる
    dbinom(m, n, p)
  }, 0, 1, rel.tol=1e-15)$value

# 確率密度関数
pdf <- function(p){
  dbinom(m, n, p)/s
}

# 累積分布関数
cdf <- function(ps){
  sapply(ps, function(p){
    integrate(function(a){
      pdf(a)
    }, 0, p)$value
  })
}

# 閾値(2分探索で計算)
threshold <- function(a){
  l <- 0
  h <- 1
  m <- 0.5
  for(i in 1:100){
    m <- (l + h)/2
    v <- cdf(m) - a
    if(v == 0){
      break;
    } else if(v < 0){
      l <- m
    } else {
      h <- m
    }
  }
  m
}

# 95%信頼区間とする
a <- 0.025
sprintf("%.4f(%.0f%%信頼区間%.5f~%.5f)", p, 100*(1-2*a), threshold(a), threshold(1-a))

f:id:uncorrelated:20200617045800p:plain

0.0010(95%信頼区間0.00031~0.00366)と計算できました。なお、正規分布からの区間推定は以下で計算できるので、mが50ぐらいのときはほぼ変化なしな事を確認すると、安心して使いやすいと思います。

se <- sqrt(v)/n # 生起確率の標準誤差
rng <- se*qnorm(a, lower.tail=FALSE)
sprintf("%.4f(%.0f%%信頼区間%.5f~%.5f)", p, 100*(1-2*a), p - rng, p + rng)

*1: 計算機の利用が一般的でなかった頃の教科書で、付属している表から紙と鉛筆で計算できる方法が紹介されたためだと思います。

*2:統計哲学的には信用区間と言った方がよいかもです。

Rと手作業で覚えるGLM

glmコマンドに頼り切って二項分布(と言うかロジスティック回帰)やポアソン回帰をしている日々ですが、速度狂は最初から何でもC++で・・・と言う話がある*1ので一般化線形回帰(GLM: Generalized Linear Model)の手順を少し詳しく見てみましょう。最尤法を使っても良いのですが、それではGLMにならないので。

1. GLMとは

GLMの概要については『統計モデル入門』の第3章と第4章と関連した付録AとBを読めば把握できますが、一応、概要を説明します。

GLMは、指数分布族の確率分布のモデルに対して使うことができる推定方法です。パラメーター \thetaのときに観測値 yが観察される確率 f(y;\theta)が、 \exp(A(y)B(θ) + C(θ) + D(y))の形式で書けるとき、 f(y;\theta)は指数分布族の関数になります。ポアソン分布、正規分布、二項分布などが含まれ、これらは A(y) = yとなるので、特に正準形と呼ばれます。

さて、対数尤度関数 \mathcal{l}(y; \beta)が与えられたときの推定を考えましょう。最適化アルゴリズムは色々とありますが、もっとも基本的なものはニュートン・ラフソン法になります。真のパラメーターを \beta、推定量 b^{(m)} mは反復回数で指数ではない)とすると、

 b^{(m)} = b^{(m-1)} - \frac{\partial \mathcal{l}^{(m-1)}}{\partial{b^{(m-1)}}\partial{b^{(m-1)}}} \frac{\partial{\mathcal{l}^{(m-1)}}}{\partial{b^{(m-1)}}}

となります。ニュートン・ラフソン法は収束が速いのですが、いたるところで2階導関数になるヘッシアンが凸を示す値を返さないと、局所最適解にはまりがちです。準ニュートン法と呼ばれる発展的な最適化アルゴリズムは、この2階導関数の計算を省く工夫していますが、GLMも同様です。

GLMの場合は、2階導関数を、ヘッシアンにマイナスをかけた値になるフィッシャーの情報行列(information matrix) \mathcal{I}で置き換えます(スコア推定)。

 \mathcal{I}^{(m-1)}b^{(m)} = \mathcal{I}^{(m-1)} b^{(m-1)} +  \frac{\partial{\mathcal{l}^{(m-1)}}}{\partial{b^{(m-1)}}}

 \mathcal{I}は、分散共分散行列の符号をマイナスにしたものになります*2。情報行列をどうやって計算すればいいのかが問題になりますが、上述の指数分布族の特性と対数尤度関数の一般的特性を前提に式を整理していくと、 \mathcal{I}のi行j列の要素を、

 \mathcal{I}_{ij} = \sum^N_{k=1} \frac{x_{ki}x_{kj}}{var(y_k)} \bigg( \frac{\partial \mu_k}{\partial \eta_k} \bigg)^2

と書くことができます。 Nはサンプルサイズです。ここで \muは観測値 yの期待値、 \etaは説明変数 Xと係数 \betaの線形結合 X\betaで、 \mu \etaはリンク関数 gで対応が定義されます。

 g(\mu) = X\beta = \eta

具体的な導出は、上述の参考文献にざっくり説明*3があるので、そちらを参照してください。推定は以上と具体的なリンク関数を整理することでできますが、もう少し式を整理すると、一般化線形回帰と言う名称の意味が分かります。

以下の要素の対角行列 Wと、

 w_{ii} = \frac{1}{var(y_k)} \bigg( \frac{\partial \mu_k}{\partial \eta_k} \bigg)^2

以下を要素に持つベクトル zを定義し、

 z_i = \sum_{k=1} x_{ik} b_k^{m-1} + (y_i - \mu_i)\bigg( \frac{\partial \mu_k}{\partial \eta_k} \bigg)

スコア法の反復式を整理すると、

 X^t W X b^m = X^t W z

となり、ウェイト付線形回帰の形式に持っていくことができます。 W z bに依存するので有り難味は微妙なのですが。

2. 実践

さて、本題に入り、べたべたとRのコードを貼っていきましょう。

2.1. ポアソン回帰

ポアソン分布の平均と分散になるλの推定は、ウェイト付線形回帰の形式で実行するとそれっぽいです。

# データセットを作成
set.seed(2022)
n <- 50 # サンプルサイズ
X <- matrix(c(rep(1, n), (1:n) %% 3), n, 2)
beta <- matrix(c(1, 2), 2, 1) # 観察できない真のパラメーター
lambda <- X %*% beta
y <- sapply(1:n, function(i){
  rpois(1, lambda=lambda[i])
})

# 推定してみる
b <- matrix(c(0, 0), 2, 1) # 初期値は(0, 0)から
for(i in 1:30){
  W <- diag(c(1/(X %*% beta)))
  z <- (X %*% beta) + (y - X %*% beta)

  L <- t(X) %*% W %*% X
  R <- t(X) %*% W %*% z

  b <- solve(L, R)
}

# 推定結果
b
# 真のパラメーターβと比較
round(b - beta, 6)

# ビルトイン関数のglmでも推定
r_glm = glm(y ~ 0 + X, family = poisson(link = "identity")) # デフォルトはlink = "log"になっているのだが、ここではlink = "identity"

# glmの結果との結果と比較
round(coef(r_glm) - b, 6)

だいたい同じ結果になります。

2.2. ロジスティック回帰

ロジスティック回帰のためだけにGLMを使っている人も多いとか何とか・・・

# データセットを作成
set.seed(2022)
n <- 100 # サンプルサイズ

# 観察できるX
X <- matrix(c(rep(1, n), runif(n, min=0, max=10)), n, 2)

# 観察できない推定パラメーター
beta <- matrix(c(-2, 1), 2, 1)
log_odds <- X %*% beta
pi <- exp(log_odds)/(1 + exp(log_odds))

# 観察できるy
y <- 0 + (runif(n) < pi)*1

# 推定してみる
b <- matrix(0, 2, 1) # 初期値
for(i in 1:30){
  t_lodds <- X %*% b
  t_pi <- exp(t_lodds)/(1 + exp(t_lodds))

  # グラディエント∂l/∂bを計算
  U <- sapply(1:ncol(X), function(i){
    sum(X[,i]*(y - t_pi))
  })
  
  I <- matrix(0, ncol(X), ncol(X))
  for(i in 1:nrow(I)){
    for(j in 1:ncol(I)){
    # log(π/(1-π)) = Xβ = ηから、せっせと微分して整理すると、∂μ/∂η = π(1-π)になる
    # var(y_k)は、平均と分散が等しいポアソン分布のような場合では無い場合、1と置いて数値演算してよい模様。結局、スコア推定の収束条件は同じになる。
      I[i,j] <- sum(X[,i]*X[,j]*t_pi*(1-t_pi))
    }
  }
  
  R <- I %*% b + U
  b <- solve(I, R)
}

# ビルトイン関数のglmでも推定
r_glm = glm(y ~ 0 + X, family=binomial(link="logit"))

# glmの結果との結果と比較
round(coef(r_glm) - b, 6)

*1: 無い。

*2:最尤法で係数の標準誤差として、ヘッシアンの逆行列の対角成分を見る理由を思い出しましょう。

*3: ルベーグ積分を未学習の人は、微分と有限ではない区間積分の交換をしているのが気になるかもです。累積分布関数なので、ルベーグの収束定理を使ってよいわけですが。

使っている名詞で2つの文の相関係数を取ってみる

Twitterで、ぬ氏が「二つの日本語の単語列(10-20文字ぐらい)があって、なるべく共通の単語が多いとか意味合い的に近いものを自動でマッチしたいときって何かいい方法ありますか?」と言うお題を出していたので、お気楽な方法を考えてみました。

わかち書きフレームワークのMeCabを利用して、2つの文を分解、単語1ダミー、単語2ダミー・・・と言う風に文Aと文Bのパラメーターをそれぞれつくって、相関係数を出して見ます。共通の単語が多ければ高い相関係数になるので、それで何かマッチング・アルゴリズムにかければ目的を達することができるでしょう。

説明すると長いわけですが、MeCabをインストールしてパスを通してコマンドラインで動くようにしてから、RMeCabをインストールして以下のように走らせたら終わりです。簡単ですね。

library(RMeCab)

# 比較する文を用意する
text <- c("牧瀬紅莉栖はオレの嫁", "牧瀬紅莉栖はツンデレ", "紅莉栖が生きている")

getWords <- function(text){
  res <- RMeCabC(text)
  value <- sapply(res, function(l){
    l[[1]]
  })
  type <- sapply(res, function(l){
    names(l[1])
  })
  value[type=="名詞"]
}

makeDummies <- function(words_all, words){
  0 + (words_all %in% words)*1
}


# 文ごとに名詞リストをつくる
words <- list()
for(i in 1:length(text)){
  words[[i]] <- getWords(text[i])
}

# 重複する単語は消して全体リストを作る
words_all <- unique(sort(unlist(words)))

# 文ごとに単語ダミーをつくる
words_dummies <- list()
for(i in 1:length(text)){
  words_dummies[[i]] <- makeDummies(words_all, words[[i]])
}

# 行列にまとめる
m <- matrix(unlist(words_dummies), length(words_all))
colnames(m) <- sprintf("text%d", 1:length(text))
rownames(m) <- words_all

# 相関係数を計算
cor(m)

ユーザー辞書を追加した、手元の環境だとこんな感じになりました。なお、ほとんどフィーリングで書いたので、RMeCabのコマンドを確認したら無駄が省ける可能性は大です。
f:id:uncorrelated:20200528183041p:plain

それでもすぐに、クラスター分析に持っていくことができます。

colnames(m) <- text
r_cluster <- hclust(dist(cor(m)), method="ward.D")
plot(r_cluster, main="MeCabで短文をクラスター分析", sub="", xlab="", ylab="", axes=FALSE)

f:id:uncorrelated:20200529112003p:plain

もちろん、具体的に何をするかは見えないので、役に立つかは謎です。文字列が長い場合は、よくあるテキストマイニング手法を真似しましょう。

同義語を揃える

使い道次第ですが、文書中に同義語(e.g. 社員募集=クルー募集)が多数考えられる場合は、同義語テーブルをつくって、一つの単語に揃えるほうが適切になるでしょう。

# 同義語テーブル(ベクターで用意/多い場合はファイルから読み込み推奨)
sv <- c("紅莉栖", "クリスティーナ", "助手", "セレブセブンティーン", "セレセブ", "蘇りし者", "ザ・ゾンビ")

# 最初の単語に揃える関数を用意(上にあわせ場合、MeCabに辞書登録している単語のみ有効)
l <- list()
for(i in 2:length(sv)){
  l[sv[i]] <- sv[1]
}

synonym <- function(words) sapply(words, function(word){
  ifelse(is.null(l[[word]]), word, l[[word]])
})

synonym(c("紅莉栖", "クリスティーナ", "セレセブ", "岡部"))

こんな感じで置換されます。

f:id:uncorrelated:20200529111723p:plain

MeCabのユーザー辞書

ユーザー辞書作成につかったcsvファイルの中身は以下です。

紅莉栖,1291,1291,5000,名詞,固有名詞,人名,名,*,*,くりす,クリス,クリス
クリスティーナ,1291,1291,5000,名詞,固有名詞,人名,名,*,*,くりすてぃーな,クリスティーナ,クリスティー
ツンデレ,1285,1285,1000,名詞,一般,*,*,*,*,*,ツンデレ,ツンデレ
セレブセブンティーン,1285,1285,5000,名詞,一般,*,*,*,*,*,セレブセブンティーン,セレブセブンティー
セレセブ,1285,1285,5000,名詞,一般,*,*,*,*,*,セレセブ,セレセブ
蘇りし者,1285,1285,5000,名詞,一般,*,*,*,*,*,ヨミガエリシモノ,ヨミガエリシモノ
ザ・ゾンビ,1285,1285,5000,名詞,一般,*,*,*,*,*,ザゾンビ,ザゾンビ

MeCabのドキュメントに詳しく書いてありますが、mecab-dict-index -d "C:\Program Files (x86)\MeCab\dic\ipadic" -u C:\path\to\animation.dic animation.csv して、res <- RMeCabC(text, dic="C:/path/to/animation.dic")と言うように使います。

トーラス世界の感染進行モデル

GitHubを使えって絶対に言うな!!!

えーっと、BCGワクチン接種が新型コロナウイルス感染症(COVID-19)に有効説検証論文に付随した話で、同一世界にワクチン接種ありとなしの集団がいたらどうなんの? — と言う御題があって、(ドラクエのマップのような縦横がつながっている)トーラス世界にマス目を区切って、そこに接種ありと接種なしの人々をランダムに配置して、隣接したマス目にいる人からのみ影響を受けるようなトイモデルで、どうなるか検討してみました。

大事なポイントなので強調しますが、以下のような感じで、乱雑にワクチン有りと無しを詰め込みます。

まぁ、ワクチンに効果があればワクチンを打つ方が感染しないですし、集団免疫獲得につれて効果量が特に落ちると言うわけでもないです。ただし、ワクチン有効率は最初が大きく中盤で落ち着きます。


しかし、詳しく特異なモデルをつくってもSIRモデルに近づくんですよね。やはり古典的な微分方程式モデルは強い。

set.seed(202004)

# 感染率β
beta_0 <- 0.50 # ワクチン非接種者
beta_1 <- 0.25 # ワクチン接種者

# 回復率γ
gamma <- 0.19

# ワクチン接種者の比率
theta <- 0.5

# トーラス状の空間を表すn行n列の行列を複数つくる
n <- 100

# 計算する期間
t <- 250

# ワクチン接種者は1で、他は0とする
shot <- numeric(n^2)
shot[order(runif(n^2))[1:as.integer(n^2*theta)]] <- 1

# ワクチン接種の有り無しの人数を数えておく
obs <- tapply(shot, shot, length)

# ワクチン接種の有無を表す行列
m_shot <- matrix(shot, n)

# 最初に感染しているのは1人
inf <- numeric(n^2)
inf[order(runif(n^2))[1]] <- 1

# 感染者を表す行列
m_inf <- matrix(inf, n)

# 暴露レベルを表す行列を返す
m_exposed <- function(m_inf){
  # 感染者の左側を1とする行列(トーラスなのに注意)
  inf <- c(m_inf)
  m_l <- matrix(c(inf[-(1:n)], inf[1:n]), n)
  
  # 感染者の右側を1とする行列
  inf <- c(m_inf)
  m_r <- matrix(c(inf[((n^2-n+1): n^2)], inf[-((n^2-n+1): n^2)]), n)
  
  # 感染者の上を1とする行列
  inf <- c(t(m_inf))
  m_u <- matrix(c(inf[-(1:n)], inf[1:n]), n, byrow=TRUE)
  
  # 感染者の下を1とする行列
  inf <- c(t(m_inf))
  m_b <- matrix(c(inf[((n^2-n+1): n^2)], inf[-((n^2-n+1): n^2)]), n, byrow=TRUE)

  # 0~4の値が返る
  m_l + m_r + m_u + m_b
}

# 完全な免疫を獲得した回復者を表す行列
m_rcv <- matrix(numeric(n^2), n)

# データフレームを作成する
df1 <- data.frame(
  t = 1:t,
  pct_infected_shot = numeric(t),
  pct_infected_no_shot = numeric(t),
  number_of_infected = numeric(t),
  VE = numeric(t)
)

for(i in 1:t){
  # 感染確率を表す行列
  # = ウイルスに暴露されている×社会的距離×{世代に応じたβ}×非回復者
  p <- (beta_0 - c(m_shot) * (beta_0 - beta_1)) * abs(c(m_rcv) - 1)
  # 周囲にいる感染者の人数に応じて、感染確率が上がる
  m_prob <- matrix(1 - (1 - p)^c(m_exposed(m_inf)), n)

  # 感染を判定し、感染者行列を更新
  m_inf[runif(n^2) < c(m_prob)] <- 1
  
  # 回復判定をし、感染者行列と免疫保持者行列を更新
  flag <- (runif(n^2) < gamma)&(m_inf==1)
  m_inf[flag] <- 0
  m_rcv[flag] <- 1

  # 感染経験者は、現在感染者もしくは回復者
  infected <- c(m_inf)|c(m_rcv)
  # ワクチン接種者と非接種者ごとに感染経験者を集計
  noi <- tapply(infected, c(m_shot), sum)
  # ワクチン接種者と非接種者ごとに感染経験率を出す
  pct <- noi/obs
  df1$number_of_infected[i] <- sum(noi)
  df1$pct_infected_no_shot[i] <- pct[1]
  df1$pct_infected_shot[i] <- pct[2]
  df1$VE[i] <- 1 - pct[2]/pct[1]
}

title <- sprintf("トーラス世界の感染進行モデル(%d行%d列,回復率%.2f,接種率%.2f)", n, n, gamma, theta)

windowsFonts(Meiryo = windowsFont("Meiryo"))
par(family="Meiryo", mar=c(5, 4, 4, 2), mgp=c(2, 1, 0), bg="white")

plot(df1$pct_infected_no_shot, col="red", type="l", lwd=2, axes=FALSE, xlab="", ylab="", lty=2, ylim=c(0, 1), main=title)
lines(df1$pct_infected_shot, col="blue", lwd=2)

y_at <- seq(0, 1, 0.25)
y_labels <- sprintf("%d%%", 100*y_at)
y_space <- abs(par()$usr[1] - par()$usr[2])/25
axis(2, at=y_at, labels=y_labels, pos=- y_space, las=2)

x_at <- seq(1, t, length.out=10)
x_labels <- sprintf("%.0f期", x_at)
x_space <- abs(par()$usr[3] - par()$usr[4])/40
axis(1, at=x_at, labels=x_labels, pos=-x_space, las=2)

legend("bottomright", lty=c(2, 1), lwd=c(2, 2), pch=c(-1, -1), col=c("red", "blue"), legend=c(sprintf("ワクチン接種無し(感染率β=%.2f)", beta_0), sprintf("ワクチン接種有り(感染率β=%.2f)", beta_1)), bty="n", y.intersp=1.25)

dev.copy(png, sprintf("%sβ0=%.2f,β1=%.2f.png", title, beta_0, beta_1), width=640, height=480, type="cairo", bg="white")
dev.off()


t2 <- as.integer(t/3)
df2 <- df1[1:t2,]
title <- sprintf("トーラス世界の感染進行モデル(%d行%d列,回復率%.2f,接種率%.2f)", n, n, gamma, theta)

windowsFonts(Meiryo = windowsFont("Meiryo"))
par(family="Meiryo", mar=c(5, 4, 4, 2), mgp=c(2, 1, 0), bg="white")

ylim <- c(0, ceiling(max(df2$pct_infected_no_shot)*10)/10)

plot(df2$pct_infected_no_shot, col="red", type="l", lwd=2, axes=FALSE, xlab="", ylab="", lty=2, ylim=ylim, main=title)
lines(df2$pct_infected_shot, col="blue", lwd=2)

y_at <- seq(0, ylim[2], length.out=5)
y_labels <- sprintf("%.0f%%", 100*y_at)
y_space <- abs(par()$usr[1] - par()$usr[2])/25
axis(2, at=y_at, labels=y_labels, pos=- y_space, las=2)

x_at <- seq(1, t2, length.out=10)
x_labels <- sprintf("%.0f期", x_at)
x_space <- abs(par()$usr[3] - par()$usr[4])/40
axis(1, at=x_at, labels=x_labels, pos=-x_space, las=2)

legend("topleft", lty=c(2, 1), lwd=c(2, 2), pch=c(-1, -1), col=c("red", "blue"), legend=c(sprintf("ワクチン接種無し(感染率β=%.2f)", beta_0), sprintf("ワクチン接種有り(感染率β=%.2f)", beta_1)), bty="n", y.intersp=1.25)

dev.copy(png, sprintf("%sβ0=%.2f,β1=%.2f,前半.png", title, beta_0, beta_1), width=640, height=480, type="cairo", bg="white")
dev.off()

Rでガブリエル比較区間を計算する方法

データを男女などで二つ以上の群にわけて、比較したいときは色々とあります。平均や分散の違いをt検定やF検定などで比較するのが教科書的かつ堅実な方法ですが、統計学に必ずしも詳しくない人々に説明する場合は、ビジュアライゼーションも考えていかないといけません。

この目的でよく使われるのは信頼区間によるエラーバーがついたグラフですが、実は往々にして誤解されています*1。2つの群のエラー・バーは重なっているのに、t検定で統計的に有意な差がでるときがあり、群を複数に分割したときの多重比較のための補正も行われません。比較対象が群ではなくて0などある一点であれば、信頼区間の厳密な説明*2とはちょっと違う使い方ではあるものの、t検定と合致した結果が得られるわけですが。

こういうわけで、群と群をグラフで比較するのには、信頼区間は使わない方が望ましいです。ガブリエル比較区間Gabriel (1978))と言う手法があるので、こちらを使うことをお勧めします。多重比較のための補正も行ってくれます。rgabrielと言うパッケージがあるので、使ってみましょう。

# 信頼区間のエラーバーが重なり、ガブリエル比較区間ではそうでないデータセット
df1 <- data.frame(category = c("C","C","C","C","C","C","T-1","T-1","T-1","T-1","T-1","T-1"), value = c(0.117,-0.718,0.7199,-0.2422,-0.3391,-0.2223,1.0094,1.4306,0.099,1.2839,0.403,0.5616))

# 有意水準(片側)
alpha <- 0.025

# ライブラリを呼び出す.インストールしていない場合は install.packages("rgabriel")を
library(rgabriel)

# 信頼区間とガブリエル比較区間を計算
attach(df1)
df2 <- data.frame(
  category=levels(category),
  count = as.numeric(tapply(value, category, length)),
  mean = as.numeric(tapply(value, category, mean)),
  sd = as.numeric(tapply(value, category, sd))
)
df2$se <- df2$sd/sqrt(df2$count-1)
df2$t_upper <- df2$mean + df2$se*qt(1 - alpha, df=df2$count-1)
df2$t_lower <- df2$mean - df2$se*qt(1 - alpha, df=df2$count-1)
df2$g_upper <- df2$mean + rgabriel(df1$value, df1$category, alpha)
df2$g_lower <- df2$mean - rgabriel(df1$value, df1$category, alpha)
detach(df1)

データフレームdf2の中のt_upperとt_lowerが信頼区間、g_upperとg_lowerがガブリエル比較区間になります。信頼区間では群Cのt_upperは群T-1のt_lowerよりも大きいですが、g_upperはg_lowerよりも小さくなっています。グラフを描くと、以下のようにガブリエル比較区間の方が縮まっているのが分かります。

f:id:uncorrelated:20190409102836p:plain
ガブリエル比較区間と信頼区間

t検定をかけて、どちらが正しいかを見てみましょう。

t.test(df1$value[df1$category=="C"], df1$value[df1$category=="T-1"])

両側検定で5%有意で帰無仮説が棄却され、群Cと群T-1に差があることが分かります。なお、多重比較のための補正は手法によって結果が変わるので、検定結果と整合的になるとは限らないことには注意してください。

グラフのソースコード

gabriel.plotがあるので、plotすればグラフになるのですがビジュアル的にイマイチだったので、以下のように描画しました。

windowsFonts(Meiryo = windowsFont("Meiryo UI"))
par(family="Meiryo", mar=c(2, 1, 1, 1), mgp=c(2, 1, 0))

unit <- 0.1
margin <- 0.05 # 信頼区間のエラーバーとガブリエル比較区間のエラーバーの隙間の幅
ylim <- c(unit*floor(min(df2$t_lower, df2$g_lower)/unit), unit*ceiling(max(df2$t_upper, df2$g_upper)/unit))
xlim <- c(c(0.5, length(df2$mean) + 0.5))
plot((1:length(df2$mean)) - margin/2, df2$mean, type="p", ylim=ylim, xlim=xlim, axes=FALSE, xlab="", ylab="", pch=19, col="red")
arrows(1:length(df2$mean) - margin/2, df2$g_lower, 1:length(df2$mean) - margin/2, df2$g_upper, code=0, col="red", lwd=2)
arrows(1:length(df2$mean) + margin, df2$t_lower, 1:length(df2$mean) + margin, df2$t_upper, code=0, lwd=2, lty=2)
points(1:length(df2$mean) + margin, df2$mean, pch=19)

abline(h=df2$t_upper[1], lty=2)
abline(h=df2$g_upper[1], col="red", lty=3)

abline(h=df2$t_lower[1], lty=2)
abline(h=df2$g_lower[1], col="red", lty=3)

axis(1, at=c(xlim[1], 1:length(df2$mean), xlim[2]), labels=c("",as.character(df2$category),""), pos=ylim[1], las=2)

legend("topleft", c(sprintf("%.0f%%信頼区間", (1-2*alpha)*100), "ガブリエル比較区間"), col=c("black", "red"), lty=c(2, 1), lwd=c(2, 2), bg='white', box.col = 'black', bty="n")

メタアナリシス

rgabriel関数のソースコードを見ていると、群を表すファクターごとの標本標準偏差のベクトル(s)と、ファクターごとのサブサンプル・サイズのベクトル(n)と、有意水準(a)があればガブリエル比較区間は計算できるので、複数の分析のデータをつなぐことも簡単にできると思います。

*1:ダメな統計学」第6章

*2: 真の平均値が観測値だとしたときに所定の観測数のデータをとると、観測される平均値が信頼区間に収まる確率が…と言うようなものです。