Google ColabでRによるベイズ推定を行う方法: 睡眠リズムの定量化を例に

こんにちは。データサイエンスチームの坂本と申します。

使い慣れた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をマウント

#@title Python3: Google Driveのマウント


# content/driveをマウント
from google.colab import drive
drive.mount('/content/drive/')

2. [ランタイム: R] ライブラリ(パッケージ)のインストール先を設定

# Rのライブラリを管理するフォルダのパスを追加 
.libPaths("/content/drive/MyDrive/library")

3. [ランタイム: R] cmdstanrパッケージのインストール

#@title Installation "cmdstanr"

# 最初の一回だけ実行
install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

4. [ランタイム: R] ライブラリの読み込みとcmdstanrの動作環境整備
少し時間がかかりますが、粘り強く待ちましょう。

# cmdstanr
library(cmdstanr)
# 使う時に実行(接続した環境にinstallする必要がある)
cmdstanr::install_cmdstan(dir = "/content/drive/MyDrive/library", overwrite = TRUE)

5. [ランタイム: R] cmdstanrの保存パス確認とパス設定
保存されたcmdstanrのパスとバージョンの確認

#@title Unpacked Directory
cmdstanr::cmdstan_path()

使用するcmdstanrのパスとバージョン設定

#@title Path Setting for CmdStanR
set_cmdstan_path(path = "/content/drive/MyDrive/library/cmdstan-2.36.0")

6. [ランタイム:R] 動作確認: 統計モデルの作成(.stanファイル)
ここでは、正規分布の平均と標準偏差を推定するサンプルモデルを例示

#@title Stan Model and Save as .stan file

ModelScript = "
 data {
   int N;
   vector[N] Y;
 }
 parameters {
   real mu;
   real<lower=0> sigma;
 }
 model {
   Y ~ normal(mu, sigma);
 }

"

# save as .stan file
writeLines(ModelScript, "作業ディレクトリのパス/normal.stan")

サンプルデータの準備とモデルのコンパイル

#@title Sample Data & Compile
datastan = list(
 N = 50,
 Y = rnorm(50, 0, 1)
)

# model compilation
cmd_model = cmdstan_model("作業ディレクトリのパス/normal.stan")

# model script
cmd_model

7. [ランタイム: R] パラメータ推定と結果の確認(収束診断と推定値)

#@title Sampling(Parameter Estimation) and Results


# パラメータ推定
fit <- cmd_model$sample(
 data = datastan,
 seed = 123,
 chains = 4,           # チェーン数
 parallel_chains = 4,  # 並列処理
 iter_sampling = 1000, # サンプリング回数
 iter_warmup = 500     # ウォームアップ期間
)

#収束診断 (Stan Recommendation Rhats < 1.05)
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%タイル点です。

#@title von-Mises Difference.stan

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)。

パラメータ推定の結果、睡眠が不規則なID=1の対象者では、睡眠中央時刻の円周上平均(μ)が05:11ごろとなり、バラつき指標(ν)は5.01(単位: hour)となりました。また、規則的な生活を送っているID=2の対象者は、睡眠中央時刻の円周上平均(μ)が02:04ごろとなり、バラつき指標(ν)が0.67(単位：hour)となりました。ID=1の対象者の方が寝ている時間が夜間遅い方向にあり、寝ている時間帯がID=2よりもバラついている、つまり「睡眠リズムがバラバラ」ということを示す結果となっています。