こんにちは。データサイエンスチームの坂本と申します。
使い慣れたRを使って、Google Colabのクラウド 環境上でベイズ 推定ができたら便利ですよね。しかしやってみると意外に環境設定手順が複雑で、悩むことになるかもしれません。
TechBlog第15回では、統計解析環境Rのユーザーが、Google Colab上でベイズ 推定を行う際の手続きを紹介します。
ベイズ 推定にはハミルトニアン ・モンテカルロ法 のNUTSアルゴリズム を使用してパラメータ推定を行うためのRパッケージ、cmdstanrを利用します。
また、せっかくですので、後半では構築した環境を使った分析例もご紹介します。ウェアラブル 端末で測定した睡眠データで「睡眠リズムの安定性(不安定性)」を定量 化してみましょう。睡眠リズムを定量 化する方法の意外な難しさと、時刻データの解析に関するちょっと面白いお話も披露できたらと思います。
※2025年3月時点で適正動作が確認できた方法を紹介しています。
Google Colabでcmdstanrのインストール
まずは環境設定です。Google Colabに、統計解析環境Rのパッケージである”cmdstanr”をインストールし、動作確認を行います。
なお、作業前に、Google Drive のマイドライブ直下に、「library」という空のフォルダを作成しておきましょう。
1. [ランタイム: python ] Google Drive をマウント
from google.colab import drive
drive.mount('/content/drive/' )
2. [ランタイム: R] ライブラリ(パッケージ)のインストール先を設定
.libPaths ( "/content/drive/MyDrive/library" )
3. [ランタイム: R] cmdstanrパッケージのインストール
install.packages ( "cmdstanr" , repos = c ( "https://mc-stan.org/r-packages/" , getOption ( "repos" )))
4. [ランタイム: R] ライブラリの読み込みとcmdstanrの動作環境整備
少し時間がかかりますが、粘り強く待ちましょう。
library ( cmdstanr)
cmdstanr:: install_cmdstan ( dir = "/content/drive/MyDrive/library" , overwrite = TRUE )
5. [ランタイム: R] cmdstanrの保存パス確認とパス設定
保存されたcmdstanrのパスとバージョンの確認
cmdstanr:: cmdstan_path ()
使用するcmdstanrのパスとバージョン設定
set_cmdstan_path ( path = "/content/drive/MyDrive/library/cmdstan-2.36.0" )
6. [ランタイム:R] 動作確認: 統計モデルの作成(.stanファイル)
ここでは、正規分布 の平均と標準偏差 を推定するサンプルモデルを例示
ModelScript = "
data {
int N;
vector[N] Y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
Y ~ normal(mu, sigma);
}
"
writeLines ( ModelScript, "作業ディレクトリのパス/normal.stan" )
サンプルデータの準備とモデルのコンパイル
datastan = list (
N = 50 ,
Y = rnorm ( 50 , 0 , 1 )
)
cmd_model = cmdstan_model ( "作業ディレクトリのパス/normal.stan" )
cmd_model
7. [ランタイム: R] パラメータ推定と結果の確認(収束診断と推定値)
fit <- cmd_model$ sample(
data = datastan,
seed = 123 ,
chains = 4 ,
parallel_chains = 4 ,
iter_sampling = 1000 ,
iter_warmup = 500
)
all ( fit$ summary()[, 'rhat' ] < 1.05 )
fit$ summary()
ここまで、問題なく実行できたでしょうか。STEP4のcmdstanrの動作環境整備では、g++のインストールに時間がかかったのではないかと思います。
現在のGoogle Colabでは、cmdstanrを利用するたびに毎回g++のインストールを実行する必要があるため、不便を感じるところかもしれません。
循環正規分布 (von Mises distribution)を用いた睡眠中央時刻の平均と分散の推定
さて、cmdstanrが使えるようになったところで、実際に推定を行ってみます。
ベイズ 推定法の利用場面に相応しい、「睡眠リズムの安定性(不安定性)」の定量 化を行ってみます。
「睡眠リズムがバラバラ」という状態は、体験的にはイメージしやすいかもしれません。しかし、いざデータ解析で「バラバラ度合い」を定量 化しようとすると、悩ましい問題が出てくることに気がつきます。
一つの問題は、「睡眠時間」「入眠時刻」「起床時刻」、あるいは「1日の睡眠の(分断)回数」などのうち、どの情報のバラバラ度合いを計算すれば良いのかという問題です。
もう一つは、23時と0時の就寝では時刻の違いは1時間ですが、数値の上では0(最小値)と23(最大値)で大きく乖離してしまう「循環データ」の問題です。0時と23時のバラバラ度合いは実際には1時間ですが、数値上は23時間のようになってしまうというものです。
これら問題にはいくつかの対処方法が考えられますが、本記事では「睡眠中央時刻のバラつき度合い」を、循環データに対応した「循環正規分布 (von Mises distribution)を用いて推定する」という方法をご紹介します。
睡眠中央時刻のバラつき
この記事では、「入眠した時刻から覚醒した時刻の中間の時刻」を睡眠中央時刻と呼びます。22時に就寝して6時に起床した場合は2時、0時に就寝して7時に起床した場合は3時半です。
規則正しい生活を送っている場合、睡眠中央時刻は特定の時刻付近に集中します。反対に、睡眠リズムが不規則な場合には、様々な時刻に睡眠中央時刻が現れます。1周が24時間の「24時間時計」を用いて、2名の成人男性の1年間の睡眠中央時刻の分布を見てみましょう(図1 )。
※いずれも測定端末はGoogle Fitbit(ID=1はCharge5、ID=2はCharge6を装着)。
図1 24時間時計(円周)上の睡眠中央時刻の分布(昼寝等を含む約1年間の睡眠データより)
循環正規分布 を用いたパラメータ推定
図1に示した2名の睡眠中央時刻データを利用し、それぞれの睡眠中央時刻の円周上の平均と分散パラメータを推定します。von Mises distributionとstanモデルの導入解説は、下記のブログ記事に大変わかりやすく整理されていますので、ご参照ください。
bayesmax.sblo.jp
睡眠中央時刻は00:00:00が0度、12:00:00が180度となるようあらかじめ時刻を角度に変換し、角度に対してπ/ 180を乗算することでラジアン も計算しておきました。stanモデルは次の通りです。なお、データブロックに記載されたradian_lowerとradian_upperは、実際の睡眠中央時刻から計算された観測データ(ラジアン )の25%タイル点と75%タイル点です。
ModelScript = "
//vonMises_difference.stan
data {
int n1;
int n2;
vector[n1] radian1;
vector[n2] radian2;
real radian1_lower;
real radian1_upper;
real radian2_lower;
real radian2_upper;
}
parameters {
real<lower=radian1_lower, upper=radian1_upper> mu1;
real<lower=radian2_lower, upper=radian2_upper> mu2;
vector<lower=0, upper=200>[2] kappa;
}
transformed parameters {
}
model {
for (n in 1:n1){
radian1[n] ~ von_mises(mu1, kappa[1]); //群1のモデル
}
for (m in 1:n2){
radian2[m] ~ von_mises(mu2, kappa[2]); //群2のモデル
}
}
generated quantities {
real mu1_angle; //mu1(ラジアン)を度数に変換
real mu2_angle; //mu2(ラジアン)を度数に変換
vector[2] a; //平均合成ベクトル長(円周分散・円周標準偏差の計算に使用)
vector[2] v; //群1と群2の円周分散
vector[2] nu; //群1と群2の円周標準偏差
vector[2] nu_angle; //nu(ラジアン)を度数に変換
real mu_angle_diff; //mu1_angleとmu2_angleの差(最短距離)
real diff_nu; //nu1とnu2の差
real diff_mu_angle; //nu1とnu2の差
// feature calculation
for (k in 1:2){
a[k] = modified_bessel_first_kind(1,kappa[k])/modified_bessel_first_kind(0,kappa[k]);
v[k] = 1 - a[k];
nu[k] = sqrt(-2 * log(a[k]));
nu_angle[k] = nu[k] / pi() * 180;
}
mu1_angle = mu1 / pi() * 180;
mu2_angle = mu2 / pi() * 180;
mu_angle_diff = mu2_angle - mu1_angle;
diff_nu = nu[2] - nu[1];
diff_mu_angle = mu1_angle - mu2_angle;
}
"
推定された結果を確認してみましょう。chains = 4、iter _sampling = 2000、iter _warmup = 1000、max_treedepth = 10、thin = 1、adapt_delta = 0.8の条件で、全てのパラメータの収束(Rhat統計量 < 1.05)を確認しています(図2 )。
図2 循環正規分布 (von Mises distribution)を用いた睡眠中央時刻の円周上平均と分散パラメータの推定値(ヒストグラム のエラーバーは±1SDを示し、棒グラフのエラーバーは分散パラメータの95%信用区間 を示す) パラメータ推定の結果、睡眠が不規則なID=1の対象者では、睡眠中央時刻の円周上平均(μ)が05:11ごろとなり、バラつき指標(ν)は5.01(単位: hour)となりました。また、規則的な生活を送っているID=2の対象者は、睡眠中央時刻の円周上平均(μ)が02:04ごろとなり、バラつき指標(ν)が0.67(単位:hour)となりました。ID=1の対象者の方が寝ている時間が夜間遅い方向にあり、寝ている時間帯がID=2よりもバラついている、つまり「睡眠リズムがバラバラ」ということを示す結果となっています。
おわりに
この記事では、「解析環境の設定」と「時刻データの統計学 的処理方法」を紹介しました。データ解析者にとっては悩ましい二つのポイントへの処方箋となれば嬉しく思います。
また、長期間にわたって記録された睡眠データの可視化と解析結果に、面白さを感じていただけたら幸いです。
書いた人:坂本