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

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

【Rで呼吸代謝分析】ミナト医科学AE-100iデータの読み込み

今回は,ミナト医科学のAE-100iで計測されたデータファイルの読み込みについてです.医療関係者と研究者しか使わない計測器だと思いますが,この計測器を使わない人はデータファイルを読むときの参考にしてください.

計測開始時刻を抽出

 AE-100iで計測されたデータファイルの最初の2行は以下のようになってます.

FileNo IDcode Days Time ...
123 hogehoge 10/28/2022 13:34:37 ...
... ... ... ... ...

この2行目に,計測日時"10/28/2022" (2022年10月28日)と,開始時刻"13:34:37"が書いてあるので,その部分を抽出して,開始日時を取り出します.そのためのRスクリプトが以下です.

# ファイル名の指定
FN.RESP <- "ファイル名.csv"
tmp <- read.csv(FN.RESP, header=FALSE, skip=1, stringsAsFactors=FALSE, nrows=1)
TIME.start <- strptime(paste(tmp$V3,tmp$V4),"%m/%d/%Y %H:%M:%S")

このスクリプトのポイントは,strptimeの使い方です.日本では,日付を表すとき,"2022/10/28"のように年月日の順番が一般的です.しかし,海外では,

 "10/28/2022"のように月日年の順番や,

 "28/10/2022"のように日月年の順番で書くこともあります.

そのようなデータを読み込むときには,strptime(..., "%m/%d/%Y %H:%M:%S")のように

%Y 年 (4桁)
%m
%d

の順番を指定して,日付を変換してください.

 上の例では日付の区切り文字が/ですが,"2022-10-28"の場合は,"%Y-%m-%d"のように指定してください.

 時刻については,

%H
%M
%S

を使って形式を指定してください.

METSデータの抽出と計測時刻の設定

 以下のスクリプトを実行すれば,呼気ガス計測で計算したMETSとその計測時刻TIMEが各変数に入ります.

# ファイル名の指定
FN.RESP <- "ファイル名.csv"
# 計測開始日時の情報を抽出
tmp <- read.csv(FN.RESP,header=FALSE,skip=1,stringsAsFactors=FALSE,nrows=1)
TIME.start <- strptime(paste(tmp$V3,tmp$V4),"%m/%d/%Y %H:%M:%S")
# 代謝等量METSの抽出
METS <- read.csv(FN.RESP,header=FALSE,skip=3,stringsAsFactors=FALSE)[,"V12"]
n.tmp <- length(METS)
# 計測時刻
TIME <- TIME.start+10*((1:n.tmp)-1)
# 結果のプロット
plot(TIME,METS,"l",col=2,xaxs="i")

まとめ

 冬になると,深夜に歩きとか自転車で,自宅に帰るとき寒いので,早めに冬着を用意しました.とはいえ,私の気持ちでは早めの準備でしたが,世間的には出遅れたようです.例えば,ダウンジャケットについては,デサントの水沢ダウンとか,ノースフェイスを買おうと思いましたが,通販ではほぼ売り切れていました.好みの色とサイズを買うには,9月には買わないといけないようです.もうちょっと待つと再入荷するかもしれませんが,とりあえず,暖かそうなダウンジャケットを一つ買いました.

 調子にのって高級なものを買っても私には無駄だと思い,ほとんどの防寒着は,昨日,近所のワークマンで買いました.ワークマンでも,私のサイズは,どれも残り1点しかありませんでした.安いので,今年の新着を一通り買いました.これで,寒さを楽しめそうです.

【フーリエ変換】フーリエ級数からフーリエ変換へ

今回はフーリエ変換の話です.ここでは横軸をxではなく,tにします.理由は,時間tとともに変化する信号x(t)を考えたいからです.

複素フーリエ級数

 周期Tの周期関数

 x(t) = x(t + T)

フーリエ級数展開は,

 \begin{eqnarray}x(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos \frac{2 \pi n  t}{T} + b_n \sin \frac{2 \pi n  t}{T} \right)
\end{eqnarray}

です.ここで,係数の\{a_0, a_1, b_1, a_2, b_2, \cdots\}は,

 \begin{eqnarray}a_n &=& \frac{2}{T} \int_{-T/2}^{T/2} \! x(t) \cos \frac{2 \pi n  t}{T} \, dt \\
b_n &=& \frac{2}{T} \int_{-T/2}^{T/2} \! x(t) \sin \frac{2 \pi n  t}{T} \, dt
\end{eqnarray}

