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

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

R言語

【長時間相関解析】DFAでホワイトノイズを分析するとどうなるか(つづき):数値実験

Detrended Fluctuation Analysis (DFA)で,ホワイトノイズを解析したらどうなるか,数値実験で確かめてみます.理論的な結果は,前の記事を参照してください. chaos-kiyono.hatenablog.com 数値時系列の生成 1次DFA 2次DFA 数値時系列の生成 Rでホワイトノ…

【長時間相関解析】DFAでホワイトノイズを分析するとどうなるか

Detrended Fluctuation Analysis (DFA)で,無相関な信号,つまり,ホワイトノイズを解析したらどうなるかのお話です.DFAに詳しい人の多くは,そんなのスケールをとして,ゆらぎ関数 (fluctuation function)は, だろというかもしれません.が大きい,漸近的…

【Rで確率モデル】疑似乱数の生成

昨日,名古屋に出張しました.新幹線は,コロナ前と比べればまだ空いていますが,海外の旅行者も戻ってきて,席が埋まってきたように感じます.ということで,今回は,Rで乱数を生成するコマンドの整理です.疑似乱数を使って下の図のような,中心極限定理の…

【Rで確率モデル】確率密度関数と確率分布関数

世の中には正規分布とか,指数分布とか,ワイブル分布とかいろいろな分布に従う現象があります.確率的な現象が従う分布の関数形を特定すると,背後にある生成メカニズムを考えるヒントになるので,分布のことを理解しておくと役に立ちます.今回は,Rで確率…

【Rで高速フーリエ変換】バタフライダイアグラムを読む

私は肉まんの白い部分が好きです.ほのかに甘い,あの白い部分だけ食べたいです. ということで,今回は最近のFFTのお話 chaos-kiyono.hatenablog.com の続きです. FFTの説明で,バタフライ演算ってのが登場します.バタフライ演算は下の図のようなバタフラ…

【Rで高速フーリエ変換】2のべき乗の長さの時系列用FFTの実装についての解説

前回,2のべき乗の長さの時系列を,高速に離散フーリエ変換 (FFT)できるアルゴリズムのRスクリプトを作りました (高速じゃないけど). chaos-kiyono.hatenablog.com このFFTアルゴリズムは,Cooley-Tukey型と呼ばれるそうです (知らんけど).今回はこのアル…

【Rで高速フーリエ変換】Rのfftは任意の長さで大丈夫だけど,学習用に基数2のFFTを自作

Rとか,Pythonとかを使って,誰でも簡単に高速フーリエ変換 (fast Fourier transform, FFT)を使える時代になりました.最近コロナで,人と会わなくなって,久しぶりに遠隔会議であった方々が老けたなーと感じることが多くなりました.同様に私も年寄りになっ…

【Rの処理時間】計算時間短縮のための工夫 ― meanよりもsumが速い ―

プログラムを作って何か処理するとき,処理時間が短い方が良いに決まってます.ということで,今回はRの処理時間についてのお話です. 処理時間を無駄に長くしないための一般論として,私が説明できることは,「forループは使わずにベクトルのまま計算した方…

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

1時間ほど前に書いた記事では,Rで文字列を時刻データに変換するのに「私はas.POSIXctを使います」と言ってましたが,訂正します. 私はas.POSIXctを2度と使いません. 今なら断言できます.Rでas.POSIXctを使ってる奴,素人です. as.POSIXctを使ってる人…

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

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

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

オムロン活動量計でHJA-750Cというのがあります.www.healthcare.omron.co.jpこれは研究や臨床検査用なので,一般には発売されていません.大学などの研究機関に所属していれば購入可能です. といことで,今回はHJA-750Cを使って研究をしたい人向けのに,R…

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

最近はYoutubeでマリマリマリーの動画を見ています.その影響で,ひ○ゆきのものマネができたらいいなーと思うようになりました.皆さんも,次のセミナーではひ○ゆき風のしゃべりでお願いします. ということで,今回は前に紹介したコサイナー解析 chaos-kiyo…

【Rでデータ処理】データテーブル操作の基礎

今回はファイルからデータテーブル (表)として読み込んで,データを処理する方法をいくつか紹介します.データ処理といっても,条件にあう列や行を抽出したりするだけです. ここでは,具体例として,ユニオンツールのmyBeatの心拍間隔データを使います.フ…

【Rで時系列解析】ハイパス,ローパス,バンドパス:バターワース (Butterworth filter)フィルタのかけ方

信号処理では,時系列に含まれる特定の周波数成分を取り出すために,ローパス (low pass)フィルタとか,ハイパス (high pass)フィルタとか,バンドパス (band pass)フィルタとかを使うことがあります.学生の皆さんには,デジタルフィルタの基礎の理論から実…

【Rでファイルの一括処理】正規表現とワイルドカードの活用

前回,フォルダ内にあるたくさんのファイルについて,フォルダ内にある全ファイル名を取得するRのコマンドとしてlist.files(...)の使い方を紹介しました.今回は,ファイル名に含まれる文字を使って,特定のファイルを見つける方法を説明します.今回,メイ…

【Rでファイルの一括処理】フォルダ内にある条件に合うファイルの取得

コンピュータを使えば,大量のデータを一括処理できます.つまり,たくさんのデータファイルがあっても,それらを自動で読み込んで,分析できちゃいます.今回は,Rで,たくさんのデータファイルを処理するときのヒントです. 事前知識 Rでフォルダ内にある…

