2022-11-10

【確率変数】変数変換したら確率密度関数はどうなるか

確率

確率変数を変数変換する話です．

変数変換した後の確率密度関数
線形変換の例
指数変換の例

変数変換した後の確率密度関数

　ここでは，確率変数 $X$ が，確率密度関数 $f_{X}(x)$ に従うとします．このとき，単調関数 (単調増加か，単調減少する関数) $g(x)$ を使って，確率変数 $Y$ が

　 $Y = g(X)$

で与えられるとします． $Y$ が従う確率密度関数 $f_{Y}(y)$ を求めたいな，というのがここでの問題です．

　変数変換しても，確率は保存しないといけません．つまり，確率密度関数の対応する部分の面積が一致して，

　 $\begin{eqnarray} P\left(a < X < b \right) &=& P\left( g(a) < Y < g(b) \right) \\ \int_{x = a}^{x=b}\!f_{X}(x)\, dx &=& \int_{y=g(a)}^{y=g(b)}\!f_{Y}(y)\, dy \end{eqnarray}$

となる必要があります．

　左辺を， $y = g(x)$ と置換した積分を考えると，

　 $\begin{eqnarray} \displaystyle \int_{x = a}^{x=b}\!f_{X}(x)\, dx = \int_{y=g(a)}^{y=g(b)}\!f_{X}\left(g^{-1}(y) \right) \frac{d \left(g^{-1}(y)\right)}{dy} \, dy \end{eqnarray}$

とできそうです． $\frac{d \left(g^{-1}(y) \right)}{dy}$ は，対応する微小面積の倍率を表しています．確率は常に0以上で，マイナスになってはダメなので，プラスマイナスは無視することにすれば，

　 $\begin{eqnarray} \displaystyle \int_{x = a}^{x=b}\!f_{X}(x)\, dx = \int_{y=g(a)}^{y=g(b)}\!f_{X}\left(g(y) \right) \left| \frac{d \left(g^{-1}(y)\right)}{dy} \right| \, dy \end{eqnarray}$

となります (一般の多重積分も，向きの違いのプラスマイナスを無視するので，ヤコビアンには絶対値がついてます)．

　ということで，確率 (面積)の保存を表す

　 $\begin{eqnarray} \int_{x = a}^{x=b}\!f_{X}(x)\, dx &=& \int_{y=g(a)}^{y=g(b)}\!f_{Y}(y)\, dy \end{eqnarray}$

と対応する部分を比べると，

　 $\begin{eqnarray} f_{Y}(y) &=& f_{X}(g^{-1}(y))\, \left| \frac{d \left(g^{-1}(y) \right)}{dy} \right| \end{eqnarray}$

となることがわかります．

　あらためて書くと，

確率変数 $X$ が確率密度関数 $f_{X}(x)$ に従うとき，単調関数を使って $Y = g(X)$ と変数変換すると， $Y$ の確率密度関数は，

　 $\begin{eqnarray} f_{Y}(y) &=& f_{X}(g^{-1}(y))\, \left| \frac{d \left(g^{-1}(y) \right)}{dy} \right| \end{eqnarray}$

ということです．

線形変換の例

　 $X$ が標準正規分布 (平均0，分散1の正規分布)に従うとします．このとき， $\sigma$ ， $\mu$ を定数として，

　 $\begin{eqnarray} f_{X}(x) &=& \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}} \end{eqnarray}$

です．

　 $\sigma$ ， $\mu$ を定数として，確率変数 $Y$ が

　 $Y = \sigma X+\mu$

で与えられる場合を考えます．

　 $\displaystyle X = \frac{Y-\mu}{\sigma}$ なので， $\displaystyle g^{-1}(y) = \frac{y-\mu}{\sigma}$ です．

　したがって，確率変数 $Y$ が従う確率密度関数は，

　 $\begin{eqnarray} f_{Y}(y) &=& f_{X}(g^{-1}(y))\, \left| \frac{d \left(g^{-1}(y) \right)}{dy} \right| \\ &=& f_{X}(g^{-1}(y))\, \frac{1}{\sigma} \\ &=& \frac{1}{\sqrt{2 \pi}} e^{-\frac{\left( \frac{y-\mu}{\sigma} \right)^2}{2}} \, \frac{1}{\sigma} \\ &=& \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(y-\mu \right)^{\,2}}{2 \sigma^2}} \end{eqnarray}$

