2022-09-29

【Rでの時刻データの扱い】as.POSIXctを使うの卒業します

R言語

1時間ほど前に書いた記事では，Rで文字列を時刻データに変換するのに「私はas.POSIXctを使います」と言ってましたが，訂正します．

　私はas.POSIXctを２度と使いません．

　今なら断言できます．Rでas.POSIXctを使ってる奴，素人です．

　as.POSIXctを使ってる人，人生を無駄にしています．

　as.POSIXctを使ってる勘違い野郎，さっきまでの私です．

　では，代わりに何を使うか．それが，fasttimeパッケージのfastPOSIXctです．さっき，偶然，インターネット検索で見つけましたが，このコマンドが私の人生を救ってくれると思います．それくらい，時系列解析，特に実世界データを分析する際には，時刻データの処理の高速化が重要です．なら，C言語使えよと言われるかもしれませんが，C言語で日付や時刻を扱うの面倒くさいです．１０年前は，C言語とRの使用時間の比は9:1くらいでしたが，今は，1:99くらいです．

fasttimeパッケージのインストール
処理時間の比較
【重要】日本時間での使い方
まとめ

`fasttime`パッケージのインストール

　一度もfasttimeパッケージを使ったことがない人は，

install.packages("fasttime")

を実行してパッケージをインストールしてください．

処理時間の比較

　fastPOSIXctのすごさを体感するために，数値実験をしてみます．

　使うのは以下のRスクリプトです．

# 日時の文字列データを生成 (テスト用に時刻データを文字列に変換)
time.chr <- as.character(seq(as.POSIXct("2022-09-20 00:00:00"),as.POSIXct("2022-09-30 00:00:00"),"1 sec"))
###############################################
# Rの関数as.POSIXctを使った場合
system.time(time <- as.POSIXct(time.chr))
###############################################
# "fasttime"パッケージのfastPOSIXctを使った場合
# パッケージインストール用
#　install.packages("fasttime")
require(fasttime)
system.time(time <- fastPOSIXct(time.chr))

　
　このRスクリプトを実行すると，処理時間が表示されます．ここで使った，system.timeは処理時間を計算してくれるコマンドです．

　as.POSIXctを使って文字列を時刻データに変換すると，

> system.time(time <- as.POSIXct(time.chr))
   ユーザ   システム       経過  
     13.67       0.00      13.67

という結果になりました．注目するのはユーザの部分の値で，処理時間が13.67秒かかったということです．

　次に，fastPOSIXctを使って文字列を時刻データに変換すると，

> system.time(time <- fastPOSIXct(time.chr))
   ユーザ   システム       経過  
      0.09       0.00       0.09

という結果になりました．処理時間は，0.09秒でした．100m走のタイムが1/100になったくらいの衝撃です．

　fastPOSIXctを使って変換しても，結果はas.POSIXctと同じです．

【重要】日本時間での使い方

　fastPOSIXctを使うときの注意は，時刻のタイムスタンプが"GMT" (グリニッジ標準時)になってしまうことです．それと，常に1970-01-01 00:00:00が基準になり，as.POSIXctにあったoriginのオプションはありません．

　タイムスタンプが"GMT"というのは不便です．現状では，日本時間で時刻を扱うときは以下のようにする必要があります．

time <- fastPOSIXct(time, tz="Japan")-9*60*60

ここでは，タイムスタンプを"JST" (日本標準時)に指定して，"GMT"から"JST"の時差として勝手に補正される9時間の時差を差し引くことで，正しい日本時間にしています．

まとめ

　そのうち，as.POSIXctは改良されて速くなると思いますが，現状では，as.POSIXctなんて使ってられません．

　私は何ヶ月にもわたる長期計測データを分析することが多いので，fastPOSIXctはむちゃくちゃ助かります．今まで時刻の処理だけで10分くらいかかっていたものが，５秒以下で終わるって，感動です．

2022-09-29

【Rで日付，時間を扱う】as.Dateを使うときはタイムゾーンを必ず指定しろ

R言語

私はサツマイモ好きなので，コンビニのサツマイモ系スイーツは必ずチェックします．今日は大学のローソンで「さつまいもあんぱんコグマパン」という菓子パンが新発売されていたのでさっそく買って食べてみました．ということで，今回はRで時刻を日付に変えるときの注意です．

`as.Date`を使うときは`tz`を指定しろ

　時刻データを扱うとき，私はas.POSIXctを使います．例えば，

　"2022/09/29 21:37:00"

という日付と時刻からなる文字列を，時刻データに変換するときは，

