餡子付゛録゛

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

モンテカルロ法で円周率を計算するベンチマーク

モンテカルロ法で円周率を計算するベンチマークで、Juliaと比較してgccが悪い言われているのですが、本当なのか検証してみましょう。

f:id:uncorrelated:20200927174314p:plain

CとJuliaでこの分岐が単純なマイクロベンチマークを実行すると、乱数生成の速度が結果を大きく左右しえます。後方互換性の問題で、Cの標準の乱数生成関数は線形合同法になっています。Juliaの方は標準がdouble precision SIMD oriented Fast Mersenne Twister pseudorandom number generator (dSFMT)で、明示的にMersenneTwisterを使った場合もdSFMTを用いているようです。Cの方もdSFMTを用いてベンチマークを取ることにします。

Cのソースコード

dSFMTがよくチューニングされているので、私には工夫の余地は無く、オーソドックスに書いています。なお、生成した実行ファイルをreadelfで調べたところ、dSFMTがヘッダーファイルでinline static宣言している関数はオンライン展開されていました。なお、10回計測して、一番短い経過時間を採用しています。

#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include "dSFMT.h"
#include <float.h>
#include <math.h>

/*
dSFMT-src-2.2.3.tgz を展開したディレクトリで、ソースコードのファイル名をpi_mc_dsfmt.cとして、

clang -o pi_mc pi_mc_dsfmt.c dSFMT-src-2.2.3/dSFMT.c -IdSFMT-src-2.2.3/ -DDSFMT_MEXP=521 -msse2 -DHAVE_SSE2 -O3 -lm

コンパイルしました。gccとclangはオプションが同じなので、gccとの比較は容易に行えます。
*/

double pi_mc(long n_iter) {
  double  x, y;
  dsfmt_t rng;

  dsfmt_init_gen_rand(&rng, time(NULL));

  int n_in = 0;
  for (long i = 1; i <= n_iter; i++) {

    x = dsfmt_genrand_open_open(&rng);
    y = dsfmt_genrand_open_open(&rng);

    if (x * x + y * y < 1.0)
      n_in++;
  }

  return (double)4.0*n_in/n_iter;
}

int main(){
  struct timespec t1, t2;
  double  dt, sum_dt, min_dt;
  double  r = 0;
  size_t  n;
  int  j, nos = 10;

  printf("log10(n)\telapsed time (seconds)\n");

  for(n=1; n<=1000000*1000; n*=10){
    sum_dt = 0;
    min_dt = DBL_MAX;
    for(j=0; j<nos; j++){
      clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t1); // clock()を使って計測しても、ほぼ同じ数字になります

      r = pi_mc(n);

      clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t2);

//      fprintf(stderr, "estimated pi: %f\n", r);

      dt = (double)(t2.tv_sec - t1.tv_sec) + (double)(t2.tv_nsec - t1.tv_nsec)/1000000000;
      
      sum_dt += dt;
      if(min_dt > dt)
        min_dt = dt;
    }

    printf("%.0f\t%.09f\n", log10(n), min_dt);
  }

  return 0;  
}

Juliaのスクリプト

Cのとだいたい同じですが、グローバル変数へのアクセスが無いように一番外側ループ部分も関数内にしまいました。

using Random
# using Random: default_rng

function pi_mc(n_iter)
n_in = 0
# rng = default_rng()
rng = MersenneTwister(time_ns())
for i in 1:n_iter
x, y = rand(rng), rand(rng)
n_in += ifelse(x^2 + y^2 ≤ 1, 1, 0)
end
4n_in/n_iter
end

function loop(s, e)
  n = s
  while(n <= e)
    sum_dt = 0;
    min_dt = 3600;
    for i in 1:10
      t1 = time_ns()
      r = pi_mc(n)
      t2 = time_ns()
      dt = (t2 - t1)/1000000000
      sum_dt = sum_dt + dt
      if(min_dt > dt)
        min_dt = dt
      end
#      println("estimated pi: ", r)
    end
    println(log10(n), "\t", min_dt)
    n = n*10
  end