となり，正規分布の式になります．

指数変換の例

　 $X$ が正規分布に従うとします．このとき， $\sigma$ ， $\mu$ を定数として，

　 $\begin{eqnarray} f_{X}(x) &=& \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(y-\mu \right)^{\,2}}{2 \sigma^2}} \end{eqnarray}$

です．

　確率変数 $Y$ が

　 $Y = \exp X$

で与えられる場合を考えます．

　 $\displaystyle X = \log Y$ なので， $\displaystyle g^{-1}(y) = \log y$ です．

　したがって，確率変数 $Y$ が従う確率密度関数は，

　 $\begin{eqnarray} f_{Y}(y) &=& f_{X}(g^{-1}(y))\, \left| \frac{d \left(g^{-1}(y) \right)}{dy} \right| \\ &=& f_{X}(g^{-1}(y))\, \frac{1}{y} \\ &=& \frac{1}{\sqrt{2 \pi} \sigma y} e^{-\frac{\left(\log y-\mu \right)^{\,2}}{2 \sigma^2}} \end{eqnarray}$

となり，対数正規分布の式になります．

2022-11-06

【データ解析】相関係数を式で考えてみる

確率統計

相関係数の話の続きです．
chaos-kiyono.hatenablog.com

　今回は確率変数を使って考えてみます．

$Y=aX+bW$ の係数 $a$ ， $b$ と， $X$ と $Y$ の散布図の関係．

$Y=aX+bW$ の係数 $a$ ， $b$ と， $X$ と $Y$ の散布図の関係．

相関って線形な関係
相関する確率変数の生成
Rで互いに相関する正規乱数の生成

相関って線形な関係

　平均0，分散 $\sigma^2$ の確率変数 $X$ と，これと独立な平均0，分散 $\sigma^2$ の確率変数 $W$ を仮定します． $X$ と $W$ の分布については，分散が有限ということだけで形は指定しません．

　 $X$ と線形な依存性がある確率変数 $Y$ を，定数 $a$ ， $b$ を使って次の形で表します．

　 $\begin{eqnarray} Y = a X + b W\end{eqnarray}$

この場合， $a \neq 0$ のとき， $Y$ には $X$ と線形な関係があり， $W$ の係数の $|b|$ が大きいほどその関係性が薄くなります．つまり， $X$ と $Y$ の相関は， $a$ の絶対値が， $b$ の絶対値よりも大きいほど強いと言えそうです．今回は実際に相関係数を計算してこのことを足しためてみます．

以下の計算で使う基本事項として， $Y$ の期待値と分散を計算しておきます． $Y$ の期待値は，

　 $\begin{eqnarray} \mathbb{E} [Y] &=& \mathbb{E} [a X + b W ] \\ &=& a \mathbb{E} [X] + b \mathbb{E} [W ] \\ &=& a \cdot 0 + b \cdot 0 \\ &=& 0 \end{eqnarray}$

$Y$ の分散は，

　 $\begin{eqnarray} {\rm Var} [Y] &=& {\rm Var} [a X + b W ] \\ &=& {\rm Var} [a X] + {\rm Var} [b W ] \\ &=& a^2 {\rm Var} [X] + b^2 {\rm Var} [ W ] \\ &=& a^2 \sigma^2 + b^2 \sigma^2 \\ &=& \left( a^2 + b^2 \right) \sigma^2 \end{eqnarray}$

です．

　上の例では， $a=0$ なら $X$ と $Y$ は独立 ( $W$ は $X$ とは独立なので)です．また， $b=0$ なら $X$ と $Y$ は単に比例しています．

　相関係数 $r$ の定義は，

　 $r = \frac{\displaystyle \mathbb{E} \big[ \left(X- \mathbb{E} \left[ X \right] \right)\left(Y- \mathbb{E} \left[ Y \right] \right)\big]}{\displaystyle \sqrt{\mathbb{E} \big[ \left(X- \mathbb{E} \left[ X \right] \right)^2\big]} \sqrt{\mathbb{E} \big[ \left(Y- \mathbb{E} \left[ Y \right] \right)^2\big]}}$

です．

　上で定義した， $X$ と $Y$ の相関係数を計算すると，