【Rで時系列解析】コサイナー解析 (cosinor analysis)

ひさびさに,共同研究でコサイナー解析 (cosinor analysis)を頼まれたので,ポイントをメモっておきます.コサイナー解析というのは,周期的に変化すると思われる時系列に,コサインの波形を当てはめてる分析法です.概日リズムの特徴付けに使われることがあ…

【信号処理の基礎】Savitzky-Golay平滑化・微分フィルタの周波数特性

今回は,Savitzky-Golayフィルタの畳込み表現の係数 (下図のピンク図形)と周波数特性の話です. Savitzky-Golay平滑化フィルタ.元の時系列(上段灰色)は,下段の破線にノイズを加えたもの.下段の灰色破線はノイズを加える前の時系列.赤実線(上下両方)…

【Rで時系列解析】ヒルベルト変換でエンベロープを計算するときの注意

前回,振動している信号のエンベロープ (包絡線)をヒルベルト変換で求める方法を紹介しました.今回は,この方法がうまくいくための注意事項をいくつか書いておきます. ポイントは, 前処理として,周波数の狭いバンドパスをかけておく あるいは 前処理とし…

【Rで時系列解析】ヒルベルト変換で振動のエンベロープを抽出

振動する時系列の中には,振動の振幅の変化に役立つ情報がある場合があります.例えば,呼吸信号(換気量や気流)は振動します.そのとき,心拍数に影響を与えるのは,振動の周期ではなく,振幅です.ということで,今回はその「振幅」の情報を取り出すため…

【Rで時系列解析】AICで多変量自己回帰過程の次数を決定

今回は,2変量の自己回帰過程 について最小2乗法でパラメタを推定し,AICで次数を決定します.お気楽な"vars"パッケージのVARなどは使いません.計算過程や式を確認できるように,Rスクリプトを書きました.Rスクリプトは一番下に掲載してあります.AICにつ…

【Rで時系列解析】AICを使った自己回帰過程の次数推定

今回は,時系列に対して自己回帰過程をフィットし,最適な次数を決定します.自己回帰過程の次数の推定には,赤池情報量規準 (AIC: Akaike Information Criterion)を使います.AICの導出をここで説明することは,私には難しいのでやりません.でも,計算過程…

【Rで時系列解析】コヒーレンス (coherence)の計算

今回は,コヒーレンス (coherence,あるいは,magnitude-squared coherence)の計算のメモです.コヒーレンスは,2つの時系列間の相関を周波数領域で見る方法です.2つの信号の周波数成分に,線形な入出力関係があれば1に近くなります.なければ0です.とはい…

【Rで時系列解析】最小2乗法で自己回帰過程の当てはめ

観測された時系列に最小2乗法で自己回帰過程を当てはめる計算のメモです.どんな形の過程を仮定しても,基本的な考え方を理解しておけば,係数を求めることができます.楽をしたければ,arとか,varsパッケージなどを使えば良いと思います.しかし,お気楽…

【Rで長時間相関】ARFIMA(0, d, 0)のAR過程近似を数値的に検証

前回に続き,ARFIMA(0, d, 0)を,無限次の自己回帰過程 (AR過程: autoregressive process)で表現する話です.今回は数値実験です. 1. 基本事項:long AR過程表現 2. 数値計算のポイント 3. 数値計算結果 Rスクリプト 前回の記事は,これです. chaos-kiyono…

【Rで時系列解析】Savitzky-Golayフィルタで平滑化

今回はギザギザ,凸凹したデータの平滑化法の一つであるSavitzky-Golayフィルタの紹介です. Savitzky-Golayフィルタの考え方 (2次でスケール5の例):部分区間に多項式をフィットして中央の値を平滑化された値として採用する.実用上は,畳み込み和を使って…

【Rで時系列解析】グレンジャー因果 (Granger causality)の周波数領域表現を数値実験で検証

グレンジャー因果 (Granger causality)のGewekeのアプローチを少しずつ理解していきます.前回の説明の続きで,今回は数値実験をします. chaos-kiyono.hatenablog.com 1. 問題設定 単純な例から考えはじめないと,私には理解できません.一般的な表現を見て…

【Rで時系列解析】Rでパワースペクトル推定:心拍変動解析での実際

Rで時系列のパワースペクトルを推定したいとき,もっともお気楽なのは,"spectrum"というコマンドを使うことだと思います.今回は,心拍変動 (RR間隔時系列)のスペクトル解析で"spectrum"コマンドを使った実例を紹介します.ちなみに,例として使用したRR間…

【Rで時系列解析】Rでspectrumを使ったパワースペクトル推定におけるspansの設定(数値デモ)

Rのコマンド"spectrum"のオプション"spans"の設定で,パワースペクトルの推定結果がどうなるかのデモです."spans"が何かは,前回の記事を参考にしてください. chaos-kiyono.hatenablog.com 時系列の生成については,以下の記事を参考にしてください. chao…

【Rで時系列解析】スペクトルウインドーを使ったパワースペクトルの平滑化:Rの"spectrum"での"spans"オプションの意味

今回は,離散フーリエ変換でえられるギザギザなパワースペクトルを,滑らかにする「平滑化」についてです.いろいろ書きたいことはありますが,まずは,Rでお世話になることが多い"spectrum"というコマンドのオプション"spans"について説明します.今回の説…