2024-02-02

【道具としての数学】マクローリン展開って何？

テイラー展開数学

　マクローリン展開について簡単に説明します．1年くらいぶりの，ひさびさの更新です．古い記事を，毎日，何百人も参考にしてくれているようで，嬉しい気持ちがあるので，リクエストにこたえて記事を書きます (わかりやすいか自信がありませんが)．

マクローリン展開以前の超基礎：方程式と恒等式
見た目が違うけど，多項式と恒等式の関係になる
マクローリン展開
マクローリン展開の例
収束半径
大学の試験で見かける自明なマクローリン展開の問題

マクローリン展開以前の超基礎：方程式と恒等式

　【注意】方程式と恒等式の違いを知っている優秀な方々は，私の記事ではなく，他の方のマクローリン展開の解説を読んだほうが役に立つと思います．

　変数 $x$ についての方程式は， $x$ に隠れた値を探すための条件式です．例えば，方程式

　 $x^2 - 5 x + 6= 0$

は解くことができて，答えは， $x = 2, 3$ です．この答えが正しいかどうかは， $x$ に値を代入すれば確かめられます．例えば， $x=2$ を代入すれば，

　 $2^2 - 5\times 2 + 6= 4 - 10 + 6 = 0$

になります．

　一方で，変数 $x$ についての恒等式は， $=$ でつながれた左右が，全く同じものになっています．つまり，「自分＝自分」，「 $x = x$ 」，「 $x^2 = x^2$ 」です．左右が全く同じものなので， $x$ にどんな値を代入しても，等しくなります．例えば，恒等式， $x^2 = x^2$ に，何を代入しても正しいし，両辺の微分をしても，両辺の逆数や対数をとっても，等しいままです．なぜなら，左右が全く同じものだからです．さらに，当たり前に感じる説明をつづけると， $a, b, c, \alpha, \beta, \gamma$ を定数とするとき，

　 $a x^2 + b x + c = \alpha x^2 + \beta x + \gamma$

が成り立つのであれば， $a = \alpha$ ， $b = \beta$ ， $c = \gamma$ です．また，

　 $a \cos x + b \sin x = \alpha \cos x + \beta \sin x$

が成り立つのであれば， $a = \alpha$ ， $b = \beta$ です．

　中高校生にするような説明をしましたが，

「マクローリン展開の式は恒等式」

です．

見た目が違うけど，多項式と恒等式の関係になる

　マクローリン展開の例として， $\cos x$ のマクローリン展開を書いてみると，

　 $\displaystyle \cos x=1-\frac{1}{2 !} x^2+\frac{1}{4 !} x^4- \frac{1}{6 !} x^6+\cdots$

になります．この式が不思議なのは，恒等式は「自分＝自分」なのに，左右で見た目が全然違うことです．左辺は $\cos x$ なのに，右辺は，無限に項が続く多項式 ( $x^n$ でできた式)です．

　具体的な数値を代入して計算する場合，無限個の項を計算できないので，右辺の多項式を途中で打ち切ったりします． $x$ の値が0に近いときは，何乗もすると，高次の項はほぼ0なので，無視しても影響は小さいです．なので，マクローリン展開は，多項式を途中で打ち切って，近似式として使うことができます．

　私は数学者ではないので，数学的に厳密なことは完全に無視して説明すると，マクローリン展開というのは，「式を多項式で表しちゃうということ」です．ちょっと，かっこつけていうと，「 $x=0$ のまわりで多項式展開する」ということです．

　自分でマクローリン展開を計算したいときは，次の関係式を使います．

マクローリン展開

　関数 $f(x)$ のマクローリン展開は，

$\begin{aligned} f(x)&=f(0)+f^{\prime}(0) x+\frac{f^{\prime \prime}(0)}{2 !} x^2+\cdots \cdots+\frac{f^{(n)}(0)}{n !} x^n+\cdots \cdots \\ &=\sum_{n=0}^{\infty} \frac{f^{(n)}(0)}{n !} x^n \end{aligned}$

です． $f(x)$ は何度でも微分できる必要があります．

　何で上の関係が成り立つのかは， $f(x)$ が多項式と恒等式になることを仮定すれば導けます．つまり， $a_0, a_1, a_2, \cdots$ を定数として，

　 $\begin{aligned} f(x)&=a_0+a_1 x+a_2 x^2+\cdots \end{aligned}$

を仮定します．これは，恒等式なので，両辺の $x$ に同じ値を代入しても，両辺を微分しても等しいはずです．

　まず， $x=0$ を代入すれば，

　 $f(0) = a_0$

になるので， $a_0$ の部分には， $f(0)$ が入ることがわかります．

　つぎに， $f(x) = a_0+a_1 x+a_2 x^2+\cdots$ の両辺を微分した，

　 $\begin{aligned} \frac{d f(x)}{dx} &= a_1 + 2 a_2 x+ 3 a_3 x^2 + \cdots \end{aligned}$

に， $x=0$ を代入すれば，

　 $\begin{aligned} a_1 &= \frac{d f(0)}{dx} \end{aligned}$

