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

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

単位ステップ関数の微分 ― 直感的理解 ―

線形システム 数学

単位ステップ関数u(t)の微分について考えてみます．公式的に、

　\displaystyle \frac{du(t)}{dt} = \delta(t)

が成り立ちます．ここで， \delta(t)は，単位インパルス関数です．今回は，何でこうなるの？というのを直感的に考えていきます．今回は，どちらかといえば，線形システムの応答を考えるためのお話です．

単位ステップ関数の定義

　単位ステップ関数の定義っていろいろあります．下の図みたいな感じです．違いは，すべて，t=0の瞬間にあります．

単位ステップ関数の定義

　上の図の(a)を，よく見る気がします．数式処理ソフトMathematicaでは，UnitStepでこの定義を使っています．

　上の図の(b)の定義もときどき見かけます．私のイメージでは，静止状態から「よーいドン」と，t=0 から走り出すように考えたいので，システムの応答を考えるときは，この定義が好みです．

　上の図の(c)の定義もときどき見かけます．ここでまとめたように，t=0 の扱いはいろんな流儀があるので，t=0 のことは考えたくない気持ちが表れています．物理的には，t=0の瞬間だけを考えても意味はないので，なくても構わない気もします．

　上の図の(d)の定義は，ヘヴィサイドの階段関数です．ディラックのデルタ関数の積分と考えられるので，物理で空間的な分布を考えるときは，この定義が良いと私は考えています．

単位ステップ関数の微分

　単位ステップ関数では，t=0 で瞬間的に値が0から1にジャンプします．システムへの入力を考えるとき，そのように値が不連続に変化するのは不自然に感じてしまいます．例えば，電流とか，水流とか，力とかは，一定値の入力を実現したくても，一定値になるまで多少時間がかかります．ということで，ここでは，単位ステップ関数をより自然な形にした下の図のような連続な関数を考えます．つまり，t=0 から入力が直線的に増加し，t=\varepsilon で 1 に達すると考えます．

有限時間 \varepsilon で一定値に達する入力

　この関数の微分 (下図左下)を考えて，\varepsilon \to 0 の極限をとると，下の図右のようになります．

単位ステップ関数の微分

　入力の微分は，t=0 の直後に，面積1をもつ関数になっているので， \varepsilon \to 0 の極限をとれば，単位インパルス関数 \delta(t) になります．ここで，

　\displaystyle \int_{0}^{\infty} \delta(t) \, dt = 1

です．単位インパルス関数については，前の記事を参考にしてください．
　
chaos-kiyono.hatenablog.com

　システムの応答を考えるとき，すべては t=0 からはじまります．でも，局在する電荷などを考える物理では，時刻 t ではなく位置 x を考え， x = 0 について対称な関数を考えた方が便利です．そのような場合は，下の図のように考えます (図では，横軸が tになってます)．

ヘビサイト関数の微分

この場合は，ヘビサイト関数を微分すると，ディラックのデルタ関数になります．ここで，

　\displaystyle \int_{-\infty}^{\infty} \delta(t) \, dt = 1

であり，上の単位インパルス関数とはちょっと違います．つまり，ヘビサイト関数を微分した関数では，

　\displaystyle \int_{0}^{\infty} \delta(t) \, dt = \frac{1}{2}

です．

まとめ

　単位ステップ関数や，単位インパルス関数を使うと，理想化された状況を扱えるので，数学的には便利です．ただし，理想化する前が，どんな状況だったかを考えないと，直感とあわなくなったり，単位ステップ関数と単位インパルス関数の定義が気持ち悪く感じたり，ということになると思います．単位ステップ関数と単位インパルス関数については，自分が考えたい問題にふさわしい定義を使ってください．

METs (メッツ)の周辺 ー 安静時代謝率と代謝等量 ー

健康を維持・増進するために，運動習慣は良い効果があるということがいろいろな研究で示されてきました．どのくらい運動したのかを知るために，運動の強度やエネルギー消費量を測る必要があります．運動の強さと量を表す単位として，METS (メッツ)というのが使われることが多いです．METsは，Metabolic equivalents (代謝等量)を略したものです．

　安静時座位の酸素消費量 (安静時代謝率)を基準として，その x 倍の酸素を消費する運動が x METsになります。今回は，METsの推定に登場する安静時代謝率 (resting metabolic rate: RMR)について調べたことをメモっておきます．インターネットで「METS」を検索すると，親切な説明がたくさんあって助かるのですが，私は職業柄，根拠となる引用文献が知りたくなるので，今回はMETsの基準である安静時代謝率に関する文献を少し調べてみました．

よく使われる安静時代謝率の値

　安静時代謝率の値として，

　3.5\ {\rm mL}\, {\rm O}_2\, {\rm kg}^{-1}\, {\rm min}^{-1}

というのが良く登場します．

これは，いつ，誰が，どのような実験で推定したのでしょうか．私は調べてみましたが，答えは分かりませんでした．知っている人は教えてください．

　私が論文を検索してみたところ，1974年に出版された

Electrocardiographic responses to atrial pacing and multistage treadmill exercise testing: Correlation with coronary arteriography - ScienceDirect

や，1976年に出版された

doi.org

