餡子付゛録゛

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

最低限知っておくべきGDBの使い方

小ネタ Linux

統合開発環境（IDE）はもちろん、emacsやvimもデバッガーを内部で呼び出して使えるので、アプリケーション開発でデバッガーを直接操作することは少ないと思いますが、使えると便利なときもあるので、最低限の使い方は覚えておきたいところです。詳しいところはマニュアルを見るか、検索するか、生成AIに訪ねたら何か答えてくれる時代ですが、簡単なウォークスルーができる程度の知識を入れておきましょう。なお、ここでの例はC言語ですが、C++やFortranで出力したコードでもデバッグできます。

デバッグ対象のコンパイル

main.cとsub.cの2つのソースコードを用意して、gdbを動かしてみましょう。
main.cは、

#include<stdio.h>
#include<stdlib.h>

int ssum(int n);

int main(int argc, char *argv[]){
        int i;

        for(i = 1; i < argc; i++){
                printf("%d\n", ssum(atoi(argv[i])));
        }

        return 0;
}

とします。
sub.cは、

int ssum(int n){
        int r, i;

        r = 0;
        for(i = 1; i <= n; i++){
                r += i*i;
        }

        return r;
}

とします。
コンパイルは、変数名などのシンボル情報を実行ファイルに含める-gオプションと、最適化が効かないようにする-O0オプションをつけて行うのが一般的です*1。

gcc -o example.exe main.c sub.c -g -O0

今回は試しに、このexample.exeをデバッガーにかけてみます。

èµ·å‹•

コマンドラインで実行ファイル指定

シェルで

gdb example.exe

とうってgdbを起動すれば、example.exeのデバッグが可能になります*2。

コアダンプを残して異常終了したときは*3、例えばコアダンプのファイル名がcore.12345のときは、

gdb example.exe core.12345

とすれば、異常終了時のプロセスの状態を調べることができます*4。

gdbのコマンドで実行ファイル指定

gdbを起動後、

gdb

gdbのコマンドでexample.exeを読み込むこともできます。

exec-file example.exe
symbol-file example.exe

シンボル情報は別に読み込ませる必要があるので、注意しましょう。

コアダンプを残して異常終了したときは、例えばコアダンプのファイル名がcore.12345のときは、

core-file core.12345

と言うように指定します。

環境の変更と確認

例えば環境変数HOMEの変更と確認は、以下のようにします。

set environment HOME = /var/tmp
show environment HOME

作業ディレクトリの変更と確認

例えば/var/tmpに移動する場合は、以下のようにします。

cd /var/tmp
pwd

cdとだけうつと、gdbの起動パスに作業ディレクトリが設定されます。

ソースコードの表示

infoコマンドで、コンパイルに使ったソースコードのファイル名が確認できます。infoはiに省略可能です。

info sources

listコマンドで、ファイル名と行番号もしくは、ファイル名と関数名を指定して、ソースコードの指定箇所前後を表示させることができます。listはlに省略可能で、表示行数を変えたい場合は、第3引数に数字を指定します。

list sub.c:1
info sources; list sub.c:1

関数名の指定は次のようにします。

list sub.c:ssum

実行を中断して変数の状態などを確認し、ステップ実行していくブレークポイントは、ファイル名と行番号もしくは、ファイル名と関数名で設定するので、これらの確認は重要です。

ブレークポイントの設定

break ファイル名:行番号、もしくは、ファイル名:関数名でブレークポイントを設定することができます。
sub.cの中のssum関数でブレークしたい場合は、

break sub.c:ssum

とします。breakはbに省略できます。if文をつけて条件付きブレークを指定することもできます。例えば、

b sub.c:6 if i==2

と言うように書きます。比較演算子はC言語と互換のものが使えます。また、文字列の比較や正規表現を使いたい場合には、gdb独自の便利関数（convenience functions）を使うことができます。

設定したブレークポイントは、

info breakpoints

で確認することができます。
なお、ブレークポイントを削除するときは、info表示の番号をdeleteで指定します。dに省略可能です。

delete 1
info breakpoints and delete a breakpoint

番号を同時に並べて、複数のブレークポイントの同時削除もできます。

実行

runで実行でき、コマンドライン引数を同様に指定できます。runはrに省略可能です。例えばexample.exe 3の実行をデバッグするときは、

run 3

となります。

run 3 > output.txt

とすると標準出力をファイルoutput.txtにパイプでつなぎます。
今回のコードでは意味がないですが、標準エラー出力をつなぐ場合は、

set args 2> output.txt
run 3

とします。

現在位置（バックトレース）の表示

ブレークポイントなどで実行停止している場合、backtrace（btに省略可能）かwhereで現在位置（バックトレース）の表示ができます。

where

例えば、以下の場合では

where

main.cの10行明から呼び出された関数の中の、sub.cの4行目であることがわかります。

バックとレースの表示に番号がついていますが、フレームの番号です。関数の呼び出しで入れ籠構造になっているわけですが、フレームはその構造のどの層かを表します。