で与えられます.

 フーリエ級数展開の各成分の周波数 f_n は,

 \displaystyle f_n = \frac{n}{T}

になってます.

 次に,x(t)や,\{a_0, a_1, b_1, a_2, b_2, \cdots\}複素数に拡張してみます.複素数の世界では,三角関数が指数関数を使って,

 \begin{eqnarray} \cos \theta &=& \frac{e^{i \theta} + e^{- i \theta}}{2} \\
\sin \theta &=& \frac{e^{i \theta} - e^{- i \theta}}{2}
\end{eqnarray}

のように与えられます.この関係を使えば,

 \begin{eqnarray}x(t) &=& \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos \frac{2 \pi n  t}{T} + b_n \sin \frac{2 \pi n t}{T} \right) \\
&=& \sum_{n=-\infty}^{\infty} c_n \exp \frac{2 \pi i n t}{T}
\end{eqnarray}

と変形できます.ここで,

 \begin{eqnarray} c_n &=& \frac{1}{T} \int_{-T/2}^{T/2} \! x(t) \exp \left(- \frac{2 \pi i n t}{T} \right) \, dt \end{eqnarray}

です.

 フーリエ級数展開では,下の図のような周期関数を想定していました.

フーリエ級数で考える周期関数の例.

 現実世界で観測される信号は周期関数ではない場合が多いです.そのような関数も周波数成分に分解して解析するために,下の図のように周期Tを伸ばして,無限の長さにします.

フーリエ変換で考える非周期関数の例.

フーリエ変換

 x(t)が周期間数のとき,上の複素フーリエ級数展開に,c_nを代入すると,

 \begin{eqnarray}x(t) &=& \sum_{n=-\infty}^{\infty} c_n \exp \frac{2 \pi i n t}{T} \\
 &=& \sum_{n=-\infty}^{\infty} \left\{\frac{1}{T} \int_{-T/2}^{T/2} \! x(t) \exp \left(- \frac{2 \pi i n t}{T} \right) \, dt \right\} \exp \frac{2 \pi i n t}{T} \\
\end{eqnarray}

となってます.時刻 t は連続なのに,周期 n/Tは飛び飛びの値 (離散)になってます (理由を考えてみてください).

 上の式で,周期Tの長さを無限にしてみます.

 周波数を,

 \displaystyle f_n = \frac{n}{T}

とすれば,この周波数は 1/T 刻みで変化しています.ということで,

 \displaystyle \Delta f = \frac{1}{T}

と置きます.上の式で, T \to \inftyをとってみます.区分求積法で,\sumが,積分になります.つまり,

 \begin{eqnarray}x(t) &=& \lim_{T \to \infty} \sum_{n=-\infty}^{\infty} \left\{\frac{1}{T} \int_{-T/2}^{T/2} \! x(t) \exp \left(- \frac{2 \pi i n t}{T} \right) \, dt \right\} \exp \frac{2 \pi i n t}{T} \\
 &=& \lim_{T \to \infty} \sum_{n=-\infty}^{\infty} \left\{\int_{-T/2}^{T/2} \! x(t) \exp \left(- 2 \pi i f_n t \right) \, dt \right\} \exp \left( 2 \pi i f_n t \right) \Delta f \\
 &=& \int_{-\infty}^{\infty} \left\{\int_{-\infty}^{\infty} \! x(t) \, e^{- 2 \pi i f t} \, dt \right\} e^{2 \pi i f t} \, df
 \end{eqnarray}

となります.ここで,上の式の中括弧の中身を X(f) と置けば,

  \begin{eqnarray}X(f) &=& \int_{-\infty}^{\infty} \! x(t) \, e^{- 2 \pi i f t} \, dt\\
x(t) &=& \int_{-\infty}^{\infty} X(f)\, e^{2 \pi i f t} \, df
 \end{eqnarray}

となります.この1番目の式は,x(t)フーリエ変換,2番目の式がフーリエ逆変換と呼ばれるものです.

【フーリエ級数】周期関数をsinとcosの和で表す

我が家にも,ついに5G (第5世代移動通信システム)の波が来ました.我が家は「おうちWi-Fi SoftBank Air」を使っています.今週,この端末のランプが「5G」の受信を教えてくれました.ということで,原稿を依頼されている「フーリエ解析」の図の準備をかねて,今回は「フーリエ級数」のお話です.ここでの目的は,説明に役立つ図を作ることで,フーリエ級数の教科書的な説明ではないので,文章の説明は期待しないでください.