　 $\begin{eqnarray} r &=& \frac{ \mathbb{E} \big[ X \left(a X + bW \right) \big]}{\sqrt{\mathbb{E} \big[ X^2 \big]} \sqrt{\mathbb{E} \big[ \left(a X + b W \right)^2\big]}} \\ &=& \frac{ \mathbb{E} \big[ a X^2 + b X W \big]}{\sqrt{{\rm Var} \big[ X \big]} \sqrt{{\rm Var} \big[ a X + b W \big]}} \\ &=& \frac{ a\, \sigma^2}{\sigma \sqrt{\left( a^2 + b^2 \right) \sigma^2}} \\ &=& \frac{a}{ \sqrt{a^2+b^2}} \end{eqnarray}$

になります．最後の式は，直角三角形をイメージしてもらって，直角をはさむ２辺の長さを $a$ ， $b$ としたときに，斜辺の長さ $\sqrt{a^2+b^2}$ と $a$ の比になっています．

　上の結果を使うと， $a=0$ のとき，つまり， $Y$ と $X$ に相関がないケースでは， $r=0$ になります．

　また， $b=0$ のとき，つまり， $Y$ が $X$ の定数倍のケースでは， $r=1$ か， $r=-1$ になります．

　そして，相関係数の範囲は，

　 $a > 0$ のとき， $0 < r \le 1$ ，

　 $a < 0$ のとき， $-1 \le r < 0$

になることが分かります．

　ということで，相関係数がとる値の範囲は， $-1 \le r \le 1$ です．

　また，「相関係数の推定量」

　 $\hat{r} = \frac{\displaystyle \frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)\left(y_i-\bar{y} \right)}{\displaystyle \sqrt{\frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)^2} \sqrt{\frac{1}{N}\sum_{i=1}^N \left(y_i-\bar{y} \right)^2}}$

についても， $-1 \le \hat{r} \le 1$ になります．上の計算はこの証明にはなっていないので，この範囲については自分で考えるか，インターネットで検索して調べてみてください．当然，いろんな教科書に載っています．

相関する確率変数の生成

　上の計算を簡単にすることを考えます．簡単化のため，互いに独立な確率変数 $X$ と $W$ は，平均0，分散1とします．そして， $Y$ の分散も1にします．このとき，

　　 $\begin{eqnarray} \mathbb{E} [Y^2] &=& \mathbb{E} [\left(a X + b W \right)^2 ] \\ &=& a^2 + b^2 = 1 \end{eqnarray}$

が成り立つ必要があるので， $b^2=1-a^2$ となります．

　以上の条件でもう一度，相関係数 $r$ を計算してみると，

　 $\begin{eqnarray} r &=& \frac{ a}{\sqrt{a^2 +b^2}} \\ &=& \frac{a}{\sqrt{a^2 + \left( 1 - a^2 \right)}} \\ &=& a \end{eqnarray}$

になります．結局， $a = r$ になりました．

　ということで，相関係数が $r$ の確率変数 $X$ と $Y$ を作りたいとき，まず，平均0，分散1の互いに独立な確率変数 $X$ と $W$ を用意して，

　 $Y = r X + \sqrt{1-r^2} W$

とすれば良いと言うことです．

Rで互いに相関する正規乱数の生成

　真の相関係数が $r$ になる正規乱数をRで生成してみます．使ったスクリプトはこれです．

#　データの長さ
N <- 1000
# 真の相関係数
r <- 0.5
# 2変量乱数の生成
x <- rnorm(N) #　正規乱数
w <- rnorm(N) #　正規乱数
y <- r*x+sqrt(1-r^2)*w
# 相関係数の推定 
r.est <- cor(x,y)
###
# 散布図を描画
y.max <- max(abs(c(x,y)))
par(pty="s",cex.main=1.5)
plot(x,y,xlim=c(-y.max*1.3,y.max*1.3),ylim=c(-y.max*1.3,y.max*1.3),pch=16,col="blue",main=sprintf("r (est) = %.2f",r.est))
if(r >= 0){
 text(y.max*1.2,-y.max*1.2,sprintf("r (true) = %.2f",r),pos=2,cex=1.4,col=2)
 text(-y.max*1.2,y.max*1.2,sprintf("N = %d",N),pos=4,cex=1.4,col=1)
}else{
 text(-y.max*1.2,-y.max*1.2,sprintf("r (true) = %.2f",r),pos=4,cex=1.4,col=2)
 text(y.max*1.2,y.max*1.2,sprintf("N = %d",N),pos=2,cex=1.4,col=1)
}

