Cox-Box変換とq対数ー足し算とかけ算のあいだー　（執筆中）

統計で使うBox-Cox変換というのがあります．データの値を $x$ として，Box-Cox変換は

　 $x^{(\lambda)} = \begin{cases}\frac{x^\lambda-1}{\lambda} & (\lambda \neq 0) \\ \log x & (\lambda=0)\end{cases}$

で定義されます．ここで， $\lambda$ はパラメタで，この値をうまく調整してやると，正規分布に従わないデータ $\{x_i\}$ を，正規分布に従うデータ $\{x^{(\lambda)}_i\}$ に変換できちゃうことがあります．

　何でこの変換でうまくいくの？という疑問に答えている説明を，統計の教科書で見たことはありませんが，一つの説明をここではします．ポイントは，

　自然界には，「足し算」と「かけ算」の間にある演算が存在するから

ということです．例えるなら，「足し算」と「かけ算」を異なる2点としたとき，その2点を結ぶ直線上にある任意の点に対応する演算が存在するということです．

　皆さんは，中心極限定理をご存じでしょうか．中心極限定理は，「ランダムな値をとる変数をたくさん足すと，足された値は正規分布に従うようになる」というやつです．

　また，世の中には対数正規分布に従う変数もたくさんあります．対数正規分布に従う変数は，対数変換すると，正規分布に従うようになります．対数正規分布が登場する理由は，「ランダムな値をとる変数をたくさんかけてある」ということです．

　対数変換は，Box-Cox変換では $\lambda=0$ の場合に対応します． $\lambda\neq0$ の場合はどうなの？というのが問題ですが，これは，かけ算と足し算をつないでその間を埋める演算が，自然界に存在すれば，説明できます．以下では，かけ算でも，足し算でもなく，その間にある演算を「ヘンテコ算」と呼ぶことにします．

　対数正規分布の場合，たくさんのかけ算を，足し算に変換する演算が，対数変換でした．

　一般に，たくさんのヘンテコ算（かけ算と足し算の間）を，足し算に変換する演算が，q対数変換です．

q代数のq対数

2023-02-25

【Rで時系列解析】弱定常過程のスペクトル表現定理

弱定常な確率過程のスペクトル表現定理が意味することを理解するために数値実験をしてみました．今回は，前回の記事
chaos-kiyono.hatenablog.com
の続きです．

　弱定常過程のパワースペクトルについて理解してほしい部分は，「パワースペクトルは，確率過程を周波数成分に分解したときの振幅を，どのような意味で指定しているのか」ということです．パワースペクトルは，周波数成分の振幅の情報をもっていますが，これは，振幅の２乗の期待値の意味であり，振幅を決定論的に決めているわけではありません．

１次自己回帰過程のパワースペクトル

　例として１次自己回帰過程を考えます．

　 $y_{n}=a y_{n-1} + w_n$

ここで， $a$ は， $-1 < a < 1$ の定数， $w_n$ は，平均0，分散 $\sigma^2$ の白色ノイズです．このとき，パワースペクトル $S(f)$ は，

　 $\displaystyle S(f)=\frac{\sigma^2}{1+a^2-2 a \cos (2 \pi f)}=\frac{\sigma^2}{(1-a)^2-2 a( \cos (2 \pi f)-1)}$

です．

　自己相関至上主義 (英: autocorrelationism)では，パワースペクトル $S(f)$ (と平均値)が分かれば，この過程のことはすべて理解できるということになります．では，このパワースペクトルは，弱定常な確率過程の何を指定しているのでしょうか？このことを前回の記事で説明したつもりですが，伝わった気がしないし，説明が不足している部分があると思いますので，数値実験して，説明を補足します．

１次自己回帰過程を周波数成分に分解

　前回の記事で説明したスペクトル表現定理からわかることは，平均0の弱定常な離散確率過程の時系列

　 $\left\{x[0], x[1], x[2], \cdots, x\left[\left(N-1\right) \right] \right\}$

を，

　 $\displaystyle x[n] = \sum_{j = 1}^{M} \left[A_{j} \cos \left(2 \pi f_j n \right) + B_{j} \sin \left(2 \pi f_j n \right) \right]$