周期関数とは

 フーリエ級数では,周期的に同じ形が現れるグラフを対象にします.つまり,周期を2Lとすれば,

 f(x) = f(x + 2 L)

となる関数です.例えば下のグラフのような関数です.

周期関数の例.

下の図のように,不連続でも大丈夫です.ただし,|f(x)|の面積が有限な場合だけです.

周期関数の例.不連続でも|f(x)|の面積が有限なら大丈夫.

フーリエ級数

 上の例のような周期2Lの関数f(x)は,

 \displaystyle f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos \frac{n \pi x}{L} + b_n \sin \frac{n \pi x}{L} \right)

の形に展開できます.この展開をフーリエ級数展開と言います.細かいことですが,関数f(x)に不連続な点がある場合は,フーリエ級数展開はその離れた2点の中央の点を通ります.

 係数の\{a_0, a_1, b_1, a_2, b_2, \cdots\}は,

 \begin{eqnarray}a_n &=& \frac{1}{L} \int_{-L}^{L} \! f(x) \cos \frac{n \pi x}{L} \, dx \\
b_n &=& \frac{1}{L} \int_{-L}^{L} \! f(x) \sin \frac{n \pi x}{L} \, dx
\end{eqnarray}

で求まります.

三角関数の直交性

 フーリエ級数展開の係数の\{a_0, a_1, b_1, a_2, b_2, \cdots\}を求める計算では,三角関数の直交性を使っています.これは,2つのベクトルが直交するとき内積が0になるってやつと同じことです.

 ここでは,関数f(x)g(x)内積

 \displaystyle  \int_{-L}^{L} \! f(x) \, g(x) \, dx

で定義します.f(x)g(x)の部分にサインとかコサインの関数を当てはめると,

 \begin{eqnarray} \int_{-L}^{L} \! \sin \frac{m \pi x}{L} \, \cos \frac{n \pi x}{L} \, dx &=& 0 \\
