はじめに
こんにちは、とある大学生です。
みなさんは睡眠をしっかりとっていますでしょうか？
今回は大学生の（自分の）就寝時間、起床時間、そして睡眠時間のデータをMarkov Chain Monte Carloを通して解析していきたいと思います。
ちなみにこれは元ネタがありまして、以下の記事に掲載されています。
また、コードはipynb形式でGitHubにあげております。汚いですが、よかったらどうぞ。データもすべて整形して準備OKなものをアップロードしているので、JupyterでRunするだけでOKです。
データについて：就寝データ、起床データ、睡眠時間データ
今回利用したデータはiPhoneのアプリ「健康リズムケア」に手動で過去二ヶ月（2019年7月から8月）入力していた自分の睡眠データです。
就寝データと起床データ、そして睡眠時間データの3データセットあります。また、大学生の理想な就寝時間は12:00、起床時間は8:00（根拠はないです）という仮定にそって話します。
就寝データ
以下が就寝データで、Indicator, time_offset, そしてtime_rangeからなります。
| Indicator | time_offset | time_range | |
|---|---|---|---|
| 0 | 0 | -60.0 | 2019-06-30 23:00:00 |
| 1 | 0 | -59.0 | 2019-06-30 23:01:00 |
| 2 | 0 | -58.0 | 2019-06-30 23:02:00 |
| 3 | 0 | -57.0 | 2019-06-30 23:03:00 |
| 4 | 0 | -56.0 | 2019-06-30 23:04:00 |
time_rangeカラムは時刻を表し、time_offsetは理想の就寝時間と時刻の差分を表します。
Indicatorは起きている間は0、寝ている間は1を表します。以下だと、2019年7月1日は12:05に就寝したので、そこから0から1に切り替わっており、そのときのtime_offsetは12:00との差分をとって5分となります。
| Indicator | time_offset | time_range | |
|---|---|---|---|
| 63 | 0 | 3.0 | 2019-07-01 00:03:00 |
| 64 | 0 | 4.0 | 2019-07-01 00:04:00 |
| 65 | 1 | 5.0 | 2019-07-01 00:05:00 |
| 66 | 1 | 6.0 | 2019-07-01 00:06:00 |
| 67 | 1 | 7.0 | 2019-07-01 00:07:00 |
起床データ
起床データも構造的に同じです。起床の理想時間を8:00に設定しているので、7:05に起きるとtime_offsetは-55分、indicatorの切り替えは7:05で行われます。
| Indicator | time_offset | time_range | |
|---|---|---|---|
| 62 | 1 | -58.0 | 2019-07-01 07:02:00 |
| 63 | 1 | -57.0 | 2019-07-01 07:03:00 |
| 64 | 1 | -56.0 | 2019-07-01 07:04:00 |
| 65 | 0 | -55.0 | 2019-07-01 07:05:00 |
| 66 | 0 | -54.0 | 2019-07-01 07:06:00 |
睡眠時間データ
睡眠時間データは構造的にすこし異なります。
| Date | Sleep | Wake | length | |
|---|---|---|---|---|
| 0 | 2019/7/1 | -5.0 | 55.0 | 7.000000 |
| 1 | 2019/7/2 | 30.0 | 20.0 | 8.166667 |
| 2 | 2019/7/3 | -95.0 | 15.0 | 6.166667 |
| 3 | 2019/7/4 | -2.0 | 10.0 | 7.800000 |
Sleep Offset と Wake Offset が表示されると同時に、lengthカラムに睡眠時間も１時間単位で表示されます。
解析開始！まずは可視化する
それではこのデータを使って分析していきます。
Scatterplotを使って就寝データを可視化
データを見ると、深夜0:00あたりからAsleepの軸が濃い青色に変化し始めていき、Awakeもどんどん淡い青色に変化していきます。逆に端っこの23:00ではフルAwakeで、午前3:30はさすがに寝ちゃっています。ちなみにグラフには表示していませんが平均的な就寝時間は12:50〜1:00です。大学生です。
Scatterplotを使って起床データを可視化
逆に起床時間はどうでしょう。6:00はほぼまだAsleepでAwakeになる気配はまったくないですが、8:00になると濃くなりはじめます。そういえば今年の夏休みは長期インターンをしていたので比較的健康的です（起床時間の平均は7月は8:20、8月は9:00ぐらいです）。
モデルを組む：シグモイド関数のパラメータを推定する
ではどのようなモデルを組むのがよいでしょうか。夜の12:00に寝ているかどうかはファジーなので夜の12:00に寝ている確率、と言い換えたほうが良さそうです。また、睡眠の状態から起きている状態への切り替えを表すモデルにしたいので、シグモイド関数を使っていくといいでしょう。
シグモイド関数
シグモイド関数はは、$ \frac{1}{1+exp^(α + βt)} $ のようにあらわせます。ここではこのαとβをMarkov Chain Monte Carloを使って推定していくことになります。異なるα,βを選ぶことで異なるsigmoid関数の分布が得られるようになります。
##マルコフチェイン
ざっくりいいますとマルコフチェインは以前の状態をもとに現在の状態を当てる方法（または現在の状態をもとに次の状態をあてる方法）です。よくあるのが天候で、高校数学の美しい物語さんがうまくそれを説明されています。
###モンテカルロ
ざっくりいいますとモンテカルロはランダムなサンプルを大量に生産し、そこから確率分布を推定する方法です。おそらくわかりやすいのが、以下のツイートです。大量の滋賀県上にランダムな座標のサンプルを生産したら、1/6の確率で琵琶湖のなかに落ちるという例です。
###マルコフチェインモンテカルロ
以上の方法を組み合わせたのがマルコフチェインモンテカルロです。分布からランダムなサンプルを大量に抽出するわけですが（モンテカルロ）、そのサンプルをとる際には現在の状態に依存します（マルコフチェイン）。難しいので今回の例にあてはめると
シグモイド関数のαとβを推定したい：
- α,βの初期値を用意
- 現在の状態から新たなαとβをランダムに用意
- データと照らし合わせ、適用可能だったら更新、不可能だったら拒否
- 2と3を繰りかえす。
3では今回はMetropolis Hastingsというアルゴリズムを使って判断します。
マルコフチェインモンテカルロを実装し、就寝データに適用
コードはほぼほぼ上記のMarkov Chain Monte Carlo in Pythonからとっています。これを就寝データ（sleep_obs）に適用します。
with pm.Model() as sleep_model:
# alphaとbetaを正規分布を事前分布として推定
alpha = pm.Normal('alpha', mu=0.0, tau=0.01, testval=0.0)
beta = pm.Normal('beta', mu=0.0, tau=0.01, testval=0.0)
# alphaとbetaをもとにシグモイド関数を推定
p = pm.Deterministic('p', 1. / (1. + tt.exp(beta * time + alpha)))
# 観測値を変換
observed = pm.Bernoulli('obs', p, observed=sleep_obs)
# Metropolis Hastings アルゴリズム
step = pm.Metropolis()
# サンプルを取る
sleep_trace = pm.sample(N_SAMPLES, step=step)
ここでαとβのトレースをみてみましょう。
図を見ると、両方ともめちゃめちゃ振動していてまったく収束する気配が見えません。収束するまで時間がかかることと初期値はだいたいまちがっていることから最初の90%の結果はあてにしないことが多いようです。ちなみに収束しているかどうかの判断も難しいそうで、割愛します。
最終的にはα = 1.35, β =-0.0270あたりになっていることがわかります。
就寝時間の確率分布を可視化
α,βを推定できたので、ここから就寝時間の分布がどうなっているかをみていきます。過去の5000サンプルを使って可視化していきましょう。
いい感じです。グラフだけみてると、0:45ぐらいに50%をこえるのがわかります。つまり2日に1回は0:45以降になってようやく寝る、みたいな現象が起きているわけです。また2:30になっても90%ぐらいの確率でまだ起きていることがあります。大学生。。。
ちなみにある時刻で就寝しているかどうかの確率も割り出すことができます。
print(f'11:30 PM probability of being asleep: {round(100 * logistic(-30, beta_est, alpha_est))} %')
print(f'00:00 PM probability of being asleep: {round(100 * logistic(0, beta_est, alpha_est))} %')
print(f'00:30 PM probability of being asleep: {round(100 * logistic(30, beta_est, alpha_est))} %')
print(f'00:50 PM probability of being asleep: {round(100 * logistic(50, beta_est, alpha_est))} %')
print(f'01:30 PM probability of being asleep: {round(100 * logistic(90, beta_est, alpha_est))} %')
#Output:
#11:30 PM probability of being asleep: 10.0 %
#00:00 PM probability of being asleep: 20.0 %
#00:30 PM probability of being asleep: 37.0 %
#00:50 PM probability of being asleep: 50.0 %
#01:30 PM probability of being asleep: 75.0 %
| Time | 寝てる確率 | |
|---|---|---|
| 0 | 23:30 | 10 % |
| 1 | 00:00 | 20 % |
| 2 | 00:30 | 37 % |
| 3 | 00:50 | 50 % |
| 4 | 01:30 | 75 % |
たしかにこの時期の平均就寝時間は0:50だったので、これは妥当だといえます。ていうか深夜1:00ぐらいでも寝てる確率50%で1:30でも75%なのか。。。
起床時間の分布を可視化
起床時間の分布も就寝時間の分布と同様計算可能です。
結果は上記のようになりました。8:00には起床する確率が25%あたりになっており、10:00にはほぼほぼ起きていますが10%ぐらいの確率でまだ寝ています。とりあえず遅寝遅起きが目立ちますね。。。
（今気づいたのですが、インターンは8月中旬で終わり、そのあとが結構堕落したからかもしれないです。）その分の堕落が関数を右におしたのかもしれません。
定量的にはかっていきましょう。
print(f"Probability of being awake at 7:30 AM: {round(100 - (100 * logistic(-30, beta=beta_est, alpha=alpha_est)))} %")
print(f"Probability of being awake at 8:00 AM: {round(100 - (100 * logistic(0, beta=beta_est, alpha=alpha_est)))} %")
print(f"Probability of being awake at 8:30 AM: {round(100 - (100 * logistic(30, beta=beta_est, alpha=alpha_est)))} %")
print(f"Probability of being awake at 9:00 AM: {round(100 - (100 * logistic(60, beta=beta_est, alpha=alpha_est)))} %")
print(f"Probability of being awake at 9:30 AM: {round(100 - (100 * logistic(90, beta=beta_est, alpha=alpha_est)))} %")
#Output
#Probability of being awake at 7:30 AM: 17.0 %
#Probability of being awake at 8:00 AM: 29.0 %
#Probability of being awake at 8:30 AM: 46.0 %
#Probability of being awake at 9:00 AM: 63.0 %
#Probability of being awake at 9:30 AM: 78.0 %
| Time | 起きてる確率 | |
|---|---|---|
| 0 | 7:30 | 17 % |
| 1 | 8:00 | 29 % |
| 2 | 8:30 | 46 % |
| 3 | 9:00 | 63 % |
| 4 | 9:30 | 78 % |
9:30でも4日に1回は9:30でもまだ寝てるのか。。。
睡眠時間の分布
最後に睡眠時間の分布です。
最頻値は7と8時間の間にきています。平均睡眠時間が7時間50分ぐらいなので、妥当です。3時間が1回しか観測されていないのはホワイトです。よくみると（あたりまえかも）早起きより遅起きのほうが多いので、right skewな確率密度関数でフィットしたほうがよさそうです。
これで寝溜めしている現象もあらわしやすくなります。
モデルにフィット
再び正規分布を事前分布として計算していきます。
結果だけお見せします。
最頻値はたしかに7時間半ぐらいに来ており、平均に集まってきています。事後分布もたしかにright skewedになっています。
duration_lis = [6.5, 8.0, 9.0]
for duration in duration_lis:
prob = round(100 * (1 - stats.skewnorm.cdf(duration, a = alpha_skew_est,
loc = mu_est, scale = 1 / tau_est)))
print(f"Probability of at least {duration} hours of sleep: "
f"{prob}")
#Output
#Probability of at least 6.5 hours of sleep: 80.0
#Probability of at least 8.0 hours of sleep: 41.0
#Probability of at least 9.0 hours of sleep: 19.0
| Duration | 実現できる確率 | |
|---|---|---|
| 0 | 6.5 時間 | 80 % |
| 1 | 8.0 時間 | 41 % |
| 2 | 9.0 時間 | 19 % |
5日に2日は必ず8時間寝ているのは健康的ですね。6.5時間も80%で実現できているのは素晴らしいです。
## 最後に
今回は正規分布を事前分布にとっていたわけですが、ポアソンとかほかの分布を使うことができたら結果も変わって面白そうです。まだまだ統計なんもわかっていないなあと未熟さも感じました。
あと感想なのですが、モデルのコードはいろんなところに落ちているのでいいとして、アプリのデータを整形するのが一番大変でした。。。