2022-11-04

【データ解析】相関係数の直感的理解

統計確率

相関係数の話の続きです．

chaos-kiyono.hatenablog.com

　今回は，相関係数が正になったり，負になったりすることを直感的に理解することが目的です．

相関係数の推定量

　対応する2変量の実現値 $\{(x_i, y_i)\}$ が $N$ 点あったとき，相関係数 $r$ の推定量は，

　 $r = \frac{\displaystyle \frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)\left(y_i-\bar{y} \right)}{\displaystyle \sqrt{\frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)^2} \sqrt{\frac{1}{N}\sum_{i=1}^N \left(y_i-\bar{y} \right)^2}}$

が一般的です．ここで， $\bar{x}$ ， $\bar{y}$ は，それぞれ， $\{x_i\}$ ， $\{y_i\}$ の標本平均です．

散布図の構造と相関係数

　上の式で計算した相関係数が， $1$ に近かったり， $0$ に近かったり， $-1$ に近かったりすることと，散布図の関係はどうなってるのでしょうか．

　標本平均 $\bar{x}$ ， $\bar{y}$ を引くのは，値のばらつきを0まわりにするためなので，以下では簡単化のため， $\{x_i\}$ ， $\{y_i\}$ の標本平均が0の場合を考えます．

　標本平均が0のとき，相関係数の推定量は

　 $r = \frac{\displaystyle \frac{1}{N}\sum_{i=1}^N x_i \, y_i}{\displaystyle \sqrt{\frac{1}{N}\sum_{i=1}^N x_i^2} \sqrt{\frac{1}{N}\sum_{i=1}^N y_i^2}}$

です．この分子

　 $\displaystyle \frac{1}{N}\sum_{i=1}^N x_i \, y_i$

の値を，散布図の象限毎に分けて考えると，相関係数の正になったり，負になったりすることの直感的な意味が見えてきます．

　下の図で，第1象限と第3象限にある点は，符号が同じなので，

　 $\displaystyle x_i \, y_i > 0$

です．それに対し，第2象限と第4象限にある点は，符号が逆なので，

　 $\displaystyle x_i \, y_i < 0$

です．

　下の図では，一番右の図が，象限毎の $\displaystyle x_i \, y_i$ の和を $N$ で割った値を示しています．

　散布図が，右上がりだと， $\displaystyle x_i \, y_i > 0$ の割合が多いので，相関係数は正になります．特に右上がり45°の直線に近いほど， $\displaystyle x_i \, y_i$ の値が大きくなります．

　逆に散布図が，右下がりだと， $\displaystyle x_i \, y_i < 0$ の割合が多いので，相関係数は負になります．

　相関係数が0に近いときは， $\displaystyle x_i \, y_i > 0$ になる領域と， $\displaystyle x_i \, y_i < 0$ になる領域で，ほぼ対称に点が分布しているということです．

今日のまとめ

　今日，Breaking Down 6を見ました．私は昔，番組を間違えてBreaking Down 3のペーパービューを買ってしまったことがあります．そのときは，せっかくなので何試合か見ましたが，試合の面白さがまったく分かりませんでした．その後，Breaking Downへの出場のためのオーディションをYoutubeで見るようになり，Breaking Downに出場して有名になりたいという人たちの振舞に，面白みを感じるようになりました．殴り合いの喧嘩って良いことではありませんが，外から見ているとドキドキする感じはします．

　私は昔，歌舞伎町の飲み屋で見ず知らずの方々と喧嘩になり，顔を激しく殴られました．今でも顔の左に痺れが残っています．あのとき脳にダメージを受けていなければ，もっと科学者として活躍できたかもしれませんが，今では命があっただけ良かったと思います．若いときのことを思い出すと，なんであんなに愚かだったのかと，恥ずかしくなります．

2022-11-03

【時系列解析】相関係数

確率統計

今回は，相関係数のお話です．相関というのは2つの変量間の直線的な関係性の程度を表しています．直線性とは，2つの変量の値を縦軸と横軸にとった散布図に現れる構造です．この構造のイメージを持ってほしいので，いくつか図を作りました．今日は眠いので，説明は省いて，図メインで載せていきます．