には，Metabolic equivalents (METs)が登場します．しかし，本文中にMETsの定義や説明はありません．2番目の論文に登場するMETsの値をみると，3.5\ {\rm mL}\, {\rm O}_2\, {\rm kg}^{-1}\, {\rm min}^{-1}くらいの値で酸素消費量を割ってMETsを計算しているように見えます．

　さらに，1972年に出版された

doi.org

に，"1 metabolic equivalent"が登場します．この場合，仰臥位安静（寝て安静）での酸素消費量が3.5\ {\rm mL}\, {\rm O}_2\, {\rm kg}^{-1}\, {\rm min}^{-1}のように書いてあります．寝てるのか，座っているのかどっちかはっきりしてほしいです．

　引用なしの記述にうんざりしていたところ，

https://doi.org/10.1152/japplphysiol.00023.2004

のイントロにやや詳しいお話を見つけました．といっても，結局よく分からない，みたいに書いてあります．なんとなくの根拠について，この論文では，体重70 kg，40歳男性1名の計測値として，3.5\ {\rm mL}\, {\rm O}_2\, {\rm kg}^{-1}\, {\rm min}^{-1}がえられたと書いてあります．現時点で，その根拠となる引用文献を入手できていないので，確証はありません．

　また，スポーツ医学でMETsの基準を提案する論文

2011 Compendium of Physical Activities: A Second Update of C... : Medicine & Science in Sports & Exercise

には，酸素消費量として，3.5\ {\rm mL}\, {\rm O}_2\, {\rm kg}^{-1}\, {\rm min}^{-1}が強く推奨されている印象がありました．根拠となる文献として，

Balke B. The effect of physical exercise on the metabolic potential, a crucial measure of physical fitness. In: Staley S, Cureton T, Huelster L, Barry AJ, editors. Exercise and Fitness. Chicago: The Athletic Institute; 1960. pp. 73–81.

が引用されていました．こんな確認しにくい本を根拠にするなよと，文句を言いたいです．しかたないので，図書館に請求して，この記事を手に入れたいと思います．

安静時代謝率は3.5より小さい傾向

　私が探した論文では，酸素消費量として，3.5\ {\rm mL}\, {\rm O}_2\, {\rm kg}^{-1}\, {\rm min}^{-1}は，大きすぎで，実際はもっと小さいという報告が多かったです．例えば，

doi.org

です．これは，2021年に出版されたシステマティックレビューなので，過去の結果をほとんど網羅していると思います．

まとめ

　私が，安静時の酸素消費量についてこだわるのは，METsを他の単位，kcalや，Wに変換したいからです．熱ストレスを評価する基準のISO7243では，代謝率の基準がWで書いてあります．今はMETsが一般的になっているので，私自身でMETsの値に直すのですが (下記参照)，安静時の酸素消費量の値が違うと，計算結果が変わってしまいます．ちなみに，熱ストレスの基準も裸の男性3人くらいの実験結果に基づいているようなので，根拠が薄いように感じます．それでも多くの人が信じているようなので，私はおとなしくしておきます．

METsとWの関係

　1 METsを，W (ワット)に変換する計算メモです．

　1 METsに相当する酸素消費量を3.5\ {\rm mL}\, {\rm O}_2\, {\rm kg}^{-1}\, {\rm min}^{-1}と仮定します．体内で酸素が使われるということは，体の中で何かが燃えているということです．

　生体内で栄養素が燃焼したとき，酸素1リットル当りの発熱量は，約5 kcal (4.7～5.0 kcal)になるそうです．引用は，孫引きですが，以下のリンクの192ページにあります．

https://www.mhlw.go.jp/bunya/shakaihosho/iryouseido01/pdf/info03k-06.pdf

したがって，3.5\ {\rm mL}の酸素の発熱量は，

　3.5 \times 10^{-3} {\rm L} \times 5 \times 10^3 {\rm cal/L} = 17.5\, {\rm cal}

です． 1\, {\rm cal} \approx 4.2 \, {\rm J} なので，これは，

　17.5\, {\rm cal} = 17.5 \times 4.2\, {\rm J} = 73.5 \, {\rm J}

です．Wの単位は，J/s なので，

　3.5\ {\rm mL}\, {\rm O}_2\, {\rm kg}^{-1}\, {\rm min}^{-1} = 73.5 \, {\rm J}\, / \, 60 \, {\rm kg}^{-1}\, {\rm s}^{-1} = 1.225 \, ({\rm J/s}) {\rm kg}^{-1} = 1.23 \, {\rm W} {\rm kg}^{-1}

となります．体重 70 kgとすると，1 METsは，

　 1.23 \, {\rm W} {\rm kg}^{-1} \times 70 \, {\rm kg} = 86 \, {\rm W}

くらいになります．

【インパルス応答】t = 0に閉じ込められた物語（その2）

数学 物理 線形システム

前の記事（以下のリンク）に続き，デルタ関数，あるいは，単位インパルス関数と呼ばれる\delta(t)を使った線形システムのインパルス応答を考えてみます．
chaos-kiyono.hatenablog.com

　今回は，\delta(t)を，以下で定義する \delta_{\varepsilon}(t) について \varepsilon \to 0 としたもので定義します．

\delta_{\varepsilon}(t)の定義．

簡単な線形システム例として，以下の微分方程式の初期値問題を考えます．

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta_{\varepsilon}(t), \quad y(0) = 0 \end{eqnarray}