end

println("log10(n)", "\t", "elapsed time (seconds)")

loop(1, 10^9)

SIMD拡張命令セットでgcc用のsum関数を書いてみる

とあるJuliaを使った教材に、JuliaやNumPyのsumはSIMD拡張命令セットを使うから速いけれども、C言語でsumを書くとそうではないから遅いような記述があって、CからSIMDを叩けば良いだろうと思ったので、私にしては珍しくC99のコードで、単浮動小数点のsumを書いてみました。gccでしかコンパイルできません。

謳い文句によるとgccは自動的にSIMD拡張命令セットを使うことになっていますが、実際のところはsumのような単純なものでも使いません。こんなにSIMD向きの単純計算でも、128bitsを32bits×4で使っても2倍ちょっとしか速くならないわけで、倍精度だとAVX-256以降でないと意味がないからでしょうか。

なお、C++のテンプレート・ライブラリEigenは行列演算などにSIMD拡張命令セットを使ってくれるので、こういう苦しい記述をしないで浮動小数点演算で速度を出したい人はC++/Eigenを使うべきです。CPUのサポートしている命令セットごとにコードを書き換えないといけないって、手間隙な上に混乱してきますからね。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <emmintrin.h>
#include <time.h>
#include <cpuid.h>

float  x87_sum(float *X, int n){
float s = 0.0;
for (size_t i = 0; i < n; ++i) {
s += X[i];
}
return s;
}

float  simd_sum(float *X, int n){
  static int ws = sizeof(__m128)/sizeof(float);
  __m128  r128, s128 = {0};
float s = 0.0, sws[ws];
  int m = m % ws;

// add by ws-float from [0 n-m]
for (size_t i = 0; i < n - m; i+=ws) {
r128 = _mm_load_ps(&X[i]);
    s128 = _mm_add_ps(r128, s128);
}
// store result of the above code.
  _mm_store_ps(sws, s128);

// sum sws up, and record it into s.
for (size_t i = 0; i < ws; i++) {
    s += sws[i];
  }

// add the rest, [n-m+1 n-1], by a float.
for (size_t i = 1; i <= m; i++) {
    s += X[n - i];
  }

return s;
}

int main(){
  float *X;
  int  n = 10000;
  int i;
  long et;
  int  ws;
  double r1, r2;
  union cpuid_sig {
    unsigned int  i;
    unsigned char  c[4];
  } sig;
  unsigned int highest, eax, ebx, ecx, edx;

  highest = __get_cpuid_max(0x8000000, &sig.i);
  if(0 < __get_cpuid(highest, &eax, &ebx, &ecx, &edx) && 0 != bit_SSE2 & edx){
    printf("SSE2 is supported.\n");
  } else {
    fprintf(stderr, "SSE2 is NOT supported!\n");
    exit(-1);
  }

  ws = sizeof(__m128)/sizeof(float);
  printf("sizeof(__m128)/sizeof(float) = %ld/%ld = %d (%ld)\n", sizeof(__m128), sizeof(float), ws, sizeof(__m128) - sizeof(float)*ws);

// setup X
  X = calloc(n, sizeof(float));
  for(i=0; i<n; i++){
    X[n] = sinf(0.314F*i);
  }

  et = clock();
  r1 = x87_sum(X, n);
  printf("elapsed time of x87 sum: %f\n", (double)(clock() - et)/CLOCKS_PER_SEC);

  et = clock();
  r2 = simd_sum(X, n);
  printf("elapsed time of SSE2 sum: %f\n", (double)(clock() - et)/CLOCKS_PER_SEC);

  printf("the difference of results: %f\n", fabs(r1 - r2));

  return 0;
}

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関数を関数を用いればよいのですが、建築や美術など見た目以外でも黄金比が実用的に使えるのを実感できます。

黄金分割探索は、微分法を使わずに一変数の最適化問題を解くことができるもので、引数の有限探索区間を3分割し、内側の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")と言うように使います。