相関係数と直線性

　相関係数の定義は，ググってください．相関係数 $r$ は， $-1$ から $1$ の値をとります．データ数が1000点ぐらいあるときは，相関係数と散布図の関係は以下のようになります．

　データ数が少なく20点ぐらいのときは，相関係数と散布図の関係は以下のようになります．

データ点数とp値

　人によっては，検定のp値は小さいほど良いとか，値が0.05以下でないと意味がないと思っているかもしれません．私は「p値」信者ではないので，そんなこと思いませんが．

　相関係数を議論するときは，相関の強さと，p値の小ささは分けて解釈する必要があります．私が論文を書くときは，「相関係数が0.4以上のときに，相関があると判断する」とか，相関係数に基準を設けます．基準は，0.4でも，0.3でもいいと思います．

　以下の例では，相関係数を0.01に設定しました．つまり，ほぼ相関はありません．それでも，データ数を増やしていけば，ときどき，pは0.05よりも小さくなります．これで，「有意な相関がある」って言える？ってことです．

2022-11-01

【Rで時系列解析】計測時刻のずれを推定して修正

複数のデバイスで同時計測をしたとき，計測時刻のズレが発生することがあります．今回は，そのズレを相互相関を使って修正する方法について説明します．

　下の図では，左上のグラフのように，赤の $x_1$ と青の $x_2$ の時系列がズレています．本当は赤と青の線がほぼ重なるはずなのに，パソコンや計測機器の時刻設定の違いでズレてしまった．ということが，結構あります．

　この時間のズレは，相互相関係数の絶対値が最大になる時間差を推定することで評価できます．上の例では，右上の相互相関のグラフが最大になる点 (lag=38)が，推定された時間差になります．相関係数というのは，2つの変量の直線的な関係性を評価していますので，右下の散布図のように，相互相関が最大になる点 (lag=38)で，2つの変数の値が最も45度の直線的に集まってきます．

　時間差が推定できれば， $x_2$ の時間を推定されたラグだけずらせば，上の左下のグラフのように，赤と青の線がほぼ重なります．

数値実験のRスクリプト

　数値的に時系列を生成して，時間差を補正するスクリプトの例が以下です．

# 真のラグを設定
lag.true <- 38
# サンプル時系列を生成
N <- 1000
time <- 1:(N+lag.true)
x.base <- dnorm(time,N/2,N/5)*sin(2*pi*time/(N/3.5))^2
#
x1 <- x.base[1:N]+rnorm(N)*0.0001
x2 <- x.base[(1+lag.true):(N+lag.true)]+rnorm(N)*0.0001
time <- 1:N
######################################
par(mfrow=c(2,2),cex=1,cex.lab=1.6)
par(pty="m",mgp=c(2.2,0.5,0))
plot(time,x1,"l",col=2,xlab="i",ylab="",yaxt="n",ylim=c(-0.0005,0.003),lwd=2)
lines(time,x2,col=4,lwd=2)
legend("topright",legend=c(expression(paste(x[1],"[i]")),expression(paste(x[2],"[i]"))),col=c(2,4),lwd=c(2,2),cex=1.2)
#################
# 相互相関を使って時間を推定
xcor <- ccf(x1,x2,lag.max=round(N/4),main="",ylab="Cross-Correlation",col="skyblue")
# 相互相関係数の絶対値が最大になるラグを検出
lag.max <- xcor$lag[which.max(abs(xcor$acf))]
# 推定されたラグを破線としてプロット
abline(v=lag.max,col=6,lty=2,lwd=2)
###
plot(time,x1,"l",col=2,xlab="i",ylab="",main=paste("lag =",lag.max),yaxt="n",ylim=c(-0.0005,0.003),lwd=2)
# x2の時刻をlag.maxだけずらす
lines(time+lag.max,x2,col=4,lwd=2)
legend("topright",legend=c(expression(paste(x[1],"[i]")),expression(paste(x[2],"[i-lag]"))),col=c(2,4),lwd=c(2,2),cex=1.2)
########
# 散布図
par(pty="s",mgp=c(0.9,0,0))
if(lag.max >=0){
　plot(x1[(1+lag.max):N],x2[1:(N-lag.max)],pch=16,xaxt="n",yaxt="n",xlab=expression(paste(x[1],"[i]")),ylab=expression(paste(x[2],"[i-lag]")),main=paste("lag =",lag.max),col=6)
}else{
　plot(x1[1:(N+lag.max)],x2[(1-lag.max):N],pch=16,xaxt="n",yaxt="n",xlab=expression(paste(x[1],"[i]")),ylab=expression(paste(x[2],"[i-lag]")),main=paste("lag =",lag.max),col=6)
}