前の記事との違いは，\delta(t) ではなく，\delta_{\varepsilon}(t)を入力にしているところです．

解いてみる

　積分因子を使って，

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta_{\varepsilon}(t), \quad y(0) = 0 \end{eqnarray}

を解いてみます (詳細は前回の記事参照)．

　微分方程式の両辺に積分因子 e^{a t}をかけて微分方程式を解きます．

　\begin{eqnarray} e^{a t} \frac{dy(t)}{dt} + a e^{a t} y(t) &=& e^{a t} \delta_{\varepsilon}(t) \\ \frac{d}{dt} \left(e^{a t} y(t) \right) &=& e^{a t} \delta_{\varepsilon}(t) \\ \int_{0}^{t} d \left(e^{a \tau} y(\tau) \right) &=& \int_{0}^{t} e^{a \tau} \delta_{\varepsilon}(\tau) \, d\tau \\ e^{a t} y(t) - e^{0} y(0) &=& \int_{0}^{t} e^{a \tau} \delta_{\varepsilon}(\tau) \, d\tau \\ e^{a t} y(t) &=& \int_{0}^{t} e^{a \tau} \delta_{\varepsilon}(\tau) \, d\tau \\ y(t) &=& e^{-a t} \int_{0}^{t} e^{a \tau} \delta_{\varepsilon}(\tau) \, d\tau \end{eqnarray}

ここで，右辺の積分に含まれる \delta_{\varepsilon}(t) は，

　\delta_{\varepsilon } (t) = \left\{\begin{array}{lc} 0 & (t < 0) \\ \displaystyle \frac{1}{\varepsilon} & (0 \le t \le \varepsilon) \\ 0 & (\varepsilon < t) \end{array}\right.

なので，微分方程式の解 (システムの応答)は，積分を実行することで，

　y(t) = \left\{\begin{array}{lc} 0 & (t < 0) \\ \displaystyle \frac{1}{a \varepsilon} \left(1-e^{-at} \right) & (0 \le t \le \varepsilon) \\ \displaystyle \frac{1}{a \varepsilon} \left(e^{a \varepsilon} - 1 \right) e^{-at} & (\varepsilon < t) \end{array}\right.

となります．

グラフを見てみる

　解y(t)をグラフにしたものが下の図です．

\delta_{\varepsilon}(t)に対するシステムの応答 (青実線)．ピンクは\delta_{\varepsilon}(t)．\varepsilonを0に近づけた場合の振舞の変化．

　上の図から分かるように，t = 0から，入力\delta_{\varepsilon}(t)が与えられ，その応答y(t) (青実線)が生じています．\varepsilonを0に近づけて極限をとれば，それが\delta(t)の応答になると考えられます．

解の正しさを検証

　いくつかの方法で，解の正しさを確かめてみます．

t=0 での右側微分係数

　まず，t=0 での，微分係数 \displaystyle \frac{dy(0)}{dt} について考えます．微分方程式を変形して，

　\begin{eqnarray} \frac{dy(t)}{dt} &=& - a \, y(t) + \delta_{\varepsilon}(t) \end{eqnarray}

となるので，t=0を代入すれば，

　\begin{eqnarray} \frac{dy(0)}{dt} &=& - a \, y(0) + \delta_{\varepsilon}(0) = \frac{1}{\varepsilon} \end{eqnarray}

となります．ここで注意してほしいことは，上のグラフからわかるように，t=0では，グラフはなめらかではないので微分はできません．ですので，上の \displaystyle \frac{dy(0)}{dt} は，右側微分係数 \displaystyle \frac{dy(+0)}{dt} を意味します．

　次に，解析解から右側微分係数 \displaystyle \frac{dy(+0)}{dt} を計算してみます．

　y(t) = \displaystyle \frac{1}{a \varepsilon} \left(1-e^{-at} \right) \quad (0 \le t \le \varepsilon)

を使って，微分してt=0を代入するだけです．ということで，

　\displaystyle \frac{dy(+0)}{dt} = \left. \frac{1}{\varepsilon} e^{-a t} \right|_{t = 0} = \frac{1}{\varepsilon}

がえられ，微分方程式の結果と一致することが確かめられます．

　さらに，\varepsilon \to 0で，\displaystyle \frac{dy(+0)}{dt} \to \inftyとなることが分かります．

\varepsilon \to 0 の解

　今度は，\varepsilon \to 0として，\delta(t)を入力とした微分方程式

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t), \quad y(0) = 0 \end{eqnarray}

を解く場合と，

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta_{\varepsilon}(t), \quad y(0) = 0 \end{eqnarray}

を解いてから，\varepsilon \to 0をとる場合を比べてみます．

　微分方程式

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t), \quad y(0) = 0 \end{eqnarray}

は，前回解きました．解は，