変数の表示

例えば変数rを表示する場合は、

print r

とします。printはpに省略可能です。配列aを、print a[1]と一部分だけ表示したり、*ptr@3でポインターptrが指し示す先の頭3つを表示したり、C言語に限らずおなじみのprintf文が使えたりします。

printf "r: %d\n", r

ローカルスコープの変数をまとめて表示する場合は、

info local

で表示されます。また、

info args

で、現在いる関数を呼び出したときの引数を表示します。frameコマンドでフレーム番号を指定するか、upもしくはdownで参照するフレームを変えることができ、info localとinfo argsの対象が変わります。

いちいちprintせずに、ステップごとに変数の状態を確認したい場合は、ウォッチポイントを設定します。

watch r

ブレークポイントと同様にinfo watchpointsで設定一覧を確認し、deleteで番号を指定して削除できます。また、使っている変数が定義されているスコープから抜けた場合、自動的にそのウォッチポイントは解除されます。

変数への代入

例えば変数rに7を代入するときは、

print r = 7

とするか、

set variable r = 7

とします。

set variable

ステップ実行

step

で、1行づつ実行できます。sに省略可能です。関数があった場合、関数の中も一行づつ実行しようとします（step-in）。

next

も、1行づつ実行するコマンドですが、関数の中身はまとめて実行します（step-over）。nに省略可能です。

watch and next

他にfinishで現在の関数の終了まで実行が出来たりします。step-overするつもりで、step-inしてしまったときに便利です。

次のブレークポイントまで実行

continue

で、次のブレークポイントまで実行するか、最後まで実行します。cに省略可能です。最後まで実行した場合、runで再実行できます。

デバッグ中のプロセスの中断

途中でデバッグしているプログラムの実行を中断することができます。

kill

gdbは停止しないので、ブレークポイントの設定などはそのまま残ります。

GDBの終了

quit

で終了できます。qに省略可能です。

まとめ

IDEからgdbを使うことが多いと思いますが、コマンドラインだと操作困難と言うわけでも無いのが分かります。gdbには他にも、停止せずにコマンド実行をさせられるtracepointsや、シグナルや例外に仕掛けられるcatchpointsが設定でき、また、機械語のレベルでもデバッグできます。歴史の長いツールなので、検索するなりすれば色々と資料が出てくるので、必要に応じて使っていきましょう。

*1:無いとgdbでデバッグ不能と言うわけではないですが、困難になります。

*2:引数は後述するrunで指定する他、gdb起動時に--argsの後に指定可能です。

*3:OSがコアダンプを残さない設定になっている場合もあり（e.g. Ubuntu で core ファイルが出力されない #gdb - Qiita ）、その変更には管理者権限が必要になります。

*4:本稿では割愛します

最低限知っておくべきviのコマンド

小ネタ Linux

Linux/UNIXでOS設定ファイルなどをいじるときや、gitなどバージョン管理システムのコミットメッセージの編集などで、なんだかんだとvi（vim-tiny）を使うときが多いわけですが*1、コマンドを忘れるので…と言うか、なんとなく使い続けていたら、バリバリvimを使うわけではないにしろ、これは知っておいた方が良かったなと言うコマンドがそこそこあったのでメモをとっておきます。

本当に最低限、Linux/UNIXの設定でしのぐのに知っておくべきコマンドには❗をつけてあり、全部で7つです。他は、コーディングなどもっと実用的に使うのに（私の主観で）最低限知っておくべきコマンドになります。

最重要

ESCキー コマンドモードに戻る ❗

viはモードエディターで、カーソル移動やコマンド入力を行うモードと、文字を挿入するモードを切り替えて使います。ESCキーでコマンド・モードに戻れることを知らないと、:q!で編集を破棄して終了もできません。

カーソル移動

実際のところviも方向キーが使える*2ので、h（←）j（↓）k（↑）l（→）は覚えなくても良いです。つまり、移動に関しては特に覚えなくてもしのげます。
ただし、

0 （空白を含めた）行頭に移動  
^ （空白を除いた）行頭に移動  
$ （空白を含めた）行末に移動  
G 行末に移動  
nG n行に移動（n∈{a∈ℕ:a≥1}）  

は、便利です。
ソースコードをいじるときは、w（次の単語の先頭）/e（次の単語の末尾）/b（前の単語の先頭） とW（空白区切りで次の文字列の先頭）/E（〃次の文字列の末尾）/B（〃次の文字列の先頭）も覚えておくと便利でしょう。

挿入モードに入る

iもしくはaだけでも良さそうですが、iを知らないと行頭に挿入が、aを知らないと行末に挿入ができません。

i カーソル位置に追加 ❗
a カーソルの次から追加 ❗
I カーソル行の行頭から追加  
A カーソル行の行末から追加  
o カーソル行の下に一行あけて、そこに挿入開始  
O カーソル行の上に一行あけて、そこに挿入開始  

