FASTGLを使ってC++で積分

数値解析における数値積分法の代表としてガウス求積と言う技法があり、それに必要なガウス点と重みの計算でJuliaがC++より桁違いに速いと言う話が流れていたのですが、Juliaで使われているFastGaussQuadrature.jlが使っているアルゴリズムのC++版を使ったところ、案の定、だいたい同じの速度になりました。

Fast computation of Gauss-Legendre Quadrature Points and Weights（FASTGL）のアルゴリズム開発者のページからfastgl.cppとfastgl.hppをダウンロードしてきて、積分に使うだけです。サブディレクトリfastglにこの2つのファイルを配置したところで、以下のソースコードのファイルを作成してコンパイルしたら簡単に実行できます。-O3 -lmオプションはいると思いますが。

#include <iostream>
#include <chrono>
#include <iomanip>
#include <functional>
#include "fastgl/fastgl.hpp"

/*
  FASTGLを使って積分をする関数
  args: ラムダ関数f(x), xの最小値, xの最大値, 節点の数
*/
double integral(std::function<double(double)> f, double s, double e, int NumberOfKnots){
  double Res = 0.0;
  double xm = 0.5 * (s + e);
  double xr = 0.5 * (e - s);

  for ( int k = 1 ; k <= NumberOfKnots ; ++k ){
    fastgl::QuadPair p = fastgl::GLPair ( NumberOfKnots, k );
    Res += p.weight * f( xm + xr * p.x());
  }

  return Res*xr;
}

int main(){
  std::chrono::system_clock::time_point stime, etime;
  stime = std::chrono::system_clock::now();

  // boost::math::quadratureあたりと異なり、節点が100万個でも大した時間にならない
  double Res = integral([](const double& t) { return 1.0 / 2 / std::sqrt(t); }, 1, 4, 1000000);

  etime = std::chrono::system_clock::now();

  std::cout << std::setprecision(16);
  std::cout << "Results: " << Res << std::endl;
  std::cout << "Elapsted Time: " << std::chrono::duration_cast<std::chrono::milliseconds>(etime - stime).count() << " ms" << std::endl;

  return 0;
}

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

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

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法の概要

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

になります。ここで、二階微分の計算があると上手く計算されないことが多々あるので、一階微分しかないフィッシャー情報行列の符号を反転させたもので代用しようと言うのが、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のですが、行列がとなる分散共分散行列なので素朴に計算できます。これを観測数で割ったもの*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. ステップ幅の計算

ステップ幅の具体的な計算方法については原論文やアルゴリズム紹介には言及が無いので、「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 = -1*solve( 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

分散共分散行列は･･･もっとサンプルサイズが要るのでしょう。

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

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

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

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

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

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

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

Wilson score intervalを使う

これが無難です。Rのbinomパッケージでサポートされています。

if(!any(suppressWarnings(library(quietly=TRUE, verbose=FALSE)$results[,"Package"] == "binom"))){
  stop("Do install.packages(\"binom\") before runnning this script.")
}
library(binom)
binom.confint(2, 1971, conf.level=0.95, methods="wilson")

シミュレーションによる信頼区間の推定

厳密な定義はややこしい信頼区間ですが、帰無仮説でパラメーターの真の値としても、統計的仮説検定で帰無仮説が棄却されないことになる値の区間になります。パラメーターから乱数でデータを生成すると、その上下の裾野が棄却域になるわけですが、裾野ギリギリのところに観測値が来るパラメーターが、信頼区間の境界点となるので、愚直に計算することも可能です。

obs <- 1971
occ <- 2
p <- occ/obs

n <- 1000
a <- 0.05

seed <- 849

# 真のパラメーターがqのときの両側a%の棄却域を求める
outskirt <- function(q){
    set.seed(seed)
    x <- rbinom(n, obs, q)
    estimated <- numeric(obs + 1)
    th0 <- 0
    th1 <- obs
    s <- 0
    for(i in 0:obs){
        estimated[i + 1] <- sum(x==i)
        s <- s + estimated[i + 1]
        if(s < n*a/2){
            th0 <- i + 1
        } else if(th1==obs && s >= n*(1-a/2)){
            th1 <- i 
        }
    }
    c(th0/obs, th1/obs)
    # qbicomを使って出した値で区間を求めると、なぜかbinom.testのP値がa/2ギリギリの値にならない
    # c(qbinom(a/2, obs, q), qbinom(1-a/2, obs, q))/obs
}

r_th0 <- optimize(function(q){
    abs(p - outskirt(q)[2])
}, c(0, p), tol=1e-12)

r_th1 <- optimize(function(q){
    abs(outskirt(q)[1] - p)
}, c(p, min(1, 5*p)), tol=1e-12)

cat(sprintf("%.0f%% CI (%.6f, %.6f)\n", 100*(1-a), r_th0$minimum, r_th1$minimum))

CI (0.000148, 0.003340)とWilson score intervalやClopper and Peason intervalと比較して広めになり、この両端の値を帰無仮説に観測値でbinom.testで統計的仮説検定をかけても棄却にはならないです。棄却域ギリギリのところまでにしているので注意してください。

図を描くとB(1971, 0.000148)の右側の裾野の境界、B(1971, 0.003340)の左側の裾野の境界に、観測数2が来ていることが分かります。B(1971, 0.000388)の0から2までの確率質量の合計は0.975以上、B(1971, 0.003340)の2から1971までの確率質量の合計も0.975以上となります。

等裾事後信用区間

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

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

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

prior <- function(p){
# Beta(1/2, 1/2)分布を事前分布に置くとJeffreys intervalという立派な名前になる
#  dbeta(p, 1/2, 1/2)
# 1を置いておくと、事前分布は一様分布[0 1]になる
    1
}

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

# 確率密度関数
pdf <- function(p){
    r <- numeric(length(p))
    for(i in 1:length(p)){
        # 事前分布にBeta(1/2, 1/2)を置いていることに注意
        r[i] <- dbinom(m, n, p[i])*prior(p[i])/s
    }
    r
}

# 累積分布関数
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))

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)