時刻を補正するRスクリプト

　時刻データを補正するサンプルスクリプトが以下です．このスクリプトではt1，t2に時刻ベクトル，x1，x2に数値データを代入しないと動きません．

# サンプリング間隔
t.samp <- 1
#####################
### 基準時系列
# 時刻
t1 <- (POSIXctの時刻)
# データ
x1 <- (時系列のベクトル)
#####################
### 時刻調整を行う時系列
# 時刻
t2 <- (POSIXctの時刻)
# データ
x2 <- (時系列のベクトル)
#####################
t.start <- max(min(t1[is.finite(x1)]),min(t2[is.finite(x2)]))
t.end <- min(max(t1[is.finite(x1)]),max(t2[is.finite(x2)]))
t.resamp <- seq(t.start,t.end,t.samp) 
###
x1.tmp <- approx(t1,x1,xout=t.resamp,method="linear")$y
x2.tmp <- approx(t2,x2,xout=t.resamp,method="linear")$y
###
par(mfrow=c(1,3),cex=1,cex.lab=1.6)
plot(t.resamp,x1.tmp,"l",col=2,xlab="Original Time",ylab="")
legend("topleft",legend=c("Reference Time Series","Target Time Series"),col=c(2,4),lwd=c(1,1),cex=0.8)
lines(t.resamp,x2.tmp,col=4)
xcor <- ccf(x1.tmp,x2.tmp,lag.max=round(length(t.resamp)/4),main="",ylab="Cross-Correlation",col="skyblue")
t.lag <- xcor$lag[which.max(abs(xcor$acf))]
abline(v=t.lag,col="blue",lwd=2,lty=2)
###
plot(t.resamp,x1.tmp,"l",col=2,xlab="Corrected Time",ylab="",main=paste("Time Lag =",Lag.sec,"sec"))
Lag.sec <- t.lag*t.samp
lines(t.resamp+Lag.sec,x2.tmp,col=4)
legend("topleft",legend=c("Reference Time Series","Target Time Series"),col=c(2,4),lwd=c(1,1),cex=0.8)

　このスクリプトを適用した例が下の図です．

まとめ

　私はウクライナの研究者との共同研究を，4年くらい前から続けています．ロシア侵攻が発生した後も，ほぼ毎週，オンラインで打ち合わせをしています．しかし，2週間前から，電力不足などの問題で，ウクライナからミーティングに接続できなくなりました．

　私は東日本大震災のときに，福島県に住んでおり，震災後はしばらく近所の小学校に避難していました．そのとき，燃料不足で暖房が使えず，とても寒かった経験があります．寒すぎて私は一睡もできませんでした．そして，そんな寒さから逃れるために家族を，福島県から避難させました．原発事故よりも，寒さの方がきつかったです．

　ウクライナでは多くの人が寒さに苦しんでいると思います．状況が少しでも改善することを祈っています．

2022-10-31

【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点しかありませんでした．安いので，今年の新着を一通り買いました．これで，寒さを楽しめそうです．

2022-10-29

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

数学時系列解析フーリエ級数

今回はフーリエ変換の話です．ここでは横軸を $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番目の式がフーリエ逆変換と呼ばれるものです．

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

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

【確率変数】変数変換したら確率密度関数はどうなるか

変数変換した後の確率密度関数

線形変換の例

指数変換の例

【データ解析】相関係数を式で考えてみる

相関って線形な関係

相関する確率変数の生成

Rで互いに相関する正規乱数の生成

【データ解析】相関係数の直感的理解

相関係数の推定量

散布図の構造と相関係数

今日のまとめ

【時系列解析】相関係数

相関係数と直線性

データ点数とp値

【Rで時系列解析】計測時刻のずれを推定して修正

数値実験のRスクリプト

時刻を補正するRスクリプト

まとめ

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

計測開始時刻を抽出

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

まとめ

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

複素フーリエ級数

フーリエ変換