大文字のアルファベットはシフトを同時に押すと理解してください。

コピー&ペースト

コピペできなくても済む場合も多いのですが、覚えておきましょう。

yy カーソル行コピー  
nyy カーソル行とその下n-1行（n∈{a∈ℕ:a≥1}）コピー（eg. 3yyだとカーソル行とその下2行をコピー）  
p カーソル行の下にペースト  
P カーソル行の上にペースト  

yiwで単語のコピー（pはカーソルの下から挿入）、yiWで次の空白か行末までコピー（〃）、:%yで全コピーができたりもできますが、設定ファイルの編集ではあまり使わないでしょう。
なお、viではコピーのことをyankといいます。

削除

xとJだけ知っておけば、だいたい何とかなります。

x カーソル位置1文字削除 ❗
X カーソル前1文字削除  
r カーソル位置1文字上書き  
J カーソル行の改行を消す ❗
dd カーソル行削除  
ndd カーソル行以下n行削除（n∈{a∈ℕ:a≥1}）  

なお、コマンドではありませんが、挿入モードにおいて、その挿入モードに入ってから入力した改行以外の文字を、BackspaceもしくはCTRL + hで1文字づつ消すことができます*3。消しても表示が残ったままで分かりづらいのですが、abcCTRL + hdと入力すると、abdが残ります。
dGでカーソル行以下をすべて削除できたり、dnGでn行以下を削除できたり、dwでカーソル下の単語を削除したりもできますが、設定ファイルの編集ではあまり使わないでしょう。

検索と置換

原理的に必要はないのですが、設定ファイルをいじるときに知らないと困るのが検索です。

/ 下方向に検索  
? 上方向に検索  
n 次を検索  
N 前を検索  

グローバル置換も覚えておくと便利です。

:%s/target/replaced/gc ファイル全体で行のすべてのtargetã‚’replacedに確認しながら置換  

なお、

  • %がないと、カーソル行だけを置換
  • gがないと、各行の最初のtargetだけを置換
  • cがないと、確認がされない

ことになります。

大文字小文字の区分け

グローバル設定を変えることで、アルファベットの大文字小文字の見分けをするか設定できます。

/set ic 大文字小文字の区分けなし  
/set noic 〃 あり  

icはignore caseの略でしょう。

やり直し

vi（vim-tiny）だと辿れる履歴が1回までなのですが、あります。

u 元に戻す  
U 最後に編集した行のすべての変更を元に戻す  

グローバル置換でミスをしたときuは重宝します。

ファイルのオープンとクローズ

設定ファイルは最初から用意されていたり、無くてもtouchコマンドで作れるので

:e filename ファイルを開く  
:r filename ファイルを読み込んで挿入  
:w! 保存 ❗
:q! 破棄して終了 ❗
:wq 保存して終了  

ZZでも保存して終了ができます。

画面分割

意外にviの最初のバージョンからある（と言われる）画面分割です。

:sp 画面上下分割  
:vs 画面左右分発  
CTRL + w + k 上の領域に移動  
CTRL + w + j 下の領域に移動  
CTRL + w + h 左の領域に移動  
CTRL + w + l 右の領域に移動  

分割後、それぞれの領域で独立して:eでファイルを読み込め、それぞれの領域で独立して保存と終了ができます。

表示

高機能な統合開発環境を使うことをお勧めしますが、何らかの事情でソースコードをいじるときは、以下も知っておいた方が便利でしょう。

:set wrap 折返し表示をする（デフォルト）  
:set nowrap 折返し表示をしない  
:set number 行番号を表示する  
:set nonumber 行番号を表示しない（デフォルト）  

*1:GUIがあればgedit、GUIがなくてもFreeBSDだとeeコマンドが、標準インストール状態から用意されている場合も多いのですが、ほとんどの場合でインストール済みなのはviなので。

*2:大昔のADM-3Aと言う方向キーがない端末で動くようになっていたわけですが、少なくともviのクローンは方向キーを受け付けます。

*3:vimだと設定で挙動が変わります

Ubuntu Linux 24のVSCodeでのfortlsの利用方法

小ネタ Fortran Linux

VSCodeでFortranを使うときの定番拡張はModern Fortranなのですが*1、使い出すとすぐに、そして繰り返し、

It is highly recommended to use the fortls to enable IDE features like hover, peeking, GoTos and many more.
For a full list of features the language server adds see: https://github.com/gnikit/fortls

とポップアップでFortranのソースコードを解析してくれるlanguage serverのfortlsをインストールすることを推奨されます。
language serverなしでもシンタックスハイライトは効くので放置していたのですが、いちいち出てくるポップアップが気に障るのでインストールしてみます。

fortlsはPythonで書かれているので、Pythonのパッケージマネージャーであるpipでインストールしますが、Ubuntu LinuxではPythonの仮想環境をつくって、Python仮想環境下でないとインストールと起動ができません*2。