\int_{-L}^{L} \! \sin \frac{m \pi x}{L} \, \sin \frac{n \pi x}{L} \, dx &=& \left\{
\begin{array}{ll}
0 & (m \neq n) \\
L & (m = n)
\end{array}
\right. \\
\int_{-L}^{L} \! \cos \frac{m \pi x}{L} \, \cos \frac{n \pi x}{L} \, dx &=&  \left\{
\begin{array}{ll}
0 & (m \neq n) \\
L & (m = n)
\end{array}
\right. 
\end{eqnarray}

となります.これらの式で,0になるってことは,ベクトルのときみたいに,お互いが直交しているってことです.

 今回は,上の積分を計算して証明するとか,よくあるお話は飛ばして,グラフを観察してみます.

グラフで三角関数の直交性を観察

 ここでは,例として

 「 f(x) = \sin \frac{3 \pi x}{L}フーリエ級数で表せ」

という,ふざけた問題を考えます.こんな感じの問題,テストで出す先生います.学生の理解度を見るにはよい問題かもしれません.答えは,もちろん,

  f(x) = \sin \frac{3 \pi x}{L}

です.

 答えは見ただけで分かるけど,あえて,
 
 \begin{eqnarray} \int_{-L}^{L} \! \sin \frac{3 \pi x}{L} \, \cos \frac{n \pi x}{L} \, dx
\end{eqnarray}

 \begin{eqnarray} \int_{-L}^{L} \! \sin \frac{3 \pi x}{L} \, \sin \frac{n \pi x}{L} \, dx
\end{eqnarray}

について考えてみます.とはいえ,計算しません.グラフを観察します.

 n=1のとき,下の図のようになります.左側が,\sin \frac{3 \pi x}{L} \, \cos \frac{n \pi x}{L}の説明用で,右側が\sin \frac{3 \pi x}{L} \, \sin \frac{n \pi x}{L}の説明用です.これらの積は,下側の実線のグラフになるので,全面積を求めると,プラス側の面積 (赤)とマイナス側の面積 (青)がちょうど同じになり,積分が0になります.左下のグラフはそんな気がする程度ですが,右下のグラフでは,赤と青が等しいことが見て分かります.

(左上)  f(x) = \sin \frac{3 \pi x}{L} (黒実線)と  f(x) = \cos \frac{\pi x}{L} (桃破線)のグラフ.(右上)  f(x) = \sin \frac{3 \pi x}{L} (黒実線)と f(x) = \sin \frac{\pi x}{L} (水破線)のグラフ.(左下)  f(x) = \sin \frac{3 \pi x}{L} f(x) = \cos \frac{\pi x}{L}の積のグラフ.(右上)  f(x) = \sin \frac{3 \pi x}{L} f(x) = \sin \frac{\pi x}{L} の積のグラフ.

 n=2のとき,下の図のようになります.左側が,\sin \frac{3 \pi x}{L} \, \cos \frac{2 \pi x}{L}の説明用で,右側が\sin \frac{3 \pi x}{L} \, \sin \frac{2 \pi x}{L}の説明用です.これらの積は,下側の実線のグラフになるので,全面積を求めると,プラス側の面積 (赤)とマイナス側の面積 (青)がちょうど同じになり,積分が0になります.これは,下のグラフをみれば0になっていることが分かります.

(左上)  f(x) = \sin \frac{3 \pi x}{L} (黒実線)と  f(x) = \cos \frac{2 \pi x}{L} (桃破線)のグラフ.(右上)  f(x) = \sin \frac{3 \pi x}{L} (黒実線)と f(x) = \sin \frac{2 \pi x}{L} (水破線)のグラフ.(左下)  f(x) = \sin \frac{3 \pi x}{L} f(x) = \cos \frac{2 \pi x}{L}の積のグラフ.(右上)  f(x) = \sin \frac{3 \pi x}{L} f(x) = \sin \frac{2 \pi x}{L} の積のグラフ.

 n=3のとき,下の図のようになります.左側が,\sin \frac{3 \pi x}{L} \, \cos \frac{3 \pi x}{L}の説明用で,右側が\sin \frac{3 \pi x}{L} \, \sin \frac{3 \pi x}{L}の説明用です.これらの積は,下側の実線のグラフになります.左下のグラフでは,プラスの面積 (赤)しかないので,0でないことが分かります.右下は,プラスとマイナスがちょうど同じになり,積分が0になることが分かります.

n3

 n=4のとき,下の図のようになります.全面積を求めると,プラス側の面積 (赤)とマイナス側の面積 (青)がちょうど同じになり,積分が0になることが分かります.

(左上)  f(x) = \sin \frac{3 \pi x}{L} (黒実線)と  f(x) = \cos \frac{4 \pi x}{L} (桃破線)のグラフ.(右上)  f(x) = \sin \frac{3 \pi x}{L} (黒実線)と f(x) = \sin \frac{4 \pi x}{L} (水破線)のグラフ.(左下)  f(x) = \sin \frac{3 \pi x}{L} f(x) = \cos \frac{4 \pi x}{L}の積のグラフ.(右上)  f(x) = \sin \frac{3 \pi x}{L} f(x) = \sin \frac{4 \pi x}{L} の積のグラフ.

まとめ

 三角関数の直交性を,何とかアニメにしたかったのですが,よいアイデアが思いつきませんでした.ホームページの良さは、無駄に長い文章でも書くことができるし,動画を載せることもできることだと思います.なるべくアニメで,数学的な概念を説明していきたいと思います.

【受動歩行の物理】リムレスホイール

Rでリムレスホイールの動画を作ったので載せておきます.これは,解析解をグラフィックとして描いただけで,数値解は計算していません.

 受動歩行を知るためにMcGeerが登場するYoutubeの動画を見てみてください.

www.youtube.com

 リムレスホイールは,McGeerの受動歩行の論文で登場します.
https://journals.sagepub.com/doi/10.1177/027836499000900206

 リムレスホイールが歩行のモデルなの?と信じられないかもしれませんが,リムレスホイールは,受動歩行の物理を理解するのに役立ちます.

水平面の場合

 水平面でリムレスホイールを転がすと,有限回回転して止まってしまいます.

水平面を転がるリムレスホイール.時間tの単位は秒です.

斜面の場合

 斜面でリムレスホイールを転がすと,一定のパターンで斜面を下る運動に漸近します.

斜面を転がるリムレスホイール. 時間tの単位は秒です.


www.youtube.com


www.youtube.com

【受動歩行の物理の基礎】角運動量の保存

受動歩行を単純化したモデルとしてリムレスホイールがあります.

受動歩行のモデルであるリムレスホイール

今回は,リムレスホイールの運動を理解するために,ポイントとなる角運動量について解説します.

ポイントは,以下の3つです.

  • 角運動量ベクトルの大きさは位置ベクトルと運動量ベクトルが作る平行四辺形の面積で,その向きは回転軸の右ネジ方向
  • 回転していない並進運動でも角運動量は保存
  • 回転の軸が変わっても角運動量は保存

質点の角運動量

 まずは,質点の場合の角運動量の定義の確認です.

質点の角運動量

 上の図のように,運動量{\bf p}をもつ質点の原点Oまわりの角運動量{\bf L}は,質点の位置ベクトルを{\bf r}として,

 {\bf L} = {\bf r} \times {\bf p}

で与えられます.この式の\times外積を表しています.

 角運動量{\bf L}の大きさ|{\bf L}|は,{\bf r}{\bf p}を隣り合う2辺とする平行四辺形の面積になります.

角運動量の大きさは位置ベクトル{\bf r}と運動量ベクトル{\bf p}が作る平行四辺形の面積.

 平行四辺形の面積は「底辺かける高さ」なので,上の図のように運動量ベクトルと平行な2直線の距離をdとすれば,角運動量の大きさ|{\bf L}|は,

 |{\bf L}| = d \, |{\bf p}|

になります.この距離dを考えると,角運動量の大きさの計算は簡単です.

 角運動量の向きは,右ネジが進む回転方向になります.昔,私が工事現場で働きはじめたとき,ネジの構造を知らなかったので,回す方向が分からず,バカにされました.ネジは右に回すと締まるのは世界の常識です.逆ネジは特別な用途のみに使われるネジです.工事現場での仕事は非常に勉強になりました.

角運動量の向きは,右ネジが進む回転方向.

並進運動でも角運動量は保存

 回ってなくても,運動量があれば任意の点における角運動量は保存します.並進運動でも,角運動量が存在していることを忘れないでください.回転軸が固定されれば,その軸まわりの角運動量を考えてください.

 ここでは例として,下の図のように,無重力中を並進運動する円盤 (奥行きのある円柱)を考えます.この円盤は,点Aに衝突し,その後,点Aを回転軸として回転運動をはじめるとします.

無重力中を並進運動して点Aに衝突する円盤.

 上の状況では点Aにぶつかる直前まで,円盤は回転していません.回転していなくても,点Aまわりの角運動量は存在し,保存しています.位置ベクトルと運動量ベクトルを描いた下の図を見てください.

並進運動の角運動量

 角運動量の大きさは,位置ベクトルと運動量ベクトルが作る平行四辺形の面積なので,上の図の下のように面積は変わらず一定になることが分かります.

 点Aに衝突した後は,この角運動量が点Aまわりの回転運動の角運動量になります.下の図のように,衝突後の点Aまわりの角速度\omegaを求めてみてください.

並進運動から回転運動への遷移.

回転の軸が変わっても角運動量は保存

 下の図のように,回転の軸が突然変わる状況を考えてみます.

回転軸が変わっても角運動量は保存するの?

 ポイントは,回転軸が変わっても角運動量は保存するということです.私はアホなので,昔は,何で保存するのか理解できていない時期がありました.上の図の状況を考えて,角運動量を計算して,保存することを確認してみます.

 計算の流れは,円盤を微小部分に分けて,各微小部分の角運動量を求め,それを円盤全体で積分します.

 下の図のように,微小部分について,運動量ベクトルと平行な直線を回転軸を通るようにとり,運動量ベクトルと平行な2つの直線 (下図の破線)の距離を考えれば,角運動量の大きさが計算できます.

微小部分の角運動量を計算.

 円盤の質量をm,半径をaとすれば,円盤の面密度\rhoは,

 \displaystyle \rho = \frac{m}{\pi a^2}

です.

 まず,微小部分について点Oまわりの角運動量の大きさdLを計算します.角運動量は,「(2本の平行線の距離)かける(密度)かける(面積)かける(速さ)」なので,

 \displaystyle  dL = r \rho (r dr d\theta) ( r \omega) = \frac{m \omega}{\pi a^2} r^3 dr d \theta
 
 円盤全体では,円盤の半径をaとして,

 \displaystyle  |L| = \int_{r = 0}^{a} \! \int_{\theta=0}^{2 \pi} \! \frac{m \omega}{\pi a^2} r^3 dr d \theta = \frac{ma^2}{2} \omega

円盤の慣性モーメントが,

 \displaystyle \frac{ma^2}{2}

になることが分かるので,これで正しいと思います.

 次に,回転の中心を点Aに変更した瞬間を考えてみます.角運動量は,「(2本の平行線の距離)かける(密度)かける(面積)かける(速さ)」なので,変わるのは,2本の平行線の距離の部分です.上の図の右のように考えれば,2本の平行線の距離は,r + a \sin \thetaになることが分かります.微小部分について点Aまわりの角運動量の大きさdLは,

 \displaystyle  dL = (r + a \sin \theta) \rho (r dr d\theta) ( r \omega) = \frac{m \omega}{\pi a^2} r^3 dr d \theta + \frac{m \omega}{\pi a} r^2 \sin \theta dr d \theta

 最右辺の第1項は,点Oまわりの計算と同じなので,第2項の積分が0になれば,回転軸が変わっても角運動量が変わらないことになります.ということで,第2項のみを考えて,

 \displaystyle  \int_{r = 0}^{a} \! \int_{\theta=0}^{2 \pi} \! \frac{m \omega}{\pi a} r^2 \sin \theta dr d \theta = 0

になることを確かめてください.

まとめ

 今回のポイントが理解できれば,リムレスホイールの運動も理解できると思います.

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

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

一様乱数に従う乱数の和の分布.赤い破線は正規分布
指数分布に従う乱数の和の分布.赤い破線は正規分布

疑似乱数生成コマンドは頭文字 r

 Rで確率分布を扱うコマンド表記には基本ルールがあります.乱数を生成するときは,頭文字はrです.コマンドの頭文字は以下のルールになってます.

頭文字 意味
r 疑似乱数の生成
d 確率密度関数
p 確率分布関数
q クォンタイル関数 (確率分布関数の逆関数)

この頭文字の後に,分布を表す文字列が続きます.例えば,正規乱数を生成したいときは,正規分布を意味する文字列normの頭にrを付けて, rnormとなります.その他の分布については,以下を参考にしてください.

分布 Rでの表記 パラメタ
正規分布 norm 平均mean標準偏差sd
一様分布 unif 最小値min,最大値max
t分布 t 自由度df,非心度ncp
F分布 f 自由度df1df2
指数分布 exp イベント発生率rate
対数正規分布 lnorm 対数平均meanlog,対数標準偏差sdlog
ワイブル分布 weibull 位置location,スケールscale
カイ2乗分布 chisq 自由度df,非心ncp
ガンマ分布 gamma 形状shape,スケールscale
ベータ分布 beta 形状shape1shape2,非心ncp
コーシー分布 cauchy 位置location,スケールscale
ロジスティック分布 logis 位置location,スケールscale

 パッケージを使って利用可能になる分布もあります.パッケージを使いたいときは,Rで一度だけ,

install.packages("パッケージ名")

を実行して,パッケージをインストールしてください.実際に使うときは,Rを起動するたびに,

require(パッケージ名)

を実行してください.

 パッケージを使って利用可能な分布は以下のようなものがあります (私が知っているものだけです).

分布 Rでの表記 パッケージのインストール
安定分布 stable install.packages("stabledist")
ラプラス分布 (両側指数) dexp install.packages("nimble")
Gumbel分布 gumbel install.packages("gumbel")
フレシェ分布 frechet install.packages("VGAM")
正規分布 (逆ガウス) invgauss install.packages("statmod")
逆ガンマ分布 invgamma install.packages("invgamma")
逆カイ2乗分布 invchisq install.packages("invgamma")
逆指数分布 invexp install.packages("invgamma")

 乱数を生成するときは,上のRでの表記の前にrをつけてください.必要なパラメタは,

?コマンド

を実行して調べてください.

中心極限定理の数値

 乱数を生成して,中心極限定理を体験してみます.中心極限定理は,分散が有限な確率変数をたくさん足すと,正規分布で近似できるようになるというものです.

 ここでは,確率変数をXの独立なコピーをX_iとして,

 \displaystyle S_N = \frac{X_1 + X_2 + \cdots + X_N}{\sqrt{N}}

で計算される,S_Nの分布を描いてみます.

一様分布の和

 X[-1, 1]区間の一様乱数で,S_Nの分布を描いてみます.下の図のように,Nが増えると,正規分布 (赤破線)の確率密度関数のグラフと良く一致することが分かります.

 >

一様乱数に従う確率変数の和の確率密度関数.赤破線は正規分布のフィット.

 サンプルのRスクリプトです.

# 和をとる変数の数 (Nの値を増やしてみてください)
N <- 1
###################
# 繰り返し回数
N.rep <- 10000
###
x.means <- c()
for(i in 1:N.rep){ # 繰り返し
###################
# 一様乱数(-1, 1)
# (この部分の分布を替えて実験してみてください)
x <- runif(N,-1, 1)
x.means[i] <- sum(x)/sqrt(N)
} # 繰り返し終わり
###################
# xの平均値x.meansの分布
hist(x.means,probability=TRUE,breaks="FD",main=paste("Sum of",N,"variables"),border=4,col=4)
# 正規分布の当てはめ
mu.est <- mean(x.means)
sig.est <- sd(x.means)
curve(dnorm(x,mu.est,sig.est),col=2,add=TRUE,lwd=2,lty=2)

ガンマ分布の和

 次は,非対称なガンマ分布 (shape=2.5)で実験した結果です.上のRスクリプトrunif(N,-1, 1)の部分を,rgamma(N,shape=2.5)に変更しました.この例でも,下の図のように,Nが増えると,正規分布 (赤破線)の確率密度関数のグラフと良く一致することが分かります.

ガンマ分布に従う確率変数の和の確率密度関数.赤破線は正規分布のフィット.

コーシー分布の和

 コーシー分布は,超裾の厚い分布です.下の図のようにコーシー分布に従う乱数の和をとっても,正規分布に近づく傾向はみられません.理由を考えてみてください.ヒントはレヴィーの安定分布です.

コーシー分布に従う確率変数の和の確率密度関数.赤破線は正規分布のフィット.

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

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

 「確率分布」という用語はその意味が曖昧なので,ここでは,「確率密度関数」,「確率分布関数」という呼び方を使います.以下では,連続的な区間の値をとる連続確率変数を考えます.

確率密度関数 (上)と確率分布関数 (下).確率密度関数の赤い部分の面積が,確率分布関数の値に対応する.

確率密度関数

 正規分布は,\mu\sigmaをパラメタとして,

  \displaystyle f(x) = \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{(x-\mu)^2}{2 \sigma^2}}

で与えられます.これは,正規分布確率密度関数です.

 確率密度関数 \displaystyle f(x)が与えられたとき,確率を求めたいのであれば区間を考える必要があります.例えば, \displaystyle f(x)に従う確率変数Xが, a < X < b区間の値をとる確率は,

  \displaystyle {\rm Prob}\, (a < X < b) = \int_a^b \!\! f(x) \, dx

で与えられます.

確率分布関数

 確率密度関数 \displaystyle f(x)とすれば,確率分布関数 \displaystyle F(x)

  \displaystyle F(x) = \int_{-\infty}^x \!\! f(\xi) \, d\xi

となります.この両辺を微分すると,

  \displaystyle f(x) = \frac{d F(x)}{dx}

となります.

Rで確率密度関数はdと確率分布関数はp

 Rで確率密度関数を扱うコマンドの最初の文字はd,確率分布関数を扱うコマンドの最初の文字はpです.その後に分布を示す文字列が付きます.

 例えば,正規分布 (normal distribution)の場合は,normを付けます.そして,

  • 確率密度関数を表すコマンドは,d+norm=dnorm
  • 確率分布関数を表すコマンドは,p+norm=pnorm

です.

 おまけですが,確率分布関数の逆関数を表すコマンドは,q+norm=qnormになります.

 以下は,確率密度関数,確率分布関数,確率分布関数の逆関数を描くRスクリプトの例です.

mu <- 1
sig <- 2
par(mfrow=c(3,1),mar=c(5,5,2,2))
curve(dnorm(x,mu,sig),xlim=c(mu-4*sig,mu+4*sig),ylab="f(x)",main="Probability Density Function",lwd=2,col=4)
curve(pnorm(x,mu,sig),xlim=c(mu-4*sig,mu+4*sig),ylab="F(x)",main="Probability Distribution Function",lwd=2,col=4)
curve(qnorm(x,mu,sig),xlim=c(0,1),ylab=expression(paste(F^{-1},"(x)")),main="Probability Distribution Function",lwd=2,col=4)

これを実行すると,以下の図が描かれます.

Rスクリプトの実行例

最後に

 次回は,他の分布や,疑似乱数の生成について説明します.