　\begin{eqnarray} y(t) =\left\{\begin{array}{lc} 0 & (t \le 0) \\ e^{-at} & (t > 0) \end{array}\right. \end{eqnarray}

です．

　次に，

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta_{\varepsilon}(t), \quad y(0) = 0 \end{eqnarray}

の解

　y(t) = \left\{\begin{array}{lc} 0 & (t < 0) \\ \displaystyle \frac{1}{a \varepsilon} \left(1-e^{-at} \right) & (0 \le t \le \varepsilon) \\ \displaystyle \frac{1}{a \varepsilon} \left(e^{a \varepsilon} - 1 \right) e^{-at} & (\varepsilon < t) \end{array}\right.

について，\varepsilon \to 0の極限をとってみます．2番目の式は潰れて無くなるので，3番目の式，

　y(t) = \displaystyle \frac{1}{a \varepsilon} \left(e^{a \varepsilon} - 1 \right) e^{-at} \quad (\varepsilon < t)

の極限を考えます．計算すると

　 \displaystyle \lim_{\varepsilon \to 0} \frac{1}{a \varepsilon} \left(e^{a \varepsilon} - 1 \right) e^{-at} = \frac{1}{a } \frac{e^{a \varepsilon} - e^{a \cdot 0}}{\varepsilon} e^{-at} = e^{-at}

となるので，

　\begin{eqnarray} y(t) =\left\{\begin{array}{lc} 0 & (t \le 0) \\ e^{-at} & (t > 0) \end{array}\right. \end{eqnarray}

と一致することが分かります．

まとめ

　今回は，\delta_{\varepsilon}(t)を考えることで，インパルスとしてぎゅーっと縮めたときの振舞を見てみました．

　今回は，t = 0での計算のすっきり感から，

　\delta_{\varepsilon } (t) = \left\{\begin{array}{lc} 0 & (t < 0) \\ \displaystyle \frac{1}{\varepsilon} & (0 \le t \le \varepsilon) \\ 0 & (\varepsilon < t) \end{array}\right.

としましたが，

　\delta_{\varepsilon } (t) = \left\{\begin{array}{lc} 0 & (t \le 0) \\ \displaystyle \frac{1}{\varepsilon} & (0 < t < \varepsilon) \\ 0 & (\varepsilon \le t) \end{array}\right.

としても，問題はないと私は思います．これは，下の図のようになります．

もう一つの\delta_{\varepsilon}(t)の定義

この場合，t=0で右側微分係数が定義できないのでダメなのか，と悩むかもしれませんが，t=0での左側微分係数を考えれば，何の矛盾もないように感じます．つまり，この場合，

　\delta_{\varepsilon } (0) = 0，さらには，\delta (0) = 0

なので，左側微分係\displaystyle \frac{dy(-0)}{dt} = 0と考えれば，微分方程式を満たしています．私としては，このように，t > 0に，\delta (t)の面積を置きたいというのが，正直な気持ちです．

　いずれにせよ，私には数学的な正しさが判断できませんので，各自で判断するか，数学に詳しい人に聞いてください．

【インパルス応答】t = 0 に閉じ込められた物語

前の記事（以下のリンク）で，デルタ関数，あるいは，単位インパルス関数がどこに立っているのかにこだわりました．今回は，線形システムのインパルス応答を考えて，私がなぜインパルスが立っている位置（時刻）にこだわるのかを説明します．
chaos-kiyono.hatenablog.com


　インパルス応答では，下の図の上段のように，t=0の瞬間にいろんなことが起こりすぎてややこしいです．私としては，まず下の図の下段のように考えて，それをt=0に閉じ込めたとみなしたいです．

インパルス応答

　以下では具体例として，

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta_{\varepsilon}(t) \end{eqnarray}

を考えます．ここで，aは正の定数，\delta_{\varepsilon}(t)は下の図のように定義します．

\delta_{\varepsilon}(t)の定義．右の図は，t>0に置いてt=0に近づける場合．

　この微分方程式でa=3として，\delta_{\varepsilon}(t)の\varepsilonを変えて，数値的に解いた結果が下の図です．ピンクの部分が\delta_{\varepsilon}(t)，青の線が応答（出力）y(t)で，\varepsilonを0に近づけています．

入力\delta_{\varepsilon}(t) (ピンク)に対するシステムの応答（数値計算結果）．\varepsilonを0に近づけた結果．

　通常は，\varepsilon \to 0とした結果を考えています．以下では，\varepsilon \to 0としたものを，\delta(t)として，
　
　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t), \quad y(0) = 0 \end{eqnarray}

を解いてみたいと思います．解き方はいくつかあります．私が思いつくのは3通りの方法です： (1) 同次方程式の解を求めてから定数変化法で特解を求める; (2) 積分因子を掛けて解く; (3) ラプラス変換で解く．

定数変化法で解く

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t) \end{eqnarray}

は，非同次（非斉次）方程式（外部からの入力を表すtの関数を含む式）です．この方程式の一般解は，同次（斉次）方程式（外部からの入力を表すtの関数がない式）

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& 0 \end{eqnarray}

の一般解y_0 (t)と，右辺に\delta(t)を含む非同次方程式の特解y_{1} (t)の和で与えられます．なぜなら，y_0 (t)とy_{1} (t)は，それぞれ，

　\left\{\begin{array}{l} \displaystyle \frac{dy_0(t)}{dt} + a \, y_0(t) &=& 0 \\ \displaystyle \frac{dy_1(t)}{dt} + a \, y_1(t) &=& \delta(t) \end{array}\right.

を満たすので，辺々足して，

　\displaystyle \frac{d}{dt} \left(y_0(t) + y_1(t) \right) + a \left(y_0(t) + y_1(t) \right) = \delta(t)

となるからです．結局，

　 y(t) = y_0(t) + y_1(t)

も元の微分方程式の解になることが分かります．

　ここから，実際に解いていきます．まず，同次方程式の一般解を求めると，

　\begin{eqnarray} \frac{dy(t)}{dt} &=& - a \, y(t) \\ \int \frac{1}{y(t)} \, dy(t) &=& - \int a \, dt \\ \log \left| y(t) \right| &=& - a t + C_0 \\ y(t) &=& \pm e^{C_0}\, e^{- a t} \\ y(t) &=& C e^{-at} \end{eqnarray}
　
となります．ここで，C_0とCは定数で，C = \pm e^{C_0}です．

　定数変化法では，上の一般解の定数部分Cを，tの関数C(t)に変えちゃいます．つまり，y(t) = C(t) e^{-at}を非同次方程式に代入し，C(t)を求めます．

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta (t) \\ e^{-at} \frac{dC(t)}{dt} - a C(t) e^{-at} + a C(t) e^{-at} &=& \delta(t) \\ \frac{dC(t)}{dt} &=& \delta(t) e^{at} \\ \int_{0}^{t} dC(\tau) &=& \int_{0}^{t} \delta(\tau) \, e^{a \tau} \, d\tau \\ C(t) - C(0) &=& \left\{\begin{array}{lc} 0 & (t \le 0) \\ 1 & (t > 0) \end{array}\right. \end{eqnarray}

私の流儀では，\delta(t)の面積を，t>0の領域に置いているので，

　\begin{eqnarray}\int_{0}^{\infty} \! \delta(t) \, dt &=& 1 \\ \int_{0}^{\infty} \! f(t) \, \delta(t) \, dt &=& f(0) \end{eqnarray}

になります．さらに，任意の\varepsilon > 0に対して，

　\begin{eqnarray}\int_{0}^{\varepsilon} \! \delta(t) \, dt &=& 1 \\ \int_{0}^{\varepsilon} \! f(t) \, \delta(t) \, dt &=& f(0) \end{eqnarray}

です．

　上のC(0)は，一般解に対応するので無視すれば，特解が

　y(t) = C(t) \, e^{-at} = e^{-at} \quad (t > 0)

と求まります．したがって，一般解は，Cを定数として，

　\begin{eqnarray} y(t) =\left\{\begin{array}{lc} C e^{at} & (t \le 0) \\ C e^{at} + e^{-at} & (t > 0) \end{array}\right. \end{eqnarray}

です．初期条件 y(0)=0を満たすようにCを設定すれば (C=0)，

　\begin{eqnarray} y(t) =\left\{\begin{array}{lc} 0 & (t \le 0) \\ e^{-at} & (t > 0) \end{array}\right. \end{eqnarray}

となります．インパルス \delta(t)に対する応答は，t > 0で生じています．t = 0にインパルスを与えたと考えると，変な感じがします．ですので，私は t > 0で，tが0に近い時刻でインパルスを入れたとみなしたいのです．

積分因子を使って解く

　\begin{eqnarray} \frac{dy(t)}{dt} + p(t) \, y(t) &=& q(t) \end{eqnarray}

の形の微分方程式は，両辺に

　\displaystyle \exp \left( \int p(t)\, dt \right)

をかけることで，解くことができます．ここで，上の不定積分の定数項は無視してください．これは，左辺を積の微分の形にするアイデアです．

　具体的に，

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t) \end{eqnarray}

を解いてみます．この場合， p(t) = aです．

　\displaystyle \exp \left( \int a \, dt \right) = e^{at}

　\displaystyle e^{at}を微分方程式の両辺にかけると，

　\begin{eqnarray} e^{at} \frac{dy(t)}{dt} + a e^{at} \, y(t) &=& \delta(t) \, e^{at} \\ \frac{d}{dt} \left( e^{at} \, y(t) \right) &=& \delta(t) \, e^{at} \end{eqnarray}

となります．左辺が，積の微分の形になっていることに気づいてください．これを解くと，

　\begin{eqnarray} \frac{d}{dt} \left( e^{at} \, y(t) \right) &=& \delta(t) \, e^{at} \\ \int_0^{t} d \left( e^{a \tau} \, y(\tau) \right) &=& \int_0^{t} \delta(\tau) \, e^{a \tau} \, dt \\ e^{a t} \, y(t) - y(0) &=& \left\{\begin{array}{lc} 0 & (t \le 0) \\ 1 & (t > 0) \end{array}\right. \\ y(t) &=& \left\{\begin{array}{lc} 0 & (t \le 0) \\ e^{-at} & (t > 0) \end{array}\right. \end{eqnarray}

となります．

ラプラス変換で解く

　ここでは，ラプラス変換については詳しく説明しません．ラプラス変換と，その逆変換の公式の導出は，自分で調べてください．

　今回の微分方程式

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t), \quad y(0) = 0 \end{eqnarray}

をラプラス変換で解いてみます．

　以下では，y(t) のラプラス変換を Y(s) とします．微分方程式の両辺をラプラス変換すると，

　\begin{eqnarray} s Y(s) - y(0) + a \, Y(s) &=& 1 \end{eqnarray}

になります．初期値 y(0)=0を使うと，

　\begin{eqnarray} s Y(s) + a \, Y(s) &=& 1 \\ Y(s) &=& \frac{1}{s+a} \end{eqnarray}

となります．両辺を逆変換すると，

　\begin{eqnarray} y(t) &=& e^{-at}\, u(t) \end{eqnarray}

がえられます．ここで，u(t)は単位ステップ関数

　\begin{eqnarray} u(t)= \begin{cases} 0 & (t \le 0) \\ 1 & (t > 0)\end{cases} \end{eqnarray}

です．

　単位ステップ関数の別の定義として，

　\begin{eqnarray} u(t)= \begin{cases} 0 & (t < 0) \\ 1 & (t \ge 0)\end{cases} \end{eqnarray}

が使われることがあります．この定義だと， y(t) = e^{-at}\, u(t)とすると，y(0) = 1になって，初期条件 y(0) = 0を満たさなくなってしまします．明らかにダメです．

まとめ

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t), \quad y(0) = 0 \end{eqnarray}

を例として解いてみました．どの解き方でも解は一致します．ここで考えたかったことは，t = 0の一瞬をどうとらえるかです．私としては，t = 0の瞬間にすべてが起きたと考えるのではなく，t=0の直後の無限に短い時間の間に応答が生じたと考えたいということです．

デルタ関数，単位インパルス関数は，どこに立っている？

数学 物理

ディラック (Dirac)のデルタ関数や，単位インパルス関数（単に，インパルス関数と呼ぶことも）は，分野によって立っている位置が違うので，一つの定義にこだわっていると，計算結果に納得できないことがあります．「デルタ関数と単位インパルス関数は同じモノだよ」と明確に書かれている教科書がありますが，分野によって定義が微妙に違い，その違いを理解することが重要というのが今回のお話です．
　
　以下では，デルタ関数，あるいは，単位インパルス関数を\delta(t)とします．\delta(t)の基本性質

　\begin{eqnarray}\int_{-\infty}^{\infty} \! \delta(t) \, dt &=& 1 \\ \int_{-\infty}^{\infty} \! f(t) \, \delta(t) \, dt &=& f(0) \\ \int_{-\infty}^{\infty} \! f(t) \, \delta(t-a) \, dt &=& f(a) \end{eqnarray}

を満たすことは，以下のどの場合も共通です．

私が最初に出会ったDiracのデルタ関数

　私は物理学科出身なので，デルタ関数を以下の図を使って理解していました．

デルタ関数の定義

　つまり，上の図のように，時刻t=0について軸対称な長方形を考えます．幅\varepsilon ，高さ\displaystyle \frac{1}{\varepsilon }とするので，面積は1です．この長方形について，\varepsilon \to 0の極限をとったものがデルタ関数\delta(t)です．

　この考えた方だと，

　\begin{eqnarray} \int_{0}^{\infty} \! \delta(t)\, dt &=& \frac{1}{2} \\ \int_{0}^{\infty} \! f(t)\, \delta(t)\, dt &=& \frac{1}{2}\, f(0) \end{eqnarray}

になりそうです．実際，いくつかの教科書では，このように書かれています．

　たまに，

　\begin{eqnarray} \delta(t)=\left\{\begin{array}{cc} 0 &(t \neq 0) \\ \infty &(t=0) \end{array}\right. \end{eqnarray}

などと，強調して書いてあることがあります．しかし，\delta(t)は，積分の中に登場して初めて意味をもつのであって，こう書いても何の意味もないし，何の役にも立たないと思います．このことは，どこかの本で読んだのか，どこかの講義で聴いたのか忘れましたが，私は経験からもそう感じるようになりました．\delta(t)は，線ではなく，幅と有限な面積をもつものです．

工学系の教科書で出会った単位インパルス関数

　私が，\delta(t)の定義に違和感を覚えたのは，線形システムのインパルス応答を計算したときです．つまり，\delta(t)のラプラス変換（片側ラプラス変換）の計算です．

　\delta(t)のラプラス変換では，

　\displaystyle \int_{0}^{\infty} \! \delta(t) \, e^{-st}\, dt

を計算します．上の図をイメージすると，t \ge 0では，デルタ関数が半分しか含まれないので，結果に\frac{1}{2}が付くよね？と考えてしまいます．

　しかし，一般的な結果は，

　\displaystyle \int_{0}^{\infty} \! \delta(t) \, e^{-st}\, dt = 1

です．

　何で？と，私はしばらく悩んだ経験があります．私は\frac{1}{2}が正しいと思い込み，さんざん世間に不平不満を言ったあげくに，\delta(t)は，上とは違う図 (t=0で軸対称ではない)で定義されているということを知りました．つまり，下の図です．

単位インパルスの定義

　この図では，t=0の右側 (t \ge 0)に長方形を置き，左端を固定して右から押しつぶしていいます．「デルタ関数と単位インパルス関数は同じモノ」とか，言わないでほしかったです．個人ごとの都合に合わせて，\delta(t)の定義は変わるよ，と最初から注意しておいてほしかったです．

　結局，この定義だと，

　\begin{eqnarray}\int_{0}^{\infty} \! \delta(t) \, dt &=& 1 \\ \int_{0}^{\infty} \! f(t) \, \delta(t) \, dt &=& f(0) \end{eqnarray}

ということです．つまり，物理よく使う\delta(t)とは違うものです．

単位インパルス関数は t > 0に立っててほしい (私の願望)

　上の単位インパルス関数\delta(t)の定義を使っても，線形システムのインパルス応答を考えるときに，私には気持ち悪さが残っています．この気持ち悪さを消すために，私としては，下の図の右のように\delta(t)を，t > 0に置いてほしいです (このように\delta(t)を定義している教科書も存在します)．

\delta(t)を，t = 0 (左)じゃなくて，t > 0 (右)に立てほしい．

　
　\delta(t)を，t > 0に置くというのは，下の図で言えば，(b)か(c)のパターンです．どの場合も，\varepsilon \to 0としたものが，\delta(t)です．

\delta(t)はどこに立っているのか？

　(b)か(c)のパターンで\delta(t)を定義すると，\delta(0)の値はどうなるの？と言われそうです．上の図の(b)，(c)では，\delta(0) = 0です．おいおい，\delta(0)=\inftyだろうと，口撃されそうです．上の図の(b)，(c)を使った定義では，\delta(t) = \inftyとなるtの値を明示できません．

　最初の方でも言いましたが，\delta(0)や\delta(t)の値なんて考えても，何の意味もないし，何の役にも立たないと思います．私は数学者ではないので詳しいことは知りませんが，\delta(t)は数学的には関数じゃないそうです (面積をもつ不思議な線)．関数っぽく書くから，\delta(0)の値を覚えて満足する人が出てくるのだと思います．

　\delta(t)が，t > 0に立っていて，面積をもっているのであれば，明確に

　\begin{eqnarray}\int_{0}^{\infty} \! \delta(t) \, dt &=& 1 \\ \int_{0}^{\infty} \! f(t) \, \delta(t) \, dt &=& f(0) \end{eqnarray}

です．単純バカの私としては，この形で話をしてほしいです．

まとめ

　\delta(t)は積分の中でこそ輝くし，\delta(t)の面積が，どこに立っているのかを意識して区別しないと直感的に理解できないというのが，私の考えです．最後に注意しておきますが，今回の説明はあくまで私の個人的感想ですので，数学的な正しさは保証できません．

【Rで最小2乗フィット】QR分解 (Householder変換)を使った最小2乗法

北川先生の本「Rによる 時系列モデリング入門」で，Householder変換を使った最小2乗フィットの話が登場したので，Rでやってみます．

　ここでは，以下の式を使ってデータを生成します．

　\displaystyle y_n = \sum_{i=1}^{2} a_i x_{ni} + \varepsilon_n

ここで，\displaystyle \varepsilon_nは，平均0，標準偏差sd.epsの正規乱数とします．

　以下のスクリプトでは，データの長さをN，パラメタa_1をa1，a_2をa2に設定します．

　\vec{y} = \begin{bmatrix} y_1 \\ y_{2} \\ \vdots\\ y_{N} \\ \end{bmatrix}, \quad Z = \begin{bmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \\ \vdots & \vdots \\ x_{N1} & x_{N2} \\ \end{bmatrix}, \quad \vec{a} = \begin{bmatrix} a_1 \\ a_2 \\ \end{bmatrix} , \quad \vec{\varepsilon} = \begin{bmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \vdots \\ \varepsilon_{N} \\ \end{bmatrix}

と表せば，

　\vec{y} = Z \vec{a} + \vec{\varepsilon}

です．ベクトルはすべて縦書きの列ベクトルとします．

逆行列を使った最小2乗法

　残差としての\{\varepsilon_{n}\} の2乗和が最小になるように \vec{a } を計算すると，

　\vec{a } = \left({}^t Z Z \right)^{-1} {}^t Z \, \vec{y}

となります．{}^t Zは，Zの転置行列を表します．説明は，過去の記事を参考にしてください．
chaos-kiyono.hatenablog.com

　以下は，N = 10，a1 = 2.4，a2 = -3.0とした例です．Rのコマンドを使って逆行列を計算して，パラメタを推定しています．

# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
##############################
# a.hat = (t(Z) Z)^{-1} t(Z) Y 
##############################
a.hat <- solve(t(Z) %*% Z) %*% (t(Z) %*% Y)
# 残差の２乗平均
sig2.hat <- sum((Y - Z %*% a.hat)^2)/nrow(Y)
# a1, a2の推定結果
a.hat[,1]
# 残差の2乗平均 sig2
sig2.hat

　このスクリプトを実行すると，

> # a1, a2の推定結果
> a.hat[,1]
[1]  2.425303 -3.634048
> # 残差の2乗平均 sig2
> sig2.hat
[1] 1.607804

のように，パラメタと残差が計算されます．

QR分解 (Householder変換)を使った最小2乗法

　次は，QR分解 (Householder変換)でパラメタを推定してみます．pracmaパッケージを使いました．最初に，このパッケージをインストールしてください．

　使ったRスクリプトは以下です．

#install.packages("pracma")
require(pracma)
# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
###
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# Householder変換
UX <- householder(X)$R
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat

　このスクリプトを実行すると，

> # a1, a2の推定結果
> a.hat
[1]  2.611575 -2.950162
> # 残差の2乗平均 sig2
> sig2.hat
[1] 0.6768868

のように，パラメタと残差が計算されます．

　pracmaパッケージを使わなくても，Rのコマンドqrを使えば，QR分解できます．中身のアルゴリズムについては，何を使っているのか分かりません．Rスクリプトは以下です．

# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
###
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# QR分解でRを求める
UX <- qr.R(qr(X))
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat

まとめ

　QR分解について，勉強してみてください．QR分解の方法には，グラム・シュミットの正規直交化法とか，Householder変換を使う方法とか，ギブンス回転を使う方法などがあるそうです．

【付録】3変数の場合

　逆行列を使う場合．

N <- 10
sig2 <- 1
eps <- rnorm(N,0,sqrt(sig2))
##############################
a1 <- 2.4
a2 <- -3.0
a3 <- -1.8
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + a3 * x3[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
x3 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+a3*x3+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N],x3[1:N]),　ncol=3)
##############################
# a.hat = (t(Z) Z)^{-1} t(Z) Y 
##############################
a.hat <- solve(t(Z) %*% Z) %*% (t(Z) %*% Y)
# 残差の２乗平均
sig2.hat <- sum((Y - Z %*% a.hat)^2)/nrow(Y)
# a1, a2の推定結果
a.hat[,1]
# 残差の2乗平均 sig2
sig2.hat

　Householder変換を使う場合．

#install.packages("pracma")
require(pracma)
###############
N <- 10
sig2 <- 1
eps <- rnorm(N,0,sqrt(sig2))
##############################
a1 <- 2.4
a2 <- -3.0
a3 <- -1.8
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + a3 * x3[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
x3 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+a3*x3+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N],x3[1:N]),　ncol=3)
##############################
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# Householder変換
UX <- householder(X)$R
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat

ロジスティック写像のカオス

今回は，決定論的カオス (deterministic chaos)の説明用の動画を作りました．決定論的カオス (以下では，カオス)というのは，時間発展を決定論的に決めるルールがあるのに，未来が予測できなくなる現象のことです．

カオスとは

　カオスを示す決定論的システムの具体例として，ロジスティック写像

　x_{n+1} = C x_n (1-x_n)

を考えます．ここで， 0 < x_n < 1で，Cは制御変数です (0 < C \le 4)．Cの値を変えると， x_nの漸近的な振舞が変化します．

　C = 4の場合が，カオスの例になっています．カオスの主な特徴は，

  • 短期予測可能
  • 長期予測不可能
  • 非周期パターン

です．下の図では，

　下の図は，x_{n+1} = 4 x_n (1-x_n)で，初期値x_1 = 0.3 (青)の場合と x_1 = 0.300000000000001 (赤)の場合を比較した結果です．初期値の差は，10^{-15}です．

ロジスティック写像のカオス．初期値x_1 = 0.3 (青)と x_1 = 0.300000000000001 (赤)の比較．

　図の左は，クモの巣図法と呼ばれる時間発展をグラフを使って調べる方法のアニメです．図の右は， x_nの時間発展のグラフです．

　 n= 40くらいまでは，青と赤の2つの点がほとんど同じ振舞を示します (2点は重なっています)．これが，短期予測可能ということです．初期値の値が近ければ，短期的には似た振舞になります．

　しかし， n= 40をこえると，青と赤の2つの点が全く違った振舞を示すようになります．これが，長期予測不可能ということです．カオスでは，初期値鋭敏性といって，ほんのわずかな差が，指数関数的に大きくなるため，どんなに小さい差でも，いつかはシステムサイズくらいの差に成長してしまいます．

周期倍分岐

　ロジスティック写像で制御変数Cの値を変えると，カオス以外にも，いろいろと面白い現象が見られます．例えば，同じ値が繰り返すパターンです．Cを0から次第に増加させると，1つの値に収束したり (1周期)，2つの値を繰り返したり (2周期)，4つの値を繰り返したり (4周期)と，周期が倍々になる周期倍分岐というのが現れます．下の図は，いくつかのCでの振舞です．

C=2.8，初期値x_1 = 0.3の場合，1周期に漸近．
C=3.2，初期値x_1 = 0.3の場合，2周期に漸近．
C=3.54，初期値x_1 = 0.3の場合，4周期に漸近．

無相関と独立性の違い

　x_{n+1} = 4 x_n (1-x_n)の振舞は，無相関と独立性の違いを考える上で良い例になっています．つまり，この写像で生成される時系列\{x_n\}は，無相関な時系列になっています．でも，\{x_n\}は，短期予測可能なので，異なる時刻の値同士は独立ではありません．相関と独立の違いを理解しておいてください．

まとめ

　私は，博士 (理学)をとるまでは，カオス力学系の研究をしていました．その後，生理学とか，医学への応用を指向しました．カオスについての読み物としては，ジェイムズ・グリック著の「カオス―新しい科学をつくる」 (新潮文庫，1991)をお勧めします．ただし，カオスは，一時，ブームになり私もその影響を受けましたが，現在は，グリックの本で書かれているような勢いはないです．私は，博士過程を修了した後，運良く日本学術振興会の特別研究員（ＰＤ）に採択されました．そのころから，カオスとか，自分の狭い専門性に頼った研究をするのではなく，世の中の役に立つことを実現するために，新しいことに積極的に挑戦するようになりました．