初めてのPythonで "error: externally-managed-environment" を解決する #Python3 - Qiitaを参考に、以下のようにインストールをして

# 仮想環境の作成（~/mypyにファイル配置）
python3 -m venv ~/mypy
# 仮想環境の起動（プロンプトの頭に(mypy)がつく）
source ~/mypy/bin/activate
# インストール
pip3 install fortls

以下のように、仮想環境を起動してから、`code .`とVSCodeを起動すると、

# 仮想環境の起動（プロンプトの頭に(mypy)がつく）
source ~/mypy/bin/activate
# VSCodeの起動
code .

fortlsが動く状態でModern Fortranを使えます。

*1: VSCodeのFortran向け拡張を整理，最新化する（2022年6月） #VSCode-Extension - Qiita

*2:Python仮想環境でないと`error: externally-managed-environment`と言うエラーメッセージが出ます。

最尤法によるパネルデータ分析

R言語

固定効果モデルや変量効果モデルをMCMCでベイズ推定をかける前の準備として、クラスターごとに異なる分散を持つことを仮定した線形モデルを最尤法で推定する関数をつくります。

よくある線形回帰モデルを考えます。

y_{it} = \alpha_i + x_{it} \beta + \tau_t + \epsilon_{it}

yが従属変数、 iが個体（もしくは介入群や対照群などの群）をあらわす添字、 tが時点をあらわす添字、 \alphaが時間に対して不変の固有効果、 xが説明変数のベクトル、 \betaが係数、 \tauが時点ダミー、 \epsilonが誤差項です。 \epsilonは正規分布を仮定して、

\epsilon_{it} \sim \mathscr N\big(0, \sigma_i^2 + \sigma_t^2 \big)

とします。 \sigma_iは個体ごとの分散、 \sigma_tは時点ごとの分散です。

推定コード

固定効果モデルにクラスター頑強標準誤差を用いるよくある手法に近いですが、同様に不均一分散に対応する一方で、系列相関は無いと見なすので劣化しています。クラスター頑強標準誤差と違って、個体の数が少ないときに過小推定バイアスが入ったり、2-wayクラスターで正定値性が保証さ れていない問題*1は起きづらいと思います。
オプションでwithin変換ができるようにして、推定する変数の数を減らせるようにしておきます。個体数が多い場合、計算量を劇的に減らすことができますし、推定結果もパラメーターの数が半分強にになって簡素になります。

fe.lm <- function(formula, cluster = within, data = NULL, within = NULL, sigma.minimum = 1e-5){

    if(is.null(cluster)){
        cluster.matrix <- matrix(1, nrow(data), 1)
    } else {
        expand.cluster <- function(i) model.matrix(as.formula(sprintf("~ 0 + %s", cluster[i])), df01)
        cluster.matrix <- expand.cluster(1)
        if(2<=length(cluster)) for(i in 2:length(cluster)){
            cluster.matrix <- cbind(cluster.matrix, expand.cluster(i))
        }
    }

    mf <- model.frame(formula, data)
    X <- model.matrix(mf, mf)
    y <- model.response(mf)

    # do within transfer
    if(!is.null(within)){
        within_transfer <- function(x){
            for(cn in within){
                i = as.integer(data[[cn]])
                m_a <- mean(x)
                m_i <- tapply(x, i, mean)
                x <- x - m_i[i] + m_a
            }
            c(x, use.names = FALSE)
        }
        X <- apply(X, 2, within_transfer)
        y <- within_transfer(y)
    }

    init.p <- c(rep(1, ncol(cluster.matrix)), rep(0, ncol(X)))

    llf <- function(p){
        sidx <- 1:ncol(cluster.matrix)
        sigma <- p[sidx]
        b <- p[-sidx]
        sum(dnorm(y - X %*% b, sd = sqrt(cluster.matrix %*% sigma^2), log = TRUE))
    }

    llfg <- function(p){
        sidx <- 1:ncol(cluster.matrix)
        sigma <- p[sidx]
        b <- p[-sidx]
        res <- y - X %*% b
        sd <- sqrt(cluster.matrix %*% sigma^2)
        dsigma <- apply(c((-1/sd + res^2/sd^3) / sd) * cluster.matrix, 2, sum)
        dbeta <- apply(X, 2, function(x) sum(x*2*res/sd^2/2))
        c(dsigma, dbeta)
    }

    r_optim <- optim(init.p, llf, llfg, method = "L-BFGS-B", 
        lower = c(rep(sigma.minimum, ncol(cluster.matrix)), rep(Inf, ncol(X))), 
        control = list(fnscale = -1), hessian = TRUE)

    r_optim$se <- sqrt(abs(diag(solve(-r_optim$hessian))))

    if(2 <= ncol(cluster.matrix)) names(r_optim$par) <- c(sprintf("sigma.%s", colnames(cluster.matrix)), colnames(X))
    else names(r_optim$par) <- c("sigma", colnames(X))
    names(r_optim$se) <- names(r_optim$par)

    r_optim
}

