【R/Rcpp】Rcppを使って項目反応理論 (IRT) のスクラッチを高速化
タイトルの通り、Rcppで項目反応理論をスクラッチ実装しました。前回記事の追記で、ほぼすべての関数をC++で書き換えて高速化しました。
sunaninattahito.hatenablog.com
【追記: 2022/02/14】RcppAramdilloを使って行列計算をちゃんとしたところ、推定時間を半分以下に抑えられるようになりました。また、ヘッダーファイルを作成するなどして、コードファイルを分散させて整理しました。作成したコードや使い方などをGitHub上で公開してますので、より詳しくみたい方はそちらもチェックしてください。
おさらい
生成モデル ~ Metropolis-Hastings Algorithmまでザっとおさらいします。より詳しいことは前回の記事に載ってるのでそちらを参照してください。
個人 $i = 1,...,I$ が 問題 $j = 1,...,J$ にて正解する ($y_{ij} = 1$; otherwise 0) 確率を以下で与えます:
\begin{align*}
p_{ij} = \text{Pr}(y_{ij}=1\ |\ \alpha_j, \beta_j, \theta_i) = \text{logit}^{-1}(\beta_j \theta_i - \alpha_j).
\end{align*}
ここで、$\alpha_j$, $\beta_j$はそれぞれ問題$j$の難易度 (difficulty) と識別度 (discrimination)、$\theta_i$は個人$i$の能力 (ability) で、$\text{logit}^{-1}$は逆ロジット関数です。Likelihood $f(\cdot\ |\ \cdot)$は
\begin{align}
f(y\ |\ \alpha, \beta, \theta) &= \prod\limits_{i=1}^{I}\prod\limits_{j=1}^J \left[p_{ij}^{y_{ij}} (1-p_{ij})^{1 - y_{ij}}\right]\\
\end{align}
です。次に、各パラメータの事後分布を考えていきます。事前分布として、以下を与えます:
\begin{align*}
\alpha_j &\sim \mathcal{N}(a_0, A_0),\\
\beta_j &\sim \mathcal{N}(b_0, B_0), \\
\theta_i &\sim \mathcal{N}(0, 1).
\end{align*}
ベイズの定理より、その他のパラメータを所与としたうえでの各パラメータのconditionalな事後分布は、それぞれ以下で与えられます:
\begin{align*}
\pi(\alpha_j\ |\ y, \beta_j, \theta_i) &\propto f(y | \alpha, \beta, \theta)\pi(\alpha)\\
&\propto \prod\limits_{i=1}^{I} \left[p_{ij}^{y_{ij}} (1-p_{ij})^{1 - y_{ij}}\right] \cdot
N(\alpha_j \ |\ a_0, A_0) \forall j = 1,\dots J,\\
\pi(\beta_j\ |\ y, \alpha_j, \theta_i) &\propto f(y | \alpha, \beta, \theta)\pi(\beta)\\
&\propto \prod\limits_{i=1}^{I}\left[p_{ij}^{y_{ij}} (1-p_{ij})^{1 - y_{ij}}\right] \cdot
N(\beta_j \ |\ b_0, B_0) \forall j=1,...J,\\
\pi(\theta_i\ |\ y, \alpha_j, \beta_j) &\propto f(y | \alpha, \beta, \theta)\pi(\theta)\\
&\propto \prod\limits_{j=1}^J \left[p_{ij}^{y_{ij}} (1-p_{ij})^{1 - y_{ij}}\right] \cdot
N(\theta_i \ |\ 0, 1) \forall i=1,...,I.
\end{align*}
次に、Metropolis-Hastings Algorithmに必要なサンプルの採択確率等を説明します。ノテーションの簡略化のために任意のパラメータを$\delta_k$と定めます ($\delta_{1} = \alpha$, $\delta_2 = \beta$, $\delta_3 = \theta$)。まず、各パラメータの候補$\delta_k^*$を以下のようなランダムウォークな分布からサンプリングします。
\[\delta_k^* \sim N(\delta_k^{(t-1)}, \kappa_\delta)\]
$\delta_k^{(t-1)}$は1イテレーション前の$\delta_k$の値であり、$\kappa_\delta$はMH法でいうところのチューニングパラメータです。また、採択確率 (acceptance probability) は以下で計算します:
\[
ap = \min\left\{\frac{f\left(\delta_k^{*}\ | \ y, \delta_{-k}^{(t-1)}\right) g\left(\delta_k^{(t-1)}\ |\ \delta_k^{*} \right)}{f\left(\delta_k^{(t-1)}\ | \ y, \delta_{-k}^{(t-1)}\right) g\left(\delta_k^*\ |\ \delta_k^{(t-1)} \right)}, 1\right\}
\]
$g(\cdot |\ \cdot)$は提案分布です。また、とあるrandom number $u$を$U(0, 1)$からサンプリングして、$u < ap$ならばパラメータのサンプル候補である$\delta_k^{*}$を採択します。そうでない場合は、$\delta_k^{(t-1)}$をサンプルとして採択します。
上記の採択確率を対数に変換します。すなわち
\begin{align*}
\log ap = \min&\{\log f\left(\delta_k^{*}\ | \ y, \delta_{-k}^{(t-1)}\right) + \log g\left(\delta_k^{(t-1)}\ |\ \delta_k^{*} \right) \\
&- \left[\log f\left(\delta_k^{(t-1)}\ | \ y, \delta_{-k}^{(t-1)}\right) + \log\left(\delta_k^*\ |\ \delta_k^{(t-1)} \right)\right], 0\}
\end{align*}
であり、$\log(u) < \log(ap)$であれば$\delta_k^*$を採択します。そうでない場合は、$\delta_k^{(t-1)}$をサンプルとして採択します。これらのステップを任意のイテレーション数まで繰り返します。
Rcppでサンプラーを書く
すべてのパラメータのサンプリングに必要な関数などをirt_sampling.cppというファイルに一まとめにしました。
#include <Rcpp.h> using namespace Rcpp; // setseed function for replicability // [[Rcpp::export]] void set_seed(int seed) { Rcpp::Environment base_env("package:base"); Rcpp::Function ss = base_env["set.seed"]; ss(seed); } // Log-likelihood funcion // [[Rcpp::export]] NumericMatrix log_likelihood(NumericMatrix Y, NumericVector alpha, NumericVector beta, NumericVector theta) { int I = theta.length(); int J = alpha.length(); NumericMatrix prob(I, J); NumericMatrix log_p(I, J); for (int j = 0; j < J; j++) { for (int i = 0; i < I; i++) { prob(i, j) = std::exp(beta[j] * theta[i] - alpha[j]) / (1 + std::exp(beta[j] * theta[i] - alpha[j])); log_p(i, j) = Y(i, j) * std::log(prob(i, j)) + (1 - Y(i, j)) * std::log((1 - prob(i, j))); } } return(log_p); } // SAMPLING ALPHA // [[Rcpp::export]] NumericVector alpha_sampling(NumericMatrix Y, NumericVector alpha_old, NumericVector beta_old, NumericVector theta_old, double a0, double A0, double MH_alpha, int seed) { set_seed(seed); int J = Y.ncol(); // # of alpha NumericMatrix alpha_star_mat(J, 1); // candidate sample of alpha NumericMatrix log_prop_star_mat(J, 1); // candidate alpha proposal density (as matrix for convenience) NumericMatrix log_prop_old_mat(J, 1); // old alpha proposal density (as matrix) // Sampling from random-walk normal prior for (int j = 0; j < J; j++) { alpha_star_mat(j, _) = Rcpp::rnorm(1, alpha_old[j], MH_alpha); } // Convert into vector NumericVector alpha_star_vec = alpha_star_mat(_, 0); // Generate proposal density for (int j = 0; j < J; j++) { log_prop_star_mat(j, _) = Rcpp::dnorm(NumericVector::create(alpha_star_vec[j]), alpha_old[j], MH_alpha, true); log_prop_old_mat(j, _) = Rcpp::dnorm(NumericVector::create(alpha_old[j]), alpha_star_vec[j], MH_alpha, true); } // Convert into vector NumericVector log_prop_star_vec = log_prop_star_mat(_, 0); NumericVector log_prop_old_vec = log_prop_old_mat(_, 0); // Posterior density NumericVector log_cc_star = colSums(log_likelihood(Y, alpha_star_vec, beta_old, theta_old), true) + Rcpp::dnorm(alpha_star_vec, a0, A0, true); NumericVector log_cc_old = colSums(log_likelihood(Y, alpha_old, beta_old, theta_old), true) + Rcpp::dnorm(alpha_old, a0, A0, true); // Calculate acceptance probability NumericVector log_ap = pmin(log_cc_star + log_prop_old_vec - log_cc_old - log_prop_star_vec, 0); // Save samples (if log(u) < log(ap), accept alpha_star, otherwise alpha_old) NumericVector log_u = log(Rcpp::runif(J, 0, 1)); NumericVector sample(J); for (int j = 0; j < J; j++) { if (log_u[j] < log_ap[j]) { sample[j] = alpha_star_vec[j]; } else { sample[j] = alpha_old[j]; } } return(sample); } // SAMPLING BETA // [[Rcpp::export]] NumericVector beta_sampling(NumericMatrix Y, NumericVector alpha_old, NumericVector beta_old, NumericVector theta_old, double b0, double B0, double MH_beta, int seed) { set_seed(seed); int J = Y.ncol(); // # of beta NumericMatrix beta_star_mat(J, 1); // candidate of beta (as matrix for convenience) NumericMatrix log_prop_star_mat(J, 1); // candidate alpha proposal density (matrix) NumericMatrix log_prop_old_mat(J, 1); // old alpha proposal density (matrix) // Sampling from random-walk normal prior for (int j = 0; j < J; j++) { beta_star_mat(j, _) = Rcpp::rnorm(1, beta_old[j], MH_beta); } // Convert into vector NumericVector beta_star_vec = beta_star_mat(_, 0); // Generate proposal density for (int j = 0; j < J; j++) { log_prop_star_mat(j, _) = Rcpp::dnorm(NumericVector::create(beta_star_vec[j]), beta_old[j], MH_beta, true); log_prop_old_mat(j, _) = Rcpp::dnorm(NumericVector::create(beta_old[j]), beta_star_vec[j], MH_beta, true); } // Convert into vector NumericVector log_prop_star_vec = log_prop_star_mat(_, 0); NumericVector log_prop_old_vec = log_prop_old_mat(_, 0); // Calculate posterior density NumericVector log_cc_star = colSums(log_likelihood(Y, alpha_old, beta_star_vec, theta_old), true) + Rcpp::dnorm(beta_star_vec, b0, B0, true); NumericVector log_cc_old = colSums(log_likelihood(Y, alpha_old, beta_old, theta_old), true) + Rcpp::dnorm(beta_old, b0, B0, true); // Calculate acceptance probability NumericVector log_ap = pmin(log_cc_star + log_prop_old_vec - log_cc_old - log_prop_star_vec, 0); // Save samples (if log(u) < log(ap), accept alpha_star, otherwise alpha_old) NumericVector log_u = log(Rcpp::runif(J, 0, 1)); NumericVector sample(J); for (int j = 0; j < J; j++) { if (log_u[j] < log_ap[j]) { sample[j] = beta_star_vec[j]; } else { sample[j] = beta_old[j]; } } return(sample); } // SAMPLING THETA // [[Rcpp::export]] NumericVector theta_sampling(NumericMatrix Y, NumericVector alpha_old, NumericVector beta_old, NumericVector theta_old, double MH_theta, int seed) { set_seed(seed); int I = Y.nrow(); // # of theta NumericMatrix theta_star_mat(I, 1); // candidate sample for theta NumericMatrix log_prop_star_mat(I, 1); // candidate proposal density NumericMatrix log_prop_old_mat(I, 1); // old proposal density // Sampling from random-walk normal prior for (int i = 0; i < I; i++) { theta_star_mat(i, _) = Rcpp::rnorm(1, theta_old[i], MH_theta); } // Convert into vector NumericVector theta_star_vec = theta_star_mat(_, 0); // Generate proposal density for (int i = 0; i < I; i++) { log_prop_star_mat(i, _) = Rcpp::dnorm(NumericVector::create(theta_star_vec[i]), theta_old[i], MH_theta, true); log_prop_old_mat(i, _) = Rcpp::dnorm(NumericVector::create(theta_old[i]), theta_star_vec[i], MH_theta, true); } // Convert into vector NumericVector log_prop_star_vec = log_prop_star_mat(_, 0); NumericVector log_prop_old_vec = log_prop_old_mat(_, 0); // Calculate posterior density NumericVector log_cc_star = rowSums(log_likelihood(Y, alpha_old, beta_old, theta_star_vec), true) + Rcpp::dnorm(theta_star_vec, 0, 1, true); NumericVector log_cc_old = rowSums(log_likelihood(Y, alpha_old, beta_old, theta_old), true) + Rcpp::dnorm(theta_old, 0, 1, true); // Calculate acceptance probability NumericVector log_ap = pmin(log_cc_star + log_prop_old_vec - log_cc_old - log_prop_star_vec, 0); // Save samples (if log(u) < log(ap), accept alpha_star, otherwise alpha_old) NumericVector log_u = log(Rcpp::runif(I, 0, 1)); NumericVector sample(I); for (int i = 0; i < I; i++) { if (log_u[i] < log_ap[i]) { sample[i] = theta_star_vec[i]; } else { sample[i] = theta_old[i]; } } return(sample); } // MH SAMPLER // [[Rcpp::export]] List mh_sampler_irt(NumericMatrix datamatrix, NumericVector alpha, NumericVector beta, NumericVector theta, double a0 = 0, double A0 = 10, double b0 = 0, double B0 = 10, double MH_alpha = 0.3, double MH_beta = 0.3, double MH_theta = 0.3, int iter = 1000, int warmup = 1000, int thin = 1, int refresh = 100, int seed = 1) { int total_iter = iter + warmup; // total iteration int sample_iter = iter / thin; // # of samples to save NumericMatrix Y = datamatrix; // rename datamatrix to Y int I = Y.nrow(); // # of individuals int J = Y.ncol(); // # of items // rename NumericVector alpha_old = alpha; NumericVector beta_old = beta; NumericVector theta_old = theta; // create storages for parameters NumericMatrix theta_store(I, sample_iter); NumericMatrix alpha_store(J, sample_iter); NumericMatrix beta_store(J, sample_iter); // WARMUP Rcout << "Warmup: " << 1 << " / " << total_iter << " [ " << 0 << "% ]\n"; for (int g = 0; g < warmup; g++) { int seed2 = seed + g; if ((g + 1) % refresh == 0) { double gg = g + 1; double ti2 = total_iter; double per = std::round((gg / ti2) * 100); Rcout << "Warmup: " << (g + 1) << " / " << total_iter << " [ " << per << "% ]\n"; } theta = theta_sampling(Y, alpha_old, beta_old, theta_old, MH_theta, seed2); theta_old = theta; alpha = alpha_sampling(Y, alpha_old, beta_old, theta_old, a0, A0, MH_alpha, seed2); alpha_old = alpha; beta = beta_sampling(Y, alpha_old, beta_old, theta_old, b0, B0, MH_beta, seed2); beta_old = beta; } // SAMPLING double gg = warmup + 1; double ti2 = total_iter; double per = std::round((gg / ti2) * 100); Rcout << "Sampling: " << gg << " / " << total_iter << " [ " << per << "% ]\n"; for (int g = warmup; g < total_iter; g++) { int seed2 = seed + g; if ((g + 1) % refresh == 0) { double gg = g + 1; double ti2 = total_iter; double per = std::round((gg / ti2) * 100); Rcout << "Sampling: " << (g + 1) << " / " << total_iter << " [ " << per << " % ]\n"; } theta = theta_sampling(Y, alpha_old, beta_old, theta_old, MH_theta, seed2); theta_old = theta; alpha = alpha_sampling(Y, alpha_old, beta_old, theta_old, a0, A0, MH_alpha, seed2); alpha_old = alpha; beta = beta_sampling(Y, alpha_old, beta_old, theta_old, b0, B0, MH_beta, seed2); beta_old = beta; if (g % thin == 0) { double th = thin; double wu = warmup; double gg = g; double ggg = (g - warmup) / thin; // Reparameterization (fix theta with mean 0 and sd 1) NumericVector theta_std = (theta_old - mean(theta_old)) / sd(theta_old); NumericVector alpha_std = beta_old * mean(theta_old) - alpha_old; NumericVector beta_std = beta_old * sd(theta_old); theta_store(_, ggg) = theta_std; alpha_store(_, ggg) = alpha_std; beta_store(_, ggg) = beta_std; } } List L = List::create(Named("alpha") = alpha_store, Named("beta") = beta_store, Named("theta") = theta_store); return(L); }
このコードの中には6つの関数があり、以下がそれぞれの説明です。
set_seed: 再現性確保のためにRからset.seed関数を持ってきて、C++内で使用できるようにする関数log_likelihood: 対数尤度を計算する関数alpha_sampling: $\alpha$をサンプリングする関数beta_sampling: $\beta$をサンプリングする関数theta_sampling: $\theta$をサンプリングする関数mh_sampler_irt: MCMCを回す関数
また、{MCMCpack}パッケージにあるMCMCirt1d関数においては、パラメータのreparameterizationが行われており、$\theta_i$の標準化 (mean = 0, sd = 1)、またそれに付随する項目パラメータの変換が為されています。今回の自作関数でも同様の計算を施します。Reparameterization後のパラメータをそれぞれ$\alpha_j^{adj}$, $\beta_j^{adj}$, $\theta_{i}^{adj}$とすると
\begin{align*}
\theta_i^{adj} &= \frac{\theta_i - \overline{\theta}}{s_\theta}\\
\beta_j^{adj} &= \beta_j s_\theta\\
\alpha_j^{adj} &= \beta_j \overline{\theta} - \alpha_j
\end{align*}
が得られます。また、先ほどの関数たちをRで使用するためにコンパイルします。
library(Rcpp) sourceCpp("cpp/irt_sampling.cpp")
さらに、事後分布の処理や諸々の利便化のために、下記のirt_cpp関数を定義します。
irt_cpp <- function(datamatrix, iter, warmup, thin, refresh, seed, init, tuning_par, prior) { cat("\n================================================================\n") cat("Run Metropolis-Hasting Sampler for 2PL item response model...\n\n") cat(" Observations:", nrow(datamatrix) * ncol(datamatrix),"\n") cat(" Number of individuals:", nrow(datamatrix), "\n") cat(" Number of items:", ncol(datamatrix), "\n") cat(" Total correct response:", sum(as.numeric(datamatrix), na.rm = TRUE), "/", nrow(datamatrix) * ncol(datamatrix), "[", round(sum(as.numeric(datamatrix), na.rm = TRUE) / (nrow(datamatrix) * ncol(datamatrix)), 2) * 100, "%]","\n\n") cat(" Priors: \n") cat(" alpha ~", paste0("N(", prior$a0, ", ", prior$A0, "),"), "beta ~", paste0("N(", prior$b0, ", ", prior$B0, "),"), "theta ~ N(0, 1).\n") cat("================================================================\n\n") # Preparation ## Measure starting time stime <- proc.time()[3] ## Set seed set.seed(seed) # Run sampler mcmc <- mh_sampler_irt(datamatrix = Y, alpha = init$alpha, beta = init$beta, theta = init$theta, a0 = prior$a0, A0 = prior$A0, b0 = prior$b0, B0 = prior$B0, MH_alpha = tuning_par$alpha, MH_beta = tuning_par$beta, MH_theta = tuning_par$theta, iter = iter, warmup = warmup, thin = thin, refresh = refresh) # Generate variable labels label_iter <- paste0("iter_", 1:(iter/thin)) alpha_lab <- paste0("alpha_", colnames(datamatrix)) beta_lab <- paste0("beta_", colnames(datamatrix)) theta_lab <- paste0("theta_", rownames(datamatrix)) colnames(mcmc$alpha) <- colnames(mcmc$beta) <- colnames(mcmc$theta) <- label_iter rownames(mcmc$alpha) <- alpha_lab rownames(mcmc$beta) <- beta_lab rownames(mcmc$theta) <- theta_lab # Redefine quantile function lwr <- function(x) quantile(x, probs = 0.025, na.rm = TRUE) upr <- function(x) quantile(x, probs = 0.975, na.rm = TRUE) mean_ <- function(x) mean(x, na.rm = TRUE) median_ <- function(x) median(x, na.rm = TRUE) # Calculate statistics alpha_post <- data.frame(parameter = alpha_lab, mean = apply(mcmc$alpha, 1, mean_), median = apply(mcmc$alpha, 1, median_), lwr = apply(mcmc$alpha, 1, lwr), upr = apply(mcmc$alpha, 1, upr)) beta_post <- data.frame(parameter = beta_lab, mean = apply(mcmc$beta, 1, mean_), median = apply(mcmc$beta, 1, median_), lwr = apply(mcmc$beta, 1, lwr), upr = apply(mcmc$beta, 1, upr)) theta_post <- data.frame(parameter = theta_lab, mean = apply(mcmc$theta, 1, mean_), median = apply(mcmc$theta, 1, median_), lwr = apply(mcmc$theta, 1, lwr), upr = apply(mcmc$theta, 1, upr)) # Aggregate result <- list(summary = list(alpha = alpha_post, beta = beta_post, theta = theta_post), sample = list(alpha = mcmc$alpha, beta = mcmc$beta, theta = mcmc$theta)) etime <- proc.time()[3] cat(crayon::yellow("Done: Total time", round(etime - stime, 1), "sec\n")) return(result) }
応用：アメリカ連邦議会上院議員の理想点推定
前回はUS Supreme Courtの判決データを使って各判事の理想点推定をしましたが、今回は{MCMCpack}パッケージに収録されているアメリカ連邦議会上院の第106議会での点呼式投票データ*1を使って各議員の理想点推定をしていこうと思います。
まずは、データを読み込みましょう。
data(Senate, package = "MCMCpack")
データを確認すると、点呼式投票 (rc~列) 以外の議員に関するデータが含まれていることがわかります。
Senate[1:10, 1:10]
id statecode state party member rc1 rc2 rc3 rc4 rc5 SESSIONS 49700 41 ALABAMA 1 SESSIONS 1 0 0 0 1 SHELBY 94659 41 ALABAMA 1 SHELBY 1 0 0 0 1 MURKOWSKI 14907 81 ALASKA 1 MURKOWSKI 1 0 0 0 1 STEVENS 12109 81 ALASKA 1 STEVENS 1 0 0 0 1 KYL 15429 61 ARIZONA 1 KYL 1 0 0 0 1 MCCAIN 15039 61 ARIZONA 1 MCCAIN 1 0 0 0 1 HUTCHINSON 29306 42 ARKANSA 1 HUTCHINSON 1 0 0 0 1 LINCOLN 29305 42 ARKANSA 0 LINCOLN 1 0 0 1 0 BOXER 15011 71 CALIFOR 0 BOXER 1 1 1 1 0 FEINSTEIN 49300 71 CALIFOR 0 FEINSTEIN 1 1 1 1 0
行名は議員の名前になっており、属する州や政党のインデックス (1 = Republican, 0 = Democrat)が入ってます。とはいえ、IRTをするうえでは不要なので全部切り捨てて、またSenateデータはデータフレームなので行列に変換します。
Y <- as.matrix(Senate[, 6:length(Senate)])
次に、初期値や事前分布、MH法のチューニングパラメータを定義します。今回も項目パラメータの事前分布を$\alpha_j\sim N(0, 10)$, $\beta_j\sim N(1, 0.2)$と定めます。
# Initial values init <- list(alpha = alpha_init, beta = beta_init, theta = theta_init) # Tuning parameters tuning_par <- list(alpha = 2, beta = 0.6, theta = 0.9) # Priors prior <- list(a0 = 0, # mean for alpha A0 = 10, # sd for alpha b0 = 1, # mean for beta B0 = 0.2) # sd for beta
これらを全て先ほど定義したirt_cpp関数にぶち込んで、IRTをしていきます。
fit <- irt_cpp(datamatrix = Y, iter = 5000, warmup = 3000, thin = 5, refresh = 500, seed = 1, init = init, tuning_par = tuning_par, prior = prior)
================================================================ Run Metropolis-Hasting Sampler for 2PL item response model... Observations: 68544 Number of individuals: 102 Number of items: 672 Total correct response: 42458 / 68544 [ 62 %] Priors: alpha ~ N(0, 10), beta ~ N(1, 0.2), theta ~ N(0, 1). ================================================================ Warmup: 1 / 8000 [ 0% ] Warmup: 500 / 8000 [ 6% ] Warmup: 1000 / 8000 [ 13% ] Warmup: 1500 / 8000 [ 19% ] Warmup: 2000 / 8000 [ 25% ] Warmup: 2500 / 8000 [ 31% ] Warmup: 3000 / 8000 [ 38% ] Sampling: 3001 / 8000 [ 38% ] Sampling: 3500 / 8000 [ 44 % ] Sampling: 4000 / 8000 [ 50 % ] Sampling: 4500 / 8000 [ 56 % ] Sampling: 5000 / 8000 [ 63 % ] Sampling: 5500 / 8000 [ 69 % ] Sampling: 6000 / 8000 [ 75 % ] Sampling: 6500 / 8000 [ 81 % ] Sampling: 7000 / 8000 [ 88 % ] Sampling: 7500 / 8000 [ 94 % ] Sampling: 8000 / 8000 [ 100 % ] Done: Total time 369.8 sec
６分ほどで終わりましたが、高速化したのにもかかわらずまだちょっと遅めですね。後程、MCMCpack::MCMCirt1d関数とスピードを比べてみようと思います。ちなみに前回のRベースの関数での実行時間は約8分50秒でした。２分半強ぐらい高速化できました。
結果を確認しましょう。結果はfit$summary$で抜け出せるようにしています。
library(tidyverse) # Extract the result theta <- fit$summary$theta %>% mutate(name = rownames(Y), # Add senators' name party = Senate$party) # Add senators' party # Plot theta %>% ggplot(aes(y = reorder(name, mean), x = mean, color = factor(party))) + geom_pointrange(aes(xmin = lwr, xmax = upr)) + theme_light() + xlab("Estimated Ideal Point") + ylab("") + ggtitle("106th US Senate Roll-Call Vote") + scale_color_discrete(limits = c("1", "0"), label = c("Republican", "Democrat")) + theme(legend.position = "bottom", legend.direction = "horizontal", legend.title = element_blank())
所属党ごとの政策選好の違いがくっきりと見て取れます。ちゃんと推定できているか不安なので、MCMCpack::MCMCirt1d関数で確認してみます。また、せっかくなので実行時間も測ってみます。
# Measure running time stime <- proc.time()[3] # Run mcmc_fit <- MCMCpack::MCMCirt1d(Y, burnin = 3000, mcmc = 5000, thin = 5, seed = 1, verbose = 1000) etime <- proc.time()[3] # Output running time etime - stime
elapsed
55.228
いや、速すぎね？？？？？？こっちもCpp使ってるんですけど？？？？？まずは結果を確認します。
theta_mcmcp <- summary(mcmc_fit)$quantiles %>% as_tibble() %>% mutate(name = rownames(Y), party = Senate$party, mean = `50%`, lwr = `2.5%`, upr = `97.5%`) %>% select(name:upr) theta_mcmcp %>% ggplot(aes(y = reorder(name, mean), x = mean, color = factor(party))) + geom_pointrange(aes(xmin = lwr, xmax = upr)) + theme_light() + xlab("Estimated Ideal Point") + ylab("") + ggtitle("106th US Senate Roll-Call Vote: Estimated with MCMCpack") + scale_color_discrete(limits = c("1", "0"), label = c("Republican", "Democrat")) + theme(legend.position = "bottom", legend.direction = "horizontal", legend.title = element_blank(), axis.text.y = element_text(size = 4))
符号が自作関数と反転してますが、大体同じような結果ではないでしょうか。符号反転の問題は後からどうにでもできることなのでどうでもいいです。いやー、しかし{MCMCpack}速すぎ問題ありますよね。どうすればあそこまで速くなるんだろう…
おわりに
今回やってわかったことはMCMCpack速すぎじゃね～～～～～？？？？？？？？？ということでした。RcppArmadilloなどのLinear Algebra系のやつも使いこなせるようになればもっと速くできるのでしょうか。精進していこうと思います。