と分解したとき ( $f_j = j/N$ ， $M = \lfloor N/2 \rfloor$ )，

　 $\begin{eqnarray} {\rm E} \left(A_{j} \right) &=& 0\\ {\rm E} \left(B_{j} \right) &=& 0\end{eqnarray}$

です．さらに， $S \left( f_j \right)$ をこの過程の片側パワースペクトルとすると，

　 $\begin{eqnarray} {\rm E} \left(A_{j} \, A_{k} \right) &=& \left\{\begin{array}{lc} \displaystyle \frac{S \left( f_j \right)}{N} & \left( j = k \right) \\ \displaystyle 0 & \left( j \neq k \right) \end{array}\right. \\ {\rm E} \left(B_{j} \, B_{k} \right) &=& \left\{\begin{array}{lc} \displaystyle \frac{S \left( f_j \right)}{N} & \left( j = k \right) \\ \displaystyle 0 & \left( j \neq k \right) \end{array}\right. \\ {\rm E} \left(A_{j} \, B_{k} \right) &=& 0 \quad \left({\rm all} \ j \ {\rm and} \ k \right)\end{eqnarray}$

です． $S \left( f_j \right)$ は，時系列から推定したパワースペクトルではなく，真のパワースペクトルということに注意してください．

　つまり，弱定常過程を

　 $\displaystyle x[n] = \sum_{j = 1}^{M} \left[A_{j} \cos \left(2 \pi f_j n \right) + B_{j} \sin \left(2 \pi f_j n \right) \right]$

と分解したとき， $A_{j}$ と $B_{j}$ は平均0，分散 $\displaystyle \frac{S \left( f_j \right)}{N}$ の正規乱数ということです．前回， $A_{j}$ と $B_{j}$ が正規分布に従う確率変数とは説明していませんでした．

　このことを数値実験で確かめてみました．以下に掲載したRスクリプトで，N.samps <- 10000とした結果が下の図です．

　上の図の一番左は，真のパワースペクトル (赤破線)と，サンプルの平均として推定したパワースペクトル (青実線)です．各段は，パワースペクトル中に縦線で示した周波数成分の結果を示してます．左から２番目の図は，縦軸に係数 $A_{j}$ ，横軸に係数 $B_{j}$ をプロットしたものです．たくさん点があるのは，10000個のサンプルの結果をプロットしたからです．図中の赤破線の円は，パワースペクトルが表す振幅を描いたものです．右側の２つの図は， $A_{j}$ と $B_{j}$ の分布を描いたものです．赤破線は，理論的に予言される正規分布です．

　上の図の左から２番目の図から分かるように， $A_{j}$ と $B_{j}$ で与えられる振幅は，赤破線の円のまわりにあるのではなく，原点を中心としてばらついているということです．右側の２つの図から，0が中心であり，正規分布にしたがってばらつくということが確認できます．

まとめ

　パワースペクトルは，パソコンを使って簡単に計算できるので，意味を理解せずに，とりあえずパワースペクトルを計算してみる人が多いように感じます．私自身の理解は，まだまだと感じていますので，これからも基礎を大事にして勉強したいと思います．

Rスクリプト

　この記事で使用したRスクリプトは以下です．

# 1次自己回帰過程 (first-order autoregressive process) 
# y[n] = a y[n-1] + w[n]
# w[n]: white noise with zero mean and variance sig^2
####################
# パラメタ (parameter)
a <- 0.9
# w[n]の標準偏差sig^2
sig2 <- 1
# 時系列の長さ
N <- 256
####
y <- c()
####################
# 平均をとるサンプル時系列数
N.samps <- 100
####################
# 最大ラグの決定
m <- N-1
####################
# 周波数の定義
freq <- (0:(N-1))/N
####################
for(ii in 1:N.samps){
 ####################
 # 白色ノイズの生成 
 w <- rnorm(N,0,sqrt(sig2))
 ####################
 # 空の変数を事前に準備
 # 初期値をy[1]=w[1]とする
 y[1] <- w[1]
 for(n in 2:N){
  y[n] <- a*y[n-1]+w[n]
 }
 ######################
 PSD <- abs(fft(y))^2/N
 Re.Y <- Re(fft(y))/sqrt(N)
 Im.Y <- Im(fft(y))/sqrt(N)
 ######################
 if(ii == 1){
  PSD.samps <- data.frame(PSD)
  Re.Y.samps <- Re.Y
  Im.Y.samps <- Im.Y
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
  Re.Y.samps <- cbind(Re.Y.samps,Re.Y)
  Im.Y.samps <- cbind(Im.Y.samps,Im.Y)
 }
}
f <- freq[freq>=0 & freq<=0.5]
Spec <- rowMeans(PSD.samps)[freq>=0 & freq<=0.5]
#################
N.plot <- 4
n.f <- length(f)
f.samp <- round(exp(seq(log(2),log(n.f-1),length=N.plot)))
#################
# 図を描く
par(mfrow=c(N.plot,4),cex.lab=1.5,las=1,cex.axis=1.3)
###
for(j in 1:N.plot){
par(pty="m")
#################
plot(f[-1],Spec[-1],"l",col=4,log="xy",main="Power spectrum density",xlim=c(f[2]/1.1,tail(f,1)*1.1),lwd=4,xlab="f",xaxt="n",yaxt="n",xaxs="i",ylab="")
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
abline(v=f[f.samp[j]],lty=2,col="lightpink",lwd=2)
########
# 解析的に計算したパワースペクトル
curve(sig2/(1-2*a*cos(2*pi*x)+a^2),add=TRUE,col=2,lwd=2,lty=2,n=1000,xlim=c(1/N,1/2))
#
points(f[f.samp[j]],sig2/(1-2*a*cos(2*pi*f[f.samp[j]])+a^2),pch=16,cex=2,col=2)
#####
par(pty="s")
r <- sqrt(sig2/(1-2*a*cos(2*pi*f[f.samp[j]])+a^2))/sqrt(2)
#colm <- densCols(Re.Y.samps[f.samp[j],],Im.Y.samps[f.samp[j],], colramp = colorRampPalette(cm.colors(16)[8:1]))
plot(Re.Y.samps[f.samp[j],],Im.Y.samps[f.samp[j],],pch=16,col=rgb(43/255,137/255,206/255,alpha=0.2),xlim=c(-r*3.5,r*3.5),ylim=c(-r*3.5,r*3.5),xlab="Re(Y) / sqrt(N)",ylab="Im(Y) / sqrt(N)",main=paste(f[f.samp[j]],"Hz"))
curve(sqrt(r^2-x^2),xlim=c(-r,r),add=TRUE,col=2,lty=2,lwd=2)
curve(-sqrt(r^2-x^2),xlim=c(-r,r),add=TRUE,col=2,lty=2,lwd=2)
par(pty="m")
r <- sqrt(sig2/(1-2*a*cos(2*pi*f[f.samp[j]])+a^2))/sqrt(2)
hist(Re.Y.samps[f.samp[j],],xlab="Re(Y) / sqrt(N)",breaks="FD",col=4,border=4,freq=FALSE,ylab="",main=paste(f[f.samp[j]],"Hz"))
curve(dnorm(x,0,r),add=TRUE,col=2,lty=2,lwd=2)
hist(Im.Y.samps[f.samp[j],],xlab="Re(Y) / sqrt(N)",breaks="FD",col=4,border=4,freq=FALSE,ylab="",main=paste(f[f.samp[j]],"Hz"))
curve(dnorm(x,0,r),add=TRUE,col=2,lty=2,lwd=2)
}

2023-02-22

弱定常過程のパワースペクトルと表現定理

なぜ，実験でえられた時系列を分析するとき，とりあえず，パワースペクトルを計算するのか？

　それは，自己相関至上主義を信仰している人が多いからです．自己相関至上主義って何？と思った皆さん，ググっても，何も出てこないと思います．ちなみにChatGPTで「自己相関至上主義とは」と入力すると，

　自己相関至上主義（じこそうかんしじょうしゅぎ、英: autocorrelationism）とは、ある時系列データが自己相関を持っていることを前提とし、その自己相関関係に基づいて分析・予測を行おうとする考え方・手法のことを指します。

　自己相関は、あるデータの過去の値と現在の値との相関関係を示す指標です。自己相関を持っている時系列データは、過去のデータから未来のデータを予測することが可能であると考えられます。そのため、自己相関至上主義者は、時系列データの分析や予測において、自己相関を重視する傾向があります。

　一方で、自己相関至上主義にはいくつかの欠点が指摘されています。例えば、自己相関関係があることを前提とした分析手法では、時系列データの性質によっては適切な結果を得られないことがある、という点が挙げられます。また、時系列データの背後にある複雑な要因を無視することにもなり得るため、解釈には注意が必要です。

と，私が想定していた意味に近い内容を答えてくれます．学生の答案であれば高得点をあげます．しかし，弱定常に関連した説明がないので，満点ではありません．注意しておきますが，自己相関至上主義とか，autocorrelationismとか，そんな用語は一般に使われていません．つまり，私が作った架空の用語です．

　自己相関至上主義という皮肉った表現を使って私が言いたいことは，弱定常な確率過程を意識して，時系列が評価されることが多いということです．弱定常では，平均値と自己共分散 (自己相関)が時間シフトに対して不変です．時間発展のルールは「自己相関」のみで完全に指定されます．これは数学的な扱いを簡単にするための仮定で，現実世界の時系列がすべて弱定常というわけではありません．

　最初の疑問に戻って，「なぜ，パワースペクトル？」という疑問に答えると，自己相関至上主義では，自己相関が最も大事なもので，自己相関とパワースペクトルには1対1の関係があり，どちらも同じ情報をもっているからです．自己相関至上主義は，パワースペクトル至上主義と言い換えることもできます (こんなこと言っている人いませんが)．観測時系列から自己相関関数の推定する場合，時間ラグが大きい領域で推定精度が急激に悪くなったり，非定常なトレンドの影響を大きく受けたり，という推定における欠点があります．それらの点については，パワースペクトルの推定の方がましですので，自己相関関数ではなく，パワースペクトルの推定が使われるとうことです．

　導入が長くなりましたが，今回は，弱定常過程のスペクトル表現定理についてのお話です．証明なんかは一切しません．スペクトル表現定理は，パワースペクトルを読み解き，確率過程の特性を解釈するときの基礎となります．パワースペクトルの解釈では，波に分解して読み解く，といった説明がされることがあると思います．私も，時間がないときは，その点だけ説明することがあります．しかし，この説明では，波形のパターンが固定された振動 (確定信号)を想定している印象があります．したがって，確率過程を考えるときには，その視点だけでは不十分です．

スペクトル表現定理

　ここでは，スペクトル表現定理を一般的に説明するというよりは，スペクトル表現定理から言える具体的な内容を説明します．

以下では，平均0の弱定常な離散確率過程の時系列を

　 $\left\{x[0], x[1], x[2], \cdots, x\left[\left(N-1\right) \right] \right\}$

とします (値は実数です)．

　一般に，

　 $\displaystyle x[n] = \sum_{j = 1}^{M} \left[A_{j} \cos \left(2 \pi f_j n \right) + B_{j} \sin \left(2 \pi f_j n \right) \right]$

と書けます．ここで， $f_j = j/N$ ， $M = \lfloor N/2 \rfloor$ です．

　スペクトル表現定理から，弱定常過程では， $A_{j}$ と $B_{j}$ は平均0の確率変数になります．つまり，

　 $\begin{eqnarray} {\rm E} \left(A_{j} \right) &=& 0\\ {\rm E} \left(B_{j} \right) &=& 0\end{eqnarray}$

です．

　さらに， $A_{j}$ と $B_{j}$ は

　 $\begin{eqnarray} {\rm E} \left(A_{j} \, A_{k} \right) &=& \left\{\begin{array}{lc} \displaystyle \sigma_j^2 & \left( j = k \right) \\ \displaystyle 0 & \left( j \neq k \right) \end{array}\right. \\ {\rm E} \left(B_{j} \, B_{k} \right) &=& \left\{\begin{array}{lc} \displaystyle \sigma_j^2 & \left( j = k \right) \\ \displaystyle 0 & \left( j \neq k \right) \end{array}\right. \\ {\rm E} \left(A_{j} \, B_{k} \right) &=& 0 \quad \left({\rm all} \ j \ {\rm and} \ k \right)\end{eqnarray}$

を満たします．

　これらの性質から

　 $\begin{eqnarray} {\rm E} \left(x[n]^2 \right) &=& \sum_{j = 1}^{M} \sigma_j^2 \\ {\rm E} \left(x[n]\, x[(n-k)] \right) &=& \sum_{j = 1}^{M} \sigma_j^2 \, \cos \left(2 \pi f_j k \right) \end{eqnarray}$

を示すことができます．

　上で現れる $\sigma_j^2$ は，この過程のパワースペクトル (両側パワースペクトル)を $S \left( f_j \right)$ とすれば，

　 $\displaystyle \sigma_j^2 = 2 \frac{S \left( f_j \right)}{N}$

で与えられます．ただし， $f_j = 1/2$ のときは，

　 $\displaystyle \sigma_j^2 = \frac{S \left( f_j \right)}{N}$

です．片側パワースペクトルであれば，後者の関係がどの周波数でも成り立ちます．片側パワースペクトルとは，パワースペクトルの負の周波数側を正の周波数側に足したものです．

　 $S \left( f_j \right)$ を片側パワースペクトルとして，上の式に代入してみれば，

　 $\begin{eqnarray} {\rm E} \left(x[n]^2 \right) &=& \sum_{j = 1}^{M} \sigma_j^2 = \sum_{j = 1}^{M} \frac{S \left( f_j \right)}{N} = \sum_{j = 1}^{M} S \left( f_j \right) \, \Delta f \end{eqnarray}$

は，パーセバルの定理，

　 $\begin{eqnarray} {\rm E} \left(x[n]\, x[(n-k)] \right) &=& \sum_{j = 1}^{M} \sigma_j^2 \, \cos \left(2 \pi f_j k \right) = \frac{1}{N} \sum_{j = 1}^{M} S \left( f_j \right) \, \cos \left(2 \pi f_j k \right) \end{eqnarray}$

は，ウィーナー・ヒンチンの定理に対応します．

　スペクトル表現定理では， $N \to \infty$ とした，より一般的な場合が扱われます．その話は，今回は面倒なので省略します．

何が言いたいの？

　ここで私が言いたいことは，パワースペクトルが，弱定常過程の特性をどのように指定しているのか，ということです．

下の図は，1次自己回帰過程のパワースペクトルを描いたものです．図中の赤い破線が，理論的に計算した真のパワースペクトルです．3つの図があるのは，パワースペクトルの平均をとるサンプル数を変えるとどうなるかを見るためです．一番左は1本，中は10本，右は100本の時系列のパワースペクトルを推定して，平均をとっています．

1次自己回帰過程のパワースペクトル．赤破線が理論的結果．青は数値的に生成した時系列からパワースペクトルを推定した結果．

　時系列が1本の場合 (上図左)は，真のパワースペクトル (赤破線)と時系列から推定したもの (青実線)とのズレが結構あります．縦軸が対数スケールなので，1000倍以上のズレがある部分もあります．しかも，時系列から推定したパワースペクトルは無茶苦茶ギザギザです．これで，真のパワースペクトル (赤破線)と合っていると言えるの？と思うかもしれません．でも，スペクトル表現定理の主張を踏まえれば，そのようになるのは当然だということが分かります．

　なぜなら，時系列を周波数成分に分解したとき，

　 $\displaystyle x[n] = \sum_{j = 1}^{M} \left[A_{j} \cos \left(2 \pi f_j n \right) + B_{j} \sin \left(2 \pi f_j n \right) \right]$

$A_{j}$ と $B_{j}$ は平均0の確率変数になり，その2乗の期待値を考えたときに，真のパワースペクトルと結びつけられます．つまり， $S \left( f_j \right)$ を片側パワースペクトルとして，

となります．隣り合う周波数でも， $A_{j}$ と $A_{j+1}$ は相関していないので，当然，時系列から推定したパワースペクトルは，滑らかにならず，ギザギザになります．

　私が言いたいことは，時系列から推定したパワースペクトルは，確率的に変動する関数だということです．

まとめ

　弱定常過程のパワースペクトルを理解するためには，スペクトル表現定理を理解することが必要だと思います．スペクトル表現定理のことをわかりやすく書いている教科書は少ないですが，

Hamilton, James Douglas. Time series analysis. Princeton university press, 2020.

に少しだけ説明があるので，興味があれば読んでみてください．私の研究室には，この本を3冊買ってあります．

2023-01-24

【Rでレーダーチャート】fmsbパッケージの利用

最近，レーダーチャートを描画する必要が生じたので，Rでレーダーチャートを描く方法を調べてみました．今回は，fmsbパッケージを使ってレーダーチャートを描くためのメモです．

これまで，fmsbパッケージを使ったことがない場合は，まず，このパッケージをインストールしてください．

install.packages("fmsb")

このパッケージを使えば，以下のような図を作れます．

使ったスクリプトは，これです．

# fmsbパッケージを使います
# install.packages("fmsb")
require(fmsb)
##################################
# タイトル
TITLE <- "TITLE"
# 変数名
LABELS <- c("A","B","C","D","E")
# 最大と最小の設定
MAX <- c(5,5,5,5,5)
MIN <- c(0,0,0,0,0)
# レーダーの目盛数
N.SEG <- 5
# 軸ラベル
AXIS.LABEL <- c("0","1","2","3","4","5")
#　得点
SCORE <- c(3,2,5,4,0)
###################
# 右まわりに並び替え
CW <- c(1,seq(length(LABELS),2,-1)) 
####################
# レーダーチャートの描画
MAX.MIN <- data.frame(matrix(c(MAX,MIN),byrow=TRUE,nrow=2))
SCORE <- data.frame(matrix(SCORE,nrow=1))
DAT.CW <- rbind(MAX.MIN,SCORE)[CW]
LABELS.CW <- LABELS[CW]
radarchart(DAT.CW, # データ（最大最小，値）のdata.frame
    vlabels=LABELS.CW, # 各値のラベル
    title=TITLE, # 上に書くタイトル
    cglcol=8, # 補助線の色 (8: 灰色)
    pcol=2, # 線の色
    pfcol=rgb(1,0,0,alpha=0.1), # 塗りつぶし 
    plwd=3, # 線の太さ
    seg=N.SEG, # 補助線の分割数
    pty=32, # マーカ (32: なし)
    axistype=1, # 軸に値を描画する場合１
    caxislabels=AXIS.LABEL, # 軸のラベル
    axislabcol=1, # 軸のラベルの色 (1: 黒)
    centerzero=TRUE # 中心を0に設定
)

私の趣味で，右まわりの順にスコアをプロットとするようにしてあります．

2023-01-08

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

線形システム数学

単位ステップ関数 $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 \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}$

です．

まとめ

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

2022-12-30

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}$

くらいになります．

2022-12-17

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

数学物理線形システム

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

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

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

　 $\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}$

の解

について， $\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.$

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

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

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

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

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

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

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

Cox-Box変換とq対数ー足し算とかけ算のあいだー　（執筆中）

q代数のq対数

【Rで時系列解析】弱定常過程のスペクトル表現定理

１次自己回帰過程のパワースペクトル

１次自己回帰過程を周波数成分に分解

まとめ

Rスクリプト

弱定常過程のパワースペクトルと表現定理

スペクトル表現定理

何が言いたいの？

まとめ

【Rでレーダーチャート】fmsbパッケージの利用

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

単位ステップ関数の定義

単位ステップ関数の微分

まとめ

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

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

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

まとめ

METsとWの関係

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

解いてみる

グラフを見てみる

解の正しさを検証

t=0 での右側微分係数

$\varepsilon \to 0$ の解

まとめ