利用例

モデル式、データフレーム、独立した分散を持つクラスターを示す変数名、within変換をかける固定効果の変数名を指定します。

テストデータ

今回もデータ生成過程を仮定して生成します。

# obs : number of observatins
# noi : number of individials
# T : maximum number of times
# nc : number of explanatory variables including the constant term
gen_d <- function(obs = 100, noi = 2, T = 3, nc = 3){
    beta <- (-1)^(1 + 1:nc)*(1:nc) # true coefficients

    # variances by cluster
    sigma2_i <- seq(1, 10, length.out = noi)
    sigma2_t <- seq(1, 10, length.out = T)

    # fixed effect
    fe_i <- seq(0, noi - 1, 1)
    fe_t <- seq(0, T - 1, 1)

    # generate explanatory variables
    X <- cbind(1, matrix(runif(obs*(nc - 1), -3, 3), obs))
    colnames(X) <- c("Const.", sprintf("x%d", 2:nc - 1))
    i <- round(runif(obs, 0.5, noi + 0.5))
    t <- round(runif(obs, 0.5, T + 0.5))

    # generate a depending variable vector
    y <- X %*% beta + fe_i[i] + fe_t[t] + rnorm(obs, sd = sqrt(sigma2_i[i] + sigma2_t[t]))

    cnames <- c(letters, apply(expand.grid(letters, letters, stringsAsFactors = FALSE)[, 2:1], 1, paste, collapse = ""))
    data.frame(y, X, i = factor(cnames[i]), t = factor(t - 1))
}

推定

推定量は小標本だと気持ちOLSより改善されていると信じたい程度ですが、ベイズ推定のための前段階なのでよしとします。個体もしくは群ごとに分散が標準誤差つきで得られるので、場合によっては便利でしょう。

# データセットを作成する
set.seed(20241209)
df01 <- gen_d(50)

# 変数iとtでwithin変換をかける（自動でcluster = c("i", "t")が指定される）
r_within <- fe.lm(y ~ x1 + x2, data = df01, within = c("i", "t"))

# 変数iでwithin変換をかけ、iとtによるクラスターにそれぞれ独立した撹乱項をつける
# r_within <- fe.lm(y ~ x1 + x2 + t, cluster = c("i", "t"), data = df01, within = c("i")) 

# with変換をかけずにダミー変数を入れ、iとtによるクラスターにそれぞれ独立した撹乱項をつける
r_dummies <- fe.lm(y ~ x1 + x2 + i + t, cluster = c("i", "t"), data = df01)

# r_no_cluster <- fe.lm(y ~ x1 + x2 + i + t, data = df01)

# 比較用のOLSによる推定
r_lm <- lm(y ~ x1 + x2 + i + t, df01)

# 比較用のestimatrパッケージによる推定（Stata互換のクラスター頑強標準誤差が計算できる）
library(estimatr)
r_estimatr <- lm_robust(y ~ x1 + x2 + i + t, data = df01, clusters = i, se_type = "stata")

# 推定結果の係数
m.coef <- matrix(round(c(1, -2, 3, 1, 1, 2, 
    coef(r_lm), coef(r_estimatr), 
    r_dummies$par[-(1:5)],
    r_within$par[6:8], NA, NA, NA), 3), 
    5, byrow = TRUE, dimnames = list(c("TRUE", "OLS", "estimatr", "ML/dummies", "ML/within"), names(coef(r_lm))))
cat("coefficients:\n")
print(m.coef)



# 推定結果の標準誤差
m.std.error <- matrix(round(c(
    coef(summary(r_lm))[,2], 
    coef(summary(r_estimatr))[,2], 
    r_dummies$se[-(1:5)],
    r_within$se[6:8], NA, NA, NA
    ), 3), 
    4, byrow = TRUE, dimnames = list(c("OLS", "estimatr", "ML/dummies", "ML/within"), names(coef(r_lm))))
cat("std.error:\n")
print(m.std.error)

*1:詳細は太田 (2013)を参照してください。なお、どちらも改善する手法が提案されており、致命的というわけではないです。

Ubuntu LinuxにTTFフォントを追加し、Rで使えるようにする

Linux R言語 小ネタ

何年か前に話題になっていたHackGenフォントã‚’Ubuntu LinuxのRで使えるようにしたので、手順を記録しておきます。

  • 最新リリースになる HackGen_v2.9.0.zip をダウンロードしてきます。
  • sudo -sをしてrootになります。
  • ディレクトリ /usr/local/share/fonts/truetype/HackGen/ を作成します*1。
  • 作成したディレクトリに、HackGen_v2.9.0.zip の中身の HackGen35ConsoleNF-Bold.ttf HackGen35ConsoleNF-Regular.ttf HackGenConsoleNF-Bold.ttf HackGenConsoleNF-Regular.ttf をコピーします。
  • fc-cache -fv でフォントリストを更新し、fc-list | grep HackGen でフォントリストに載ったかを確認します。