最高事後密度信用区間

最高事後密度信用区間（HPD: highest posterior density）で計算しろって言われたので追記。広く使われている両側の裾野の面積が同じになる信頼区間の代わりにしたら混乱がありそうですが。

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

# 確率密度関数
ppdf <- function(p){
    r <- numeric(length(p))
    for(i in 1:length(p)){
        r[i] <- dbinom(m, n, p[i])/s
    }
    r
}

# 累積分布関数
# minpdf <= pdf(x)のpdf(x)の値だけ積分する
cdf <- function(minpdf){
    integrate(function(a){
        v <- ppdf(a)
        ifelse(minpdf <=v, v, 0)
    }, 0, 1)$value
}

# Highest Posterior Densityで95%信用区間のpdfの閾値を求める
# pdf(x)が一定以上のxの点だけHPDの領域に入る
a <- 0.025
r_optimize <- optimize(function(minpdf){ abs(cdf(minpdf) - (1 - 2*a)) }, c(0, ppdf(p)))
mpdf <- r_optimize$minimum

# Highest Posterior Densityは多峰に対応するというか、区間が分割される可能性がある
# chgpの奇数アドレスに開始点を、偶数アドレスに終了点を詰んでいく
xlim <- c(0, 6*p)
chgp <- numeric()
x <- seq(xlim[1], xlim[2], length.out=100) # プロットのピンク部が粗い場合はグリッドサイズを増やす
y <- ppdf(x)
flag <- FALSE
for(i in 1:length(x)){
  if(flag){
    if(y[i] <= mpdf){
      flag <- FALSE
      chgp <- c(chgp, i)
      # グリッドから計算すると低精度なので、境界値を高精度に求める
      if(1<i){
        r_optimize <- optimize(function(p){ abs(mpdf - ppdf(p)) }, c(x[i-1], x[i]), tol=1e-12)
        # x[i]とy[i]を更新
        x[i] <- r_optimize$minimum
        y[i] <- ppdf(x[i])
      }
    }
  } else {
    if(y[i] >= mpdf){
      flag <- TRUE
      chgp <- c(chgp, i)
      # グリッドから計算すると低精度なので、境界値を高精度に求める
      if(1<i){
        r_optimize <- optimize(function(p){ abs(mpdf - ppdf(p)) }, c(x[i-1], x[i]), tol=1e-12)
        # x[i]とy[i]を更新
        x[i] <- r_optimize$minimum
        y[i] <- ppdf(x[i])
      }
    }
  }
}

j <- 1
while(j < length(chgp)){
  th0 <- x[chgp[j]]
  th1 <- x[chgp[j + 1]]
  cat(sprintf("%.0f%% HPD CI (%.5f, %.5f)\n", 100*(1-2*a), th0, th1))
  j <- j + 2
}

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

*2:Bernstein-von Mises定理からだいたい同じものと捉えていますが、厳密には信用区間と呼ぶべきです。

Rと手作業で覚えるGLM

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

1. GLMとは

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

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

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

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

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

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

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

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

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

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

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

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

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