であることがわかります．

　これを繰り返せば，

　 $\displaystyle a_n = \frac{f^{(n)}(0)}{n !}$

であることが導けます．

マクローリン展開の例

　ここでは，実例として，

　 $\displaystyle \frac{1}{1-x}$

をマクロ―リン展開してみます．つまり，上の公式で $\displaystyle f(x)=\frac{1}{1-x}$ の場合を考えます．

　これを，繰り返し微分してみると，

　 $\begin{aligned}\frac{d f(x)}{dx} &= \frac{d}{dx} \frac{1}{1-x} = \frac{1}{\left(1-x\right)^2}\\ \frac{d^2 f(x)}{dx^2} &= \frac{d^2}{dx^2} \frac{1}{\left(1-x\right)^2} = \frac{1}{\left(1-x\right)^3}\\ \frac{d^3 f(x)}{dx^3} &= \frac{d^3}{dx^3} \frac{1}{\left(1-x\right)^3} = \frac{1}{\left(1-x\right)^4}\\ \end{aligned}$

になるので，なんとなく，

　 $\displaystyle \frac{d^n f(x)}{dx^n} = \frac{1}{\left(1-x\right)^{n+1}}$

になりそうです (細かいことは気にしないので証明なんてしません)．これに， $x=0$ を代入すれば，

　 $\displaystyle \frac{d^n f(0)}{dx^n} = \frac{1}{\left(1-0\right)^{n+1}} = 1$

になるので，マクローリン展開の式に代入すれば，

　 $\displaystyle \frac{1}{1-x} = 1 + x + x^2 + x^3 + \cdots$

になります．実は，この式は，すべての $x$ で成り立ちません． $-1 < x < 1$ の範囲で成り立ちます．この領域を収束半径と言います． $x=1$ を代入すれば，発散するので，これがダメそうな雰囲気がすると思います．

収束半径

　マクローリン展開が，すべての $x$ で成り立たずに， $x=0$ のまわりの領域でしか正しくないときがあります．この領域を収束半径と言います．収束半径の，説明は私には面倒なので，皆さんが複素解析を勉強したときに自分で理解してください．

大学の試験で見かける自明なマクローリン展開の問題

　私は数学者ではありませんが，ある私立大学で6年間数学を教えていたことがあります．私は出題したことはありませんが，

　「 $x^2+2 x + 3$ をマクローリン展開せよ．」

みたいな問題を出す先生がいました．これは， $f(x)=x^2+2 x + 3$ として，マクローリン展開の公式に当てはめて．．．，とか考えると，時間の無駄です．

　マクローリン展開した結果は，元の式と恒等式の関係になっているので，微分とか計算する必要はありません．答えとして，次数の順に並べ替えたものを書くなら，答えは，

　 $x^2+2 x + 3 = 3 + 2 x + x^2$

です．フーリエ変換でも，似た感じの問題を出す先生がときどきいます．

2023-07-03

【Rでカオス】ロジスティック写像の分岐図

今年度は研究室のセミナーで，カオスの勉強をしています．今回は，ロジスティック写像の分岐図を描くRスクリプトを載せておきます．まずは，ロジスティック写像について簡単に説明します．

ロジスティック写像 (logistic map)

　ロジスティック写像は，以下の差分方程式で定義されます．

　 $x_{n+1} = r x_n \left(1 - x_n \right)$

ここで， $0 < x_n < 1$ ， $r$ は制御変数で， $0 < r \le 4$ の範囲で考えることが多いです．R.L. デバネーの「カオス力学系入門」では， $r > 4$ で，不安定なカオス軌道の存在を示したりしますが，通常は $r \le 4$ を考えます．

　ロジスティック写像の振る舞いについては，私が作った以下の動画を参考にしてください．検索すると，もっと詳しい説明を見つけられると思います．
youtu.be

youtu.be

Rで分岐図

以下のRスクリプトを実行すると，ロジスティック写像の分岐図を描くことができます．

###############
# 制御変数の範囲
r.range <- c(2.8, 4)
# 変化させるrの総数
n.r <- 500
###############
r.list <- seq(r.range[1],r.range[2],length=n.r)
###
# xの初期値
x0 <- 0.5
###
plot(c(),c(),xlim=r.range,ylim=c(0,1),xaxs="i",yaxs="i",xlab="r",ylab=expression({x}[n]))
###
# 過渡状態ステップ数
n.trans <- 500
# 周期の確認ステップ数
n.per <- 256
# 非周期の場合の追加ステップ数
n.plt <- 1000
for(r in r.list){
 x.tmp <- x0
 for(i in 1:n.trans){
  x.tmp <- r*x.tmp*(1-x.tmp)
 }
 xn <- c()
 xn[1] <- x.tmp
 for(i in 2:n.per){
  xn[i] <- r*xn[i-1]*(1-xn[i-1])
 }
 xn.u <- unique(xn)
 n.xn <- length(xn.u)
 if(n.xn <= 64){ # 64周期以下であることを判定
  points(rep(r,n.xn),xn.u,"p",pch=16,cex=0.05,col=rgb(0,0.5,1,alpha=1))
 }else{
  for(i in (n.per+1):(n.per+n.plt)){
   xn[i] <- r*xn[i-1]*(1-xn[i-1])
  }
  n.xn <- length(xn)
  points(rep(r,n.xn),xn,"p",pch=16,cex=0.05,col=rgb(0,0.5,1,alpha=0.05))
 }
 # 最後の値を初期値に設定
 x0 <- tail(xn,1)
}