これでpngなどでHackGenが使えるようになりました。Linux/UNIXのRではpar(family = "HackGen35 Console NF")という風に指定できます*2。

EPSで使えないと不便なので、さらにghostscriptのcidfmapを編集します。Ubuntu 24のaptでインストールしたバージョンでは、/var/lib/ghostscript/fonts/cidfmap に

/HackGen35 << /FileType /TrueType /Path (/usr/local/share/fonts/truetype/HackGen/HackGen35ConsoleNF-Regular.ttf) /SubfontID 0 /CSI [(Japan1) 5] >> ;
/HackGen35-Bold << /FileType /TrueType /Path (/usr/local/share/fonts/truetype/HackGen/HackGen36ConsoleNF-Bold.ttf) /SubfontID 0 /CSI [(Japan1) 5] >> ;

の2行を追加します。(Japan1) 5の5は、既存の(Japan1) 数字と被らないようにしているだけなので、適時、変更してください。これで、

if(!"HackGen35" %in% names(postscriptFonts())){
    postscriptFonts(HackGen35 = CIDFont("HackGen35", "UniJIS-UCS2-H", "UCS-2BE"))
}
fname <- "example_HackGen35.eps"
postscript(fname, family = "HackGen35") # parでは指定できない
n <- 100
x <- runif(n)
y <- 1 + x + rnorm(n)
plot(y ~ x, main = "An example plot")
dev.off()
# フォントを埋め込んだファイルをつくる
embedFonts(fname,
    outfile = gsub(".eps$", "-embedded.eps", fname),
    options = "-c \"<</NeverEmbed []>> setdistillerparams\" -f ")

が出来るようになりました。

*1:~/.fonts/HackGen/ でも良く、この場合はrootにならなくてよいです

*2:WindowsのRの場合は、windowsFontsの登録が事前にいります

最尤法による欠損値の補定と回帰モデルの推定の同時実行

R言語

ここ10年間ぐらいで多重代入法を用いた欠損データ処理が一般化しましたが、多重代入法を使わなくても多重代入法と同様の結果が得ることもできます。Fully Bayesian Imputation Methodと呼ばれている手法ですが*1、最尤法による単一代入（imputation model）と回帰分析（analysis model）の同時推定するだけで、頻度主義的にも扱えるので試してみましょう*2。

データセット

Listwise手法ではバイアスが入るように、MARなデータセットを作成します。

gen_data <- function(n, mu = c(2, -3), S = matrix(c(2, -1.5, -1.5, 2), 2, 2)){
    library(mvtnorm)
    X <- rmvnorm(n, mu, S)
    beta <- c(1, -1, 1)
    y <- cbind(1, X) %*% beta + rnorm(n)
    x <- X[,1]
    z <- X[,2]
    p <- 1/(1  + exp(-(4 + y)))
    mx <- x 
    mx[runif(n) < p] <- NA
    data.frame(y, x, mx, z)
}

set.seed(1029)
df01 <- gen_data(50)

yが従属変数、xとzが欠損値なしの説明変数です。cov(x, z)=-1.5≠0なので、zからx、xからzが推定できます。xに欠損値をつくったのがmxになり、今回、推定で用いるのはこれです。xはベンチマークとしてのみ使います。切片項、x、zの係数はそれぞれ1、-1、1になります。

データ整理

データフレームから、欠損値も入った行列X_aと、欠損値を除いた行列X_c、欠損値だけの行列X_mの三つを作成します。

frml <- y ~ mx + z
df02 <- model.frame(frml, data = df01, na.action = na.pass)
y <- model.response(df02)
X_a <- model.matrix(frml, df02)
X_c <- X_a[complete.cases(X_a), ]
X_m <- X_a[!complete.cases(X_a), ]

推定

最尤法をかけるわけですが、

  • 線形回帰のための最尤法の目的関数の中で、補定処理のための最尤法を行う入れ籠構造になっている
  • 補定は多変量正規分布を仮定した最尤法（実は複数の欠損に対応）
  • 線形回帰のパラメーターになるσとβの他、補定に用いるパラメーターμ（多変量正規分布の平均）を同時推定
  • 補定に用いるΣ（多変量正規分布の共分散行列）は、μとデータから作成（concentrated loglikelihood function）
  • 1è¡Œ1行は補定せず、全部の行をまとめて補定している（→大きなデータセットには対応できない）

ことに注意してください。