as.POSIXct("2022-09-29 21:37:00")

とします．Rの時刻データは，

　"2022-09-29 21:37:00 JST"

のように表示されます．"JST"の部分はタイムゾーンを表し，日本の時刻ということです．日本語環境だとデフォルトで"JST"になると思います．

　このような時刻データから，日付だけを取り出したいとき，as.Dateを使います．例えば，

time <- as.POSIXct("2022-09-29 21:37:00")
as.Date(time, tz="Japan")

とします．今回のポイントはas.Dateでは，as.Date(... , tz="Japan")のように必ずタイムゾーンを指定しろということです．

　そのうち修正が入るかもしれませんが，as.POSIXctは気を利かせてタイムゾーンを私が住んでいる日本にしてくれますが，as.Dateは，デフォルトでタイムゾーンがグリニッジ標準時 (Greenwich Mean Time: GMT)になります．日本とは9時間の時差があります (日本は西側の国より未来にいます)．

　つまり，as.Date(time)として，タイムゾーンtzの指定を省略すると，日本とは違う日付になることがあります．

日付がずれることの確認

　結果はわかっていますが，as.Dateでタイムゾーンtzの指定を省略すると，日付がずれることの確認です．以下のスクリプトを実行すると，

# 日時データを生成
date_time <- seq(as.POSIXct("2022-09-28 00:00:00"),as.POSIXct("2022-09-30 00:00:00"),"4 hours")
date_time <- as.POSIXct(date_time,tz="GMT")
# as.Dateを使って日付を抽出
date <- as.Date(date_time)
# 曜日を抽出
weekday <- weekdays(date_time,abbreviate=TRUE)
# as.Dateの挙動を確認
data.frame(date_time,weekday,as.Date=date,date=format(date_time,"%Y-%m-%d"))

以下のようになります．

             date_time weekday    as.Date       date
1  2022-09-28 00:00:00      水 2022-09-27 2022-09-28
2  2022-09-28 04:00:00      水 2022-09-27 2022-09-28
3  2022-09-28 08:00:00      水 2022-09-27 2022-09-28
4  2022-09-28 12:00:00      水 2022-09-28 2022-09-28
5  2022-09-28 16:00:00      水 2022-09-28 2022-09-28
6  2022-09-28 20:00:00      水 2022-09-28 2022-09-28
7  2022-09-29 00:00:00      木 2022-09-28 2022-09-29
8  2022-09-29 04:00:00      木 2022-09-28 2022-09-29
9  2022-09-29 08:00:00      木 2022-09-28 2022-09-29
10 2022-09-29 12:00:00      木 2022-09-29 2022-09-29
11 2022-09-29 16:00:00      木 2022-09-29 2022-09-29
12 2022-09-29 20:00:00      木 2022-09-29 2022-09-29
13 2022-09-30 00:00:00      金 2022-09-29 2022-09-30

となります．上の出力の，as.Dateの列がas.Dateの結果です．dateの列が正しい日付なので，as.Dateの7から9行目と13行の結果は，日にちがずれてます．

まとめ

　大規模な生体データを分析するときは，日付や時間を変数として使いこなすことが必要です．これから時々，時刻データを扱うテクニックを紹介していきます．

2022-09-28

【Rで活動量分析】オムロン研究用活動量計HJA-750Cのデータファイルの読み込み

R言語

オムロン活動量計でHJA-750Cというのがあります．

www.healthcare.omron.co.jp

これは研究や臨床検査用なので，一般には発売されていません．大学などの研究機関に所属していれば購入可能です．

　といことで，今回はHJA-750Cを使って研究をしたい人向けのに，Rを使ったデータファイルの読み込み方の説明です．

活動量評価のデータファイルの中身

　HJA-750Cの専用アプリで，計測結果 (集計ファイル作成)を書き出すと，

　(被験者ID)10SECMETS_(計測日).csv

のような名前がついたファイルが作成されます．ファイル名の，(被験者ID)と(計測日)の部分は，それぞれ，各自が設定したID，計測日になります．今回はこのファイルを読み込んでみます．

　計測結果のファイルの中身は，

(被験者ID)
2022/09/28

時刻,活動強度,運動種別
00:00:00,0,計測なし
00:00:10,0,計測なし
00:00:20,0,計測なし

のような感じです．この2行目が計測日です．4行目が各列の見出しで，5行目以下に計測結果が10秒ごとに記録されています．

　HJA-750Cの計測データで注意してほしい点は，記録される時刻が必ず00:00:00からはじまることです．例えば，その日の夕方17:00:00に計測しても，データは00:00:00から記録されるので，00:00:00 ～ 17:00:00までの時間帯は，活動強度がずーっと「0」，運動種別がずーっと「計測なし」になります．

活動量ファイルの読み込み

　Rでファイルを読み込みたいとき，私であれば以下のようにします．

# データがあるフォルダの指定
DIR <- "ここにデータファイルがあるフォルダを指定"
# データファイルの指定
FN <- "10SECMETSファイルの名前を指定"
####
# フォルダの移動
setwd(DIR)
# ファイルから日付の読み込み
YMD <- read.csv(FN,header=FALSE,skip=1,stringsAsFactors=FALSE,nrows=1)[1,"V1"]
# ファイルからデータの読み込み
DAT <- read.csv(FN,header=FALSE,skip=4,stringsAsFactors=FALSE)
# DATの各列の変数名の設定
colnames(DAT) <- c("time","METs","activity_type")
# timeを時刻データに変換 (日付＋時刻)
DAT$time <- as.POSIXct(paste(YMD,DAT$time),origin="1970-1-1 00:00:00")
# グラフを描画
plot(DAT$time,DAT$METs,"l",col=2,xlab="Time",ylab="METs")

Windowsで"... invalid multibyte character ..."というエラーが出たら，上のスクリプトの11行目を

DAT <- read.csv(FN,header=FALSE,skip=4,stringsAsFactors=FALSE,fileEncoding="CP932")

に修正してください．

　上のRスクリプトを実行すると，DATの中身は以下のようになります．

time	METs	activity_type
2022-07-02 09:48:30	6.3	生活活動
2022-07-02 09:48:40	2.6	生活活動
2022-07-02 09:49:00	1.0	生活活動
2022-07-02 09:49:10	1.0	生活活動
2022-07-02 09:49:20	1.0	生活活動
2022-07-02 09:49:30	1.0	生活活動

　さらに，「計測なし」の行を除きたいときは，以下のようにします．

# データがあるフォルダの指定
DIR <- "ここにデータファイルがあるフォルダを指定"
# データファイルの指定
FN <- "10SECMETSファイルの名前を指定"
####
# フォルダの移動
setwd(DIR)
# ファイルから日付の読み込み
YMD <- read.csv(FN,header=FALSE,skip=1,stringsAsFactors=FALSE,nrows=1)[1,"V1"]
# ファイルからデータの読み込み
DAT <- read.csv(FN,header=FALSE,skip=4,stringsAsFactors=FALSE)
# 「計測なし」の行を削除
DAT <- DAT[DAT$V3!="計測なし",]
# DATの各列の変数名の設定
colnames(DAT) <- c("time","METs","activity_type")
# timeを時刻データに変換 (日付＋時刻)
DAT$time <- as.POSIXct(paste(YMD,DAT$time),origin="1970-1-1 00:00:00")

消費カロリーの計算

　推定されたMETsを使って消費カロリーを求めたいときは，

　エネルギー消費量 [kcal] = METs × (時間 [h]) × (体重 [kg])

を使います．インターネットで「METs 消費カロリー」検索すると，

　エネルギー消費量 [kcal] = METs × (時間 [h]) × (体重 [kg]) × 1.05

という式が多く見つかると思います．最初の式と「× 1.05」の部分が違います．

　どちらも間違いではなく，どちらも近似的な式で，どちらもほとんど値は同じです．「× 1.05」の部分は，体表面積などの影響で調整するための値を，標準的値で近似したものです．体型によって体表面積は異なるので，本当は個人個人で変えた方が良いかもしれません．とはいえ，年齢，性別など他の因子にも依存するので，しょせん近似です．

　運動生理やスポーツ医学の業界では，「× 1.05」は省略してもいいでしょ，ということにしたようです．この論文で， [Ainsworth et al. Medicine & science in sports & exercise 43 (2011): 1575-1581]，「× 1.05」なしの式が推奨されています．

　ということで，METsデータから消費カロリーを計算してみましょう，と言うのが宿題です．

2022-09-26

充電用USB type-Cケーブル

スマホやパソコンを充電したり，周辺機器をつないだりするケーブルとして，USB type-Cが一般的になっています．USB type-Cのケーブルって見た目は同じなのに，機能や性能が違うことを最近知りました．ということで，今回は「充電用のUSB type-Cケーブル」について調べてみました．

Type-Cで2m超の長いケーブルはUSB 2.0規格
データ通信はUSB2.0だけど，急速充電可能な充電用USB type-Cケーブルがある
60Wケーブルと100Wケーブルの違い
まとめ

Type-Cで2m超の長いケーブルはUSB 2.0規格

　データ転送の早さが売りのUSB type-Cケーブルは，規格上，長さの上限があります．例えば，「USB 3.1 Gen1」では2ｍまで，「USB 3.1 Gen2」や「USB 3.2」では1ｍまでとなっています．ですので，2mを超えるUSB type-Cケーブルは，今となっては古い「USB 2.0」規格の可能性大です．つまり，データの転送速度は威張れるほどではありません．でも，充電の性能は悪くないかもしれないというのが今回のポイントです．

　USB type-Aからtype-Cへのケーブルだと，USB type-Aのコネクタの内部が黒くなっていれば，USB 2.0規格です．昔ながらのUSB 2.0規格準拠だと，最大5V×0.5A = 2.5 Wになるので，最新の急速充電はできません．が，しかし，USB type-Aからtype-Cへのケーブルでも，"USB Type-C Current"という規格で充電できるようで，5V×3A = 15 Wくらいいけるそうです．

一般的なUSB type-Aからtype-Cへのケーブルを，USB CABLE CHECKER 2 (Bit Trade One製)で評価．ディスプレイに「SOURCE 0.5A」とあるように，0.5Aの電流が上限として設定されるので，供給電圧が5VのUSB2.0規格では，5V×0.5A = 2.5 Wの給電が最大になる．と想像しましたが，Anker PowerPort Atom III 65W SlimでUSB type-Aからtype-c接続でスマホを充電して，計ってみると2A近く (10Wくらい)流れていました．ケーブルの説明を読むと，「3A急速充電」ということで，最大3Aまで大丈夫と書いてありました．想像では，"USB Type-C Current"という規格なんだと思います．

データ通信はUSB2.0だけど，急速充電可能な充電用USB type-Cケーブルがある

　私は，ノートパソコンやスマホへの充電に，AnkerとMcdodoのケーブルを使っています．Ankerだと「PowerLine III Flow USB-C & USB-C ケーブル」，Mcdodoだと「Type C Type Cケーブル」（コネクタに電力が表示されるやつ）を主に使っています．理由は，「USB Power Delivery (略してUSB PD)」という急速充電の規格に対応していると明記してあることと， 100Wまでいけると給電電力も明記してあるので，なんとなくよさそうと感じたからです．これって，私の想像です．

Anker PowerLine III Flow USB-C & USB-C ケーブルをUSB CABLE CHECKER 2 (Bit Trade One製)で評価．ケーブル内の配線はほぼUSB2.0だけと，Configration Channel (CC)の配線があり，E-Markedケーブルになっている．ディスプレイ「E-Marked」が表示されているので，eMarkerチップがケーブルに内蔵されていることが分かる．E-Markedケーブルにしてある理由は，最大電流5Aの100Wケーブルにするためです．

Mcdodo Type C Type CケーブルをUSB CABLE CHECKER 2 (Bit Trade One製)で評価．これも，ケーブル内の配線はほぼUSB2.0だけと，Configration Channel (CC)の配線があり，E-Markedケーブルになっている．

　今回調べてわかったことは，私が充電に使っているUSB type-Cは，データ通信については，USB2.0規格だということです (上の写真)．しかし，給電については，「USB 3.1 Gen. 1」以上の規格 (あるいは，USB PDの規格)で定められた，eMarkerというチップが内蔵されています．eMarkerチップが内蔵されたケーブルは，「E-Markedケーブル」と呼ばれます．ここで紹介したケーブルでは，eMarkerチップに書かれた情報により，最大5Aまで流せるケーブルとして認識されます．現行のUSB PD規格では，最大電圧が20Vなので，5A×20V=100Wの供給が可能になり，高速充電が実現できるようです．なんか，充電に特化したケーブルのようです．下の写真のように，Mcdodo Type C Type Cケーブルを使ってノートパソコンを充電してみると，USB PDに対応して60W供給しているようなので，古いUSB2.0の2.5Wよりは高速充電できてます．今回，ノートパソコンの充電には，最大65W出力のUSB PD充電器Anker Nano II 65Wを使いました．

USD PD充電器Anker Nano II 65WとMcdodo Type C Type Cケーブルでパソコンを充電．Mcdodoのケーブルのコネクタ部分に供給されている電力が表示される．

　PD対応で100W急速充電可能なType-Cケーブルは，USB PD充電器を使った充電用のケーブルです．しかし，ディスプレイの接続には使えないし，データ通信の速度も最高速ではないということです．とはいえ， USB2.0と同様にデータ通信はできます．

　USB3.x規格では，ケーブルが短く，充電用途には不便なので，USB2.0のtype-Cケーブルが存在するのだと思います．これは，私の想像です．

USB 3.1 Gen2規格のケーブルをUSB CABLE CHECKER 2 (Bit Trade One製)で評価するとこんな感じになる．

60Wケーブルと100Wケーブルの違い

　両端がType-Cのケーブルは，規格上，60Wの給電に対応しているそうです．つまり，規格を無視した粗悪品は除いて，どのType-Cケーブルでも，最大3Aの電流が流せて，USB PD規格で最大20Vをかけたとすれば，3A×20V=60Wの供給が可能ということです．

　60Wを超えて100Wの供給が可能なType-Cケーブルでは，eMarkerチップに「最大5Aまで流せるよ」という情報が書いてあるそうです．

　USB2.0規格のType-Cケーブルでは，eMarkerチップがそもそも内蔵されていないものがあります．その場合は，充電時に最大電流3Aの60Wケーブルと判定されます．ただし，60Wケーブルでも，Type-Cの仕様であるConfigration Channel (CC)という配線が存在する必要があって，その有無がUSB CABLE CHECKER 2のCCの点灯で分かります．両端がType-Cでないと，CCの配線は作れません．

　上で紹介したE-Markedケーブルでは，CCもあり，eMarkerチップも内蔵されていて，充電時に100Wケーブルと判定されるということです．

まとめ

　今回は気になっていたUSB type-Cことを調べてみました．分かったことは，充電用の長めのUSB type-Cケーブル (USB 2.0)と，データ転送用のUSB3.x仕様のUSB type-Cケーブルを両方用意して，用途に応じて使い分けた方がよいかもということです．とはいえ，データ保存に使う外付けSSDを買うと，短いUSB type-Cがついているので，データ転送用のケーブルをわざわざ買う必要性は低いです．

　私は，USBケーブルと充電の専門家ではないので，間違いがあると思います．間違いがあれば教えてください．

2022-09-24

【Rで時系列解析】cosinorパッケージでコサイナー解析，cosinor2も必要

R言語

最近はYoutubeでマリマリマリーの動画を見ています．その影響で，ひ○ゆきのものマネができたらいいなーと思うようになりました．皆さんも，次のセミナーではひ○ゆき風のしゃべりでお願いします．

　ということで，今回は前に紹介したコサイナー解析
chaos-kiyono.hatenablog.com
の続編です．以前の私の記事ではRのパッケージは使いませんでしたが，今回はcosinorパッケージを使います．ただし，cosinorパッケージだけでは，使い勝手が良くないので，cosinor2パッケージも一緒に使います．

周期的リズムのモデル

　ここでは，周期 $T$ を既知として，以下のモデルを仮定します．

　 $\displaystyle y(t) = \mu + A \cos \left(\frac{2 \pi t}{T} + \phi \right) + \varepsilon(t)$

ここで， $\mu$ は，MESOR (midline statistic of rhythm)， $A$ は振幅， $\phi$ はacrophase (頂点位相)と呼ばれるパラメタです． $\varepsilon(t)$ は，ランダムな要因やその他の変動を表します．前の記事では $\phi$ を $-\phi$ としましたが，今回は [Bingham et al. (1982)] に合わせて $+\phi$ )としてあります．

　コサイナー解析では，このモデルのパラメタを推定します．

cosinorとcosinor2パッケージでコサイナー解析

　cosinorパッケージを使ったコサイナー解析の問題点は，acrophase (頂点位相) $\phi$ が自動的に計算できないことです．これは， $\phi$ の推定にアークタンジェントを使っているからだと思います．数学では，アークタンジェントの値域は， $-\pi/2 < \phi < \pi/2$ (-90°から90°)です．しかし，これでは，225°が45°になる問題があるということです．

　cosinorパッケージの関数cosinor.lm(...)の結果をみれば，角度を正しく修正できますが，それでは面倒です．なので，私はcosinorパッケージを使うことはありませんでした．しかし，2022年にcosinor2というパッケージが登場して，この面倒をなくしてくれました．

　使い方の例として，Rスクリプトを載せておきます．

# cosinorパッケージを使った解析
#　【注意】　事前にcosinor, cosinor2をインストール
# install.packages("cosinor")
# install.packages("cosinor2")
#############################
# パッケージを読み込み
require(cosinor)
require(cosinor2)
#############################
# 数値実験
###
# 周期
T <- 24
# テスト用の数値データの生成
time <-  seq(0,24,1)
mesor <- 1
amp <- 3
phase <- 2.4
eps <- rnorm(length(time))
y.obs <- mesor+amp*cos(2*pi*time/T+phase) + eps
#####################
# プロット
par(mfrow=c(1,1),pty="m")
plot(time,y.obs,xlab="time",col=4,cex=1.2,xaxs="i")
curve(mesor+amp*cos(2*pi*x/T+phase),add=TRUE,col=8,xlim=c(0,24),lwd=2)
fit.cos <- cosinor.lm(y~time(t), period = 24, data=data.frame(t=time,y=y.obs))
para <- coef(fit.cos)
# paraの中身は
# mesor: para[[1]]
# amp: para[[2]]
# acr: para[[3]]
########################################
# 【重要】 cosinor2でアクロフェイズを補正
para[3] <- correct.acrophase(fit.cos) 
########################################
curve(para[[2]]*cos(2*pi*x/24+para[[3]])+para[[1]],xlim=c(min(time-12),max(time+12)),add=TRUE,col=2,lwd=3,lty=2)
para

2022-09-21

【ラグオペレータで高校数学を解いてみる】3項間漸化式と連立漸化式

最近，高校の数学で出てくる漸化式

　 $a_{n+2}=p\, a_{n+1}+q\, a_n$

や，

　 $\left\{ \begin{array}{l} a_{n+1} = p\, a_n + q\, b_n \\ b_{n+1} = r\, a_n + s\, b_n \end{array} \right.$

の解き方を教える機会がありました．私にとってはかなり昔の記憶ですが，高校の数学では，公式的なパターンに変形したり，行列の形にして対角化するとかありました．ラグオペレータでも解けるので，今回はラグオペレータで解いてみます．

ラグオペレータの基本

　ラグオペレータを $L$ とします． $L$ を左からかけると，時間が1だけ戻ります．つまり，

　 $L a_n = a_{n-1}$

となります．さらに，

　 $L^m a_n = a_{n-m}, \quad L^{-1} a_n = a_{n+1}$

となります．

３項間漸化式の例題

　高校数学で登場するような問題を考えます．例題として，

　 $a_{n+2}=3\, a_{n+1}-2\, a_n$

で， $a_{1}=2, a_{2}=3$ となる $a_n$ を求めてみます．

　 $n$ を2減らして，

　 $a_{n}=3\, a_{n-1}-2\, a_{n-2}$

としてから，ラグオペレータを使えば，

　 $a_{n}=3 L a_{n} - 2 L^2 a_{n}$

です．さらに，変形すると，

　 $\begin{eqnarray} a_{n} - 3 L a_{n} + 2 L^2 a_{n} &=& 0 \\ \left(1-3L+2L^2 \right) a_{n} &=& 0 \\ \left(1-L\right)\left(1-2L\right) a_{n} &=& 0 \end{eqnarray}$

となります．左辺が0になるとすれば，　

　 $\left\{ \begin{array}{l} \left(1-L\right) a_{n} &=& 0 \\ \left(1-2L\right) a_{n} &=& 0 \end{array} \right.$

であれば良いことが分かります．ラグオペレータを使わないで表すと，この式は，

　 $\left\{ \begin{array}{l} a_{n} - a_{n-1} &=& 0 \\ a_{n} - 2 a_{n-1} &=& 0 \end{array} \right.$

ってことです．ここから，

　 $\left\{ \begin{array}{l} a_{n} &=& a_{n-1} = a_1 = 2 \\ a_{n} &=& 2 a_{n-1} = 2^{n-1} a_1 = 2^{n-1} \cdot 2 = 2^n \end{array} \right.$

と求まります．「解が2つ！」とびっくりしないでください．まだ，解ではありません．線形性から (過去の記事参照)，一般解は，2つの解の線形結合になります．つまり， $A$ と $B$ を定数として，一般解は

　 $a_n=A \cdot 2 + B\cdot 2^{n}$

の形になります．

　初期条件 $a_{1}=2, a_{2}=3$ を満たすように， $A, B$ を求めると，

　 $\displaystyle A=\frac{1}{2}, \quad B=\frac{1}{2}$

となります．ということで，解は，

　 $a_n=1+2^{n-1}$

です．

連立漸化式の例題

　例題として，

　 $\left\{ \begin{array}{l} a_{n+1} = 3\, a_n + b_n \\ b_{n+1} = 2\, a_n + 2\, b_n \end{array} \right.$

で， $a_{1}=1, b_{1}=-1$ となる $a_n$ を求めてみます．

　 $n$ を1減らしてから，ラグオペレータを使えば，

　 $\left\{ \begin{array}{l} a_{n} = 3 L a_n + L b_n \\ b_{n} = 2 L a_n + 2 L b_n \end{array} \right.$

　 $\left\{ \begin{array}{l} (1-3L)a_{n} - L b_n = 0 \\ 2 L a_n + (2 L -1 ) b_n = 0 \end{array} \right.$

となります．連立方程式を解いて， $b_n$ を消すと，

　 $\begin{eqnarray} \left(4L^2 - 5 L + 1 \right) a_{n} &=& 0 \\ \left(4L - 1 \right)\left(L - 1 \right) a_{n} &=& 0 \end{eqnarray}$

となります．左辺が0になるとすれば，　

　 $\left\{ \begin{array}{l} \left(4 L -1\right) a_{n} &=& 0 \\ \left(L - 1 \right) a_{n} &=& 0 \end{array} \right.$

であれば良いことが分かります．ラグオペレータを使わないで表すと，この式は，

　 $\left\{ \begin{array}{l} a_{n} &=& 4 a_{n-1} \\ a_{n} &=& a_{n-1} \end{array} \right.$

ってことです．ここから，

　 $\left\{ \begin{array}{l} a_{n} &=& 4 a_{n-1} = 4^{n-1} a_1 = 4^{n-1} \\ a_{n} &=& a_{n-1} = a_1 = 1 \end{array} \right.$

と求まります． $A$ と $B$ を定数として，一般解は

　 $a_n=A \cdot 4^{n-1} + B \cdot 1$

の形になります．

　初期条件 $a_{1}=1, a_{2}=3a_1+b_1 = 3 - 1 = 2$ を満たすように， $A, B$ を求めると，

　 $\displaystyle A=\frac{1}{3}, \quad B=\frac{2}{3}$

となります．ということで， $a_{n}$ は，

　 $\displaystyle a_n= \frac{4^{n-1}+2}{3}$

です． $b_{n}$ は，

　 $\begin{eqnarray} b_n &=& a_{n+1} - 3 a_{n} \\ &=& \frac{4^{n}+2}{3} - 3 \cdot \frac{4^{n-1}+2}{3} \\ &=& \frac{4^{n-1}-4}{3} \end{eqnarray}$

まとめ

　ラグオペレータを使って高校数学の漸化式の問題を解いてみました．計算は自信がないので，間違いを見つけたら教えてください．

2022-09-20

【心拍変動解析の基礎】RR間隔時系列の修正

心拍変動解析では，通常の正常洞調律 (normal sinus rhythm)の間隔時系列を分析します．正常洞調律についての説明は，以前の記事を参照してください．
chaos-kiyono.hatenablog.com

　心拍変動では，正常洞調律を分析するというのが建前とはいえ，多くの場合，R波の検出不良があったり，洞調律ではない不整脈が含まれていたりします．心拍変動解析では，そーいった非洞調律によるRR間隔データを，正常洞調律っぽく修正してから分析します．今回は，心拍変動解析のための基礎知識として，心電図波形とRR間隔の関係について説明します．

まず，正常洞調律

　下の図上段のような正常洞調律の心電図 (ECG)波形を見ると，同じ波形パターンが繰り返し現れます．でも，R波の出現間隔 (図下段のRR間隔)は，メトロノームのように一定ではなく，ある程度ばらつきます．

短長パターンの典型の心室性期外収縮

　心室性期外収縮 (ventricular premature contraction：VPC)は，心室が調和を乱して，自分勝手に発火 (電気発生)してしまうことで生じる不整脈です．下の図上段の心室性期外収縮 (VPC)の部分にみられるパターンです．RR間隔が自動検出されたデータでは，心室性期外収縮 (VPC)の部分のRR間隔は，異常に短くなります．その直後のRR間隔は，異常に長くなります．ここで，異常とは，正常洞調律ではありえないということです．

　心室性期外収縮 (VPC)の場合，心室は拍子を外していますが，洞結節からの電気の流れは正常な間隔を反映しています．ですので，心室性期外収縮 (VPC)の前後の，短長の２つのRR間隔を足したものが，正常洞調律のRR間隔2拍分に等しいと考えられます．

　ということで，心室性期外収縮 (VPC)のRR間隔の修正は，

　 $\displaystyle {\bf 修正後の}{\rm RR}_3 = {\bf 修正後の}{\rm RR}_4 = \frac{{\rm RR}_3+{\rm RR}_4}{2}$

とすることが多いです．

ノイズを誤検出

　下の図は，ノイズで発生したピークを誤検出し，真のR波を見逃した例です．R波の自動検出アルゴリズムでは，心筋の不応期 (発火すると300msくらい休む特性)を考慮して，R波の間隔が近くならないようにしてあることが多いです．ですので，ノイズピークの後のR波は飛ばされる可能性があります．

　上の図の場合も，心室性期外収縮と同じようにRR間隔を修正すれば大丈夫です．

　ただし，通常のR波とR波の間にノイズのピークが検出された可能性がある場合は，２つのRR間隔を足して，１つ分のRR間隔にしてください．

R波の検出漏れ

　心電計測で，R波の高さが低くなったり，ノイズに隠れて検出できないことがあります．下の図の例では，２個のR波が飛ばされて言います．

　この場合は，直前の ${\rm RR}_2$ と ${\rm RR}_3$ の比から，おおよそ， $1:3$ とわかるので， ${\rm RR}_3$ はRR間隔３個分だと予想できます．ということで，修正は， ${\rm RR}_3$ を3で割ったものを， ${\rm RR}_3$ の代わりに3つ挿入すれば大丈夫です．

上室性期外収縮

　心房側が勝手なふるまいをするのが，上室性期外収縮 (supraventricular premature contraction: SVPC)です．下の図のように上室性期外収縮 (SVPC)は，修正が難しいケースです．なぜなら，心室性期外収縮 (VPC)のように足して２で割る方法がイマイチだからです．それでも，足して２で割る修正が，選択肢の一つです．

修正不可能な心房細動と心房粗動

　下の図は，上段が心房細動，下段が心房粗動の例です．

　心房細動では，1.5秒を超えるような長いRR間隔がむちゃくちゃたくさん含まれていて，非常にノイジーなグラフになります．心房粗動は，逆に結構周期的に見えて，RR間隔のグラフに数本の筋が見えたりします．患者さんのデータを分析するときは，心房細動，心房粗動，ペースメーカなどのデータが含まれていることがあるので，洞調律とは分けで分析するようにしてください．

学生のための宿題

　以上の知識を踏まえて，RR間隔データを自動修正するプログラムを各自で作ってください．

ケィオスの時系列解析メモランダム

時系列解析，生体情報学，数学・物理などの解説です．

【Rでの時刻データの扱い】as.POSIXctを使うの卒業します

`fasttime`パッケージのインストール

処理時間の比較

【重要】日本時間での使い方

まとめ

【Rで日付，時間を扱う】as.Dateを使うときはタイムゾーンを必ず指定しろ

`as.Date`を使うときは`tz`を指定しろ

日付がずれることの確認

まとめ

【Rで活動量分析】オムロン研究用活動量計HJA-750Cのデータファイルの読み込み

活動量評価のデータファイルの中身

活動量ファイルの読み込み

消費カロリーの計算

充電用USB type-Cケーブル

Type-Cで2m超の長いケーブルはUSB 2.0規格

データ通信はUSB2.0だけど，急速充電可能な充電用USB type-Cケーブルがある

60Wケーブルと100Wケーブルの違い

まとめ

【Rで時系列解析】cosinorパッケージでコサイナー解析，cosinor2も必要

周期的リズムのモデル

cosinorとcosinor2パッケージでコサイナー解析

【ラグオペレータで高校数学を解いてみる】3項間漸化式と連立漸化式

ラグオペレータの基本

３項間漸化式の例題

連立漸化式の例題

まとめ

【心拍変動解析の基礎】RR間隔時系列の修正

まず，正常洞調律

短長パターンの典型の心室性期外収縮

ノイズを誤検出

R波の検出漏れ

上室性期外収縮

修正不可能な心房細動と心房粗動

学生のための宿題

fasttimeパッケージのインストール

処理時間の比較

【重要】日本時間での使い方

まとめ

as.Dateを使うときはtzを指定しろ

日付がずれることの確認

まとめ

活動量評価のデータファイルの中身

活動量ファイルの読み込み

消費カロリーの計算

Type-Cで2m超の長いケーブルはUSB 2.0規格

データ通信はUSB2.0だけど，急速充電可能な充電用USB type-Cケーブルがある

60Wケーブルと100Wケーブルの違い

まとめ

周期的リズムのモデル

cosinorとcosinor2パッケージでコサイナー解析

ラグオペレータの基本

３項間漸化式の例題

連立漸化式の例題

まとめ

まず，正常洞調律

短長パターンの典型の心室性期外収縮

ノイズを誤検出

R波の検出漏れ

上室性期外収縮

修正不可能な心房細動と心房粗動

学生のための宿題

`fasttime`パッケージのインストール

`as.Date`を使うときは`tz`を指定しろ