このスクリプトを，少し変えて描いた分岐図の例が以下です．

分岐図を描く領域を変えて，分岐図に見られる自己相似性を調べてみてください．

2023-04-22

【RでWorm plot】正規性を可視化して確認

データが正規分布に従うかどうかを確認するときにQQ plotっていうのを使うことがあります (QQ plotは正規分布以外でも使えますが)．最近私は，身長，体重，BMIを分析する研究をしていています．これは，痩せすぎとか，太りすぎとか，そういった人を判定することを目的にしています．この研究に関連したWHOの報告書などで正規性を可視化するのに"Worm plot"っていうのが使われていることがあります．私もWorm plotを描きたいと思い，パッケージに頼らずにWorm plotを描いてみました．

　ということで，Worm plotを描くRスクリプトはこれです．ここでは，xにデータベクトルが入っているとします．

# xにデータを格納（ここでは例として，正規乱数を生成しています）
x <- rnorm(1000)
# データの長さ
n <- length(x)
# 標準化
x.stdzd <- scale(x)
#####################
# プロット位置の分位数を計算する関数
qt.plot <- function(k,n){
 return((k-0.5)/n)
}
#####################
# normal quantileの位置
qt <- qnorm(p=qt.plot(1:n,n))
# 小さい順に並び替え
x.stdzd.o <- sort(x.stdzd)
##############
# zスコアとnから標準誤差を計算
z.SE <- function(z, n) {
  sqrt(pnorm(z)*(1-pnorm(z))/n)/dnorm(z)
}
###
x.abs.max <- ceiling(max(abs(x.stdzd),na.rm=TRUE)/0.5)*0.5
x.range <- c(-x.abs.max,x.abs.max)
z.plot <- seq(-x.abs.max,x.abs.max,length.out = 500)
sd.x <- sd(x.stdzd)
SE <- sd.x*z.SE(z.plot, n)
###
# トレンドを除く
x.detrend <- x.stdzd.o - sd.x*qt
###
par(mar=c(5,5,2,2),las=1,cex.lab=1.5,cex.axis=1.3)
plot(qt,x.detrend,pch=16,col=4,lwd=1.5,xlab="Unit normal quantile",ylab="",xlim=x.range,xaxs="i",yaxs="i",ylim=c(-0.5,0.5))
mtext("Deviation",side=2,line=3.5,cex=1.6,las=3)
lines(z.plot,1.96*SE,col="gray10",lwd=1.5,lty=2)
lines(z.plot,-1.96*SE,col="gray10",lwd=1.5,lty=2)
abline(h=0,col="gray50",lwd=1.5,lty=2)
abline(v=0,col="gray50",lwd=1.5,lty=2)

　このスクリプトを実行すると，下のような図が描かれます．

　この図の上下の曲線が正規分布の95%信頼区間を表しています．したがって，グラフがこの信頼区間の間に収まっていれば，正規分布かもしれないということです．

　Worm plotは，QQ plotから傾きを除いて水平にし，正規分布の95%信頼区間を書いただけの図です．技術的に難しいことはありませんが，あまり普及していないので，日本語で検索しても説明が出てこないと思います．

Worm plotの読み方

　元のデータが正規分布だと，ヒストグラム法で推定した確率密度関数と，worm plotは下の図のようになります．

　正規分布と比べて左右のバランスが崩れた歪みがあるデータ

x <- -(rnorm(1000)+5)^0.2

だと下の図のようになります．

　正規分布と比べて中心が尖って，裾が厚い分布

x <- rt(1000,df=4)

だと下の図のようになります．

　今まで世界の代表的研究者が，BMIを，べき乗変換 (Box-Cox変換)すると正規分布で近似できると考えているようですが，子供についてはそうではないと思います．日本人の大規模データをみると，明らかに非正規です．近いうちに論文を発表したいと思います．

参考資料

　以下のページを参考にさせていただきました．ありがとうございました．大変勉強になりました．
www.r-bloggers.com

2023-03-09

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に設定
)

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

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

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

【道具としての数学】マクローリン展開って何？

マクローリン展開以前の超基礎：方程式と恒等式

見た目が違うけど，多項式と恒等式の関係になる

マクローリン展開

マクローリン展開の例

収束半径

大学の試験で見かける自明なマクローリン展開の問題

【Rでカオス】ロジスティック写像の分岐図

ロジスティック写像 (logistic map)

Rで分岐図

【RでWorm plot】正規性を可視化して確認

Worm plotの読み方

参考資料

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

q代数のq対数

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

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

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

まとめ

Rスクリプト

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

スペクトル表現定理

何が言いたいの？

まとめ

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