objf <- function(p){
    sigma <- p[1]
    beta <- p[2:4]
    mu <- p[5:6]
    X0 <- apply(X_c[, -1], 1, \(x) x - mu) # 切片項は補定しないので、除外する
    Sigma <- X0 %*% t(X0) / (ncol(X0) - 1)
    invSigma <- chol2inv(chol(Sigma))
    init.p <- rep(0, sum(is.na(X_a)))
    r_optim <- optim(init.p, function(p){
        X_a[is.na(X_a)] <- p
        sum(dmvnorm(X_a[, -1], mu, Sigma, log = TRUE))
    }, function(p){ # グラディエントは無くても動く
        na <- is.na(X_m)
        index <- 1:length(X_m)
        X_m[na] <- p
        Y <- X_m[, -1] - rep(mu, each = nrow(X_m))
        g <- t(-invSigma %*% t(Y))
        g[index[na] - nrow(X_m)] # gは切片項を除外した行列からの計算結果なので、-nrow(X_m)で位置をあわせる
    }, method = "BFGS", control = list(fnscale = -1))
    X_a[is.na(X_a)] <- r_optim$par
    r <- sum(dmvnorm(X_c[,2:3], mu, Sigma, log = TRUE))
    r <- r + sum(dnorm(y - X_a %*% beta, sd = sigma, log = TRUE))
    cat(sprintf("convergence: %d, loglikelihood: %f", r_optim$convergence, r), "\n")
    r
}

init_p <- c(1, rep(1, 3), apply(df02, 2, mean, na.rm = TRUE)[2:3])
r_mnorm <- optim(init_p, objf, method = "L-BFGS-B", lower = c(1e-1, rep(-Inf, 5)), control = list(fnscale = -1, maxit = 1000), hessian = TRUE)

そこそこ待たされます。

比較のための推定

OLSとmiceでも推定しておきます。

r_lm_c <- lm(y ~ x + z, df01)
r_lm_m <- lm(y ~ mx + z, df01)

library("mice")
df_mice <- mice(df01[,-2], m = 20, maxit = 50, method = "pmm", seed = 2024)
miced_lm <- with(df_mice, lm(y ~ mx + z))
r_mice <- pool(miced_lm)

推定結果

表にまとめます。

推定量の比較

m <- matrix(c(coef(r_lm_c), coef(r_lm_m), r_mnorm$par[2:4], summary(r_mice)[, 2]), 4, byrow = TRUE)
rownames(m) <- c("OLS/no missing", "OLS/list-wise", "ML/imputed", "mice/imputed")
colnames(m) <- sprintf("beta %d", (1:ncol(m))-1)
print(m)
推定量の比較

今回の手法（ML/imputed）が、欠損値のある行を削除するOLS/Listwise手法よりも良い推定値を出していることが分かります。

完全データのOLSとの差異

d <- t(sapply(2:nrow(m), \(i) m[i, ] - m[1, ]))
rownames(d) <- rownames(m)[-1]
print(d)
欠損値無しデータのOLS推定量との差異

完全データの推定量との差分をとってmiceの結果と比較すると、切片項以外は今回の手法（ML/imputed）の方がmice推定量よりも完全データの推定量に近い推定量となっています。

標準誤差の比較

m.se <- matrix(c(
    coef(summary(r_lm_c))[,2],
    coef(summary(r_lm_m))[,2],
    sqrt(diag(solve(-r_mnorm$hessian)))[2:4],
    summary(r_mice)[,3]
    ), 4, byrow = TRUE)
colnames(m.se) <- colnames(m) 
rownames(m.se) <- rownames(m)
print(m.se)
標準誤差の比較

単一代入だと標準誤差が過少推定される問題があるのですが、今回の手法（ML/imputed）の方がmice推定量よりも大きいので、過剰に帰無仮説を棄却する恐れは少ないことが分かります。

なぜ最尤法を使うのか？

多くの場合はmiceで十分間に合い、信頼のおけるパッケージと利用実績からmiceを使う方が望ましいのですが、この手法にも利点があります。ベイズ統計学への応用を考えると多重代入法は煩雑になりますが、尤度関数ですべて処理できれば簡単になります。また、補定処理に使うパラメーターも点推定されるので、実際に最尤推定量で代入された値が特定できるのは、プロットに使うときなどで便利なことが多いと思います。
なお、処理時間がかかりすぎるので使う気になれない場合は、補定処理を多変量正規分布の最尤推定から線形回帰の予測になどに書き換えれば時間短縮になります。複数の欠損値に対応できなくなりますが、実用上は特定列の欠損値を埋められれば足りることが多いでしょう。

*1:Best and Mason (2013) "Bayesian methods for missing data: part 1 — Key Concepts"を参考にしました。

*2:2年前に試しに書いたコード中でも使っていたりするのですが、今回は欠損値補定の部分をもっと一般化してみます。

R/FortranでHamiltonian Monte Carlo法によるベイズ推定を書いてみたよ！

R言語 Fortran

マルコフ連鎖モンテカルロ法（MCMC）には様々なアルゴリズムがあり、メトロポリス・ヘイスティングス法は尤度関数があれば推定ができる使い勝手のよい方法だと知られています。しかし、高次元で効率が悪いことが知られており、また多峰の確率分布には上手く対応できないそうで、最近のベイズ推定ではHamiltonian Monte Carlo法や、その発展であるNo-U-Turn Samplerが使われています。

Hamiltonian Monte Carlo法は教科書的な手法のようで、Rにもhmclearnと言うHamiltonian Monte Carlo法の教育用パッケージがあるのですが、参考にしつつFortranで書いてみました。推定モデルにLinear Mixed EffectsがなくOrdered Logitがあったり、初期値とステップ幅を最尤法の結果から設定するようにするなど*1、細部では変更があります。

https://github.com/uncorrelated/HMCFortran

Linux環境でgit cloneしたディレクトリでR CMD INSTALL . でRにインストールして使うようになっていますが、Fortranのライブラリを2つ（と片方に内蔵されたCのライブラリを1つ）をリンクします。つまり、lbfgsbとdSFMT_F03_interfaceを、どこかにディレクトリー（e.g. /var/tmp）を作ってgit cloneして、MakefileのSHSRC=/var/tmpを適時書き換える必要があります*2。

cd /var/tmp
git clone https://github.com/uncorrelated/HMCFortran.git
git clone https://github.com/jacobwilliams/lbfgsb.git
git clone https://github.com/jsspencer/dSFMT_F03_interface.git
cd HMCFortran
R CMD INSTALL .

使い方は以下のように、なるべくR的にしました。モデル式を書いて、summaryやplotなどのジェネリック関数で結果が見られるようにしています。

library(hmcfortran)

# data generator for practice
mkdata <- function(nr, nc, sd = 1){
	X <- matrix(c(rep(1, nr), runif(nr*(nc-1), -2, 2)), nr)
	colnames(X) <- c("Const.", letters[1:(nc - 1)])
	beta <- c(0.5, rep(1, nc - 1))*(-1)^(1:nc)
	y <- X %*% beta + rnorm(nr, sd = sd)
	cbind(y, as.data.frame(X))
}

np <- 4 # number of parameters
data <- mkdata(150, np - 1)

# prior/hyper parameters
ig.alpha <- 1
ig.gamma <- 1
beta.mean <- rep(0, np - 1)
beta.Sigma <- diag(1e+10, np - 1, np - 1)

# do Hamiltonian Monte Carlo method to estimate a liner model
r_hmc_lm <- hmc.blm(y ~ a + b, ig.alpha, ig.gamma, beta.mean, beta.Sigma, data = data, N = 3000, nchains = 3)

# get the marginal likelihood
mlf_1 <- log_ml(r_hmc_lm)

# alternative model: third parameter, a, is zero.
r_hmc_lm_0 <- hmc.blm(y ~ b, 1, 1, rep(0, np - 2), diag(1e+10, np - 2, np - 2), data = data, N = 3000)

# get the marginal likelihood of alternative model
mlf_0 <- log_ml(r_hmc_lm_0)

cat("Laplace Approximated BF_{10}:", exp(mlf_1 - mlf_0), "\n")

BF <- exp(log_BF(r_hmc_lm, 3, 0))
cat("IWAMDE Savage-Dickey ratio BF_{01}:", BF, " BF_{10}:", 1/BF, "\n")

# show the prediction, where explanatory variables are the mean of X.
prdct <- predict(r_hmc_lm, X = apply(r_hmc_lm$input$X, 2, mean))
hist(prdct, breaks = 20)

# call coda::plot
plot(r_hmc_lm)

# calcurate R-hat
gelman.diag(r_hmc_lm)

# call other functions of coda package
# gelman.plot(r_hmc_lm$mcmc.list)
# geweke.plot(r_hmc_lm$mcmc.list)
# crosscorr.plot(r_hmc_lm$mcmc.list)
# autocorr.diag(r_hmc_lm$mcmc.list)
# autocorr.plot(r_hmc_lm$mcmc.list)

# call coda::summary
summary(r_hmc_lm)

# call coda::HPDinterval
HPDinterval(r_hmc_lm)

# statistics of the sample of MCMC
mean(r_hmc_lm)
median(r_hmc_lm)
quantile(r_hmc_lm, 0.5)


線形回帰、二項ロジット、順序ロジット、混合ロジット、ポアソン回帰の他、Rのユーザー定義の目的関数を使えるようにしています。プリセットのモデルでは周辺尤度のラプラス近似と、MCMCサンプルから制約モデルのベイズファクター（IWMDE/Chen (1992) Savage-Dickey density ratio）が計算でき、codaパッケージを利用することでMCMCの収束判定も教科書的に可能なので、MCMCpackと同程度に気軽にベイズ推定ができます。ドキュメントがまったく整っていないので、何とかしたいですね！

今後、気力と体力があればLinear Mixed Effectsと欠損値処理は実装したいのですが、どうも無さそうなので期待しないでください('-' )＼(--；)BAKI

*1:Find Reasonable EpsilonやDual Averagingは試したのですが、かなりの頻度で異常な推定結果がもたらされるので諦めました。

*2:ライセンス的には、まとめて配布しても問題ないのですが、まだまだ開発途上なので分けています。