統計的な相関まとめ
平均
揺らぎ
分散(標準偏差の二乗 or 揺らぎの二乗の平均)
統計学における分散と不偏分散 例題でわかりやすく解説 | 全人類がわかる統計学
不偏分散(アンサンブル平均ではない平均で取る)
統計学における分散と不偏分散 例題でわかりやすく解説 | 全人類がわかる統計学
共分散(分散の一般化)
共分散 | 統計用語集 | 統計WEB
相関係数
26-3. 相関係数 | 統計学の時間 | 統計WEB
スペクトル解析のコヒーレンスの形ほぼそのままであり、コヒーレンスが統計量に強く結びついているのがわかる。
koideforest.hatenadiary.com
偏相関係数(疑似相関が疑われる場合)
26-4. 偏相関係数 | 統計学の時間 | 統計WEB
因子の条件を除いた時にとの間で本当に相関があるのかをチェック出来る(らしい)。
相関係数とは何か。その求め方・公式・使い方と3つの注意点 | アタリマエ!
これらはあくまで直線的な相関を調べることしか出来ない(二次以上の依存関係は極値の場所に依って解析不能)。
また特定の分布を仮定するため、異常なデータに弱い。
スペクトル解析のコヒーレンスとアンサンブル平均
以前に相関関数をいろいろまとめた。
koideforest.hatenadiary.com
koideforest.hatenadiary.com
koideforest.hatenadiary.com
しかし気になったのは、「今のコヒーレンスの定義では 1 にしかならないのではないか?」という点である。
広島大の講義ノートにクロススペクトルの意味についての記述を見つけた。
https://home.hiroshima-u.ac.jp/nishino/Classforundergraduate/Keisoku/keisoku2shou.pdf
しかし、フーリエ変換しているのになぜか変換前の引数に依存しているし、百歩譲ってそれは良しとしても、平面波から勝手に関数が出て来ているし、意味不明である。
ポイントは、アンサンブル平均にある。
何回か同じ測定をして、データのセットが個あるとすると、
簡単のため、今として各スペクトルのアンサンブル平均を求めると、
これらを用いてコヒーレンスを計算すれば、
これで分母と分子に違いが出た。クロススペクトルによって、データセット間のもしくは同士の関係を考慮することが出来る。
したがって、このコヒーレンスはそもそも統計量に関するものであり、「入力と出力の波形がどれくらい似ているか」ではなく、「複数回測定しても強く再現される信号はどの波数成分か」を表すものと言える。(その意味で「相関」と言えなくもないが、相関関数の相関とは別な意味だと思う)
そのため、全てのデータセットが全く同じものの場合(もしくは一度しか測定を行わなかった場合)、コヒーレンスは常に 1 になる。
三角関数の相関関数
以前の記事で、相関関数、及びそのガウス関数の時の振舞について調べた。
koideforest.hatenadiary.com
koideforest.hatenadiary.com
ここでは三角関数()を使って相関関数の振舞を調べる。
フーリエ変換の定義
三角関数( )
自己相関関数の定義
スペクトルと自己相関関数
デルタ関数の二乗は以下のサイトを参照。
やっと「δ関数の2乗」:T_NAKAの阿房ブログ
無限空間の積分なので、自己相関は発散します。
相互相関関数の定義
クロススペクトルと相互相関関数
コヒーレンスと位相
やはりズラした程度ではコヒーレンスは破れず、一方で位相は並行移動の情報を含んでいる。
追記:
コヒーレンスについて再考
koideforest.hatenadiary.com
ガウス関数の相関関数
以前の記事で相関関数をまとめた。
koideforest.hatenadiary.com
ここではガウス関数を使って、振舞を調べる。
フーリエ変換の定義
自己相関関数
スペクトルと自己相関関数
で、がちゃんと求まっていることがわかる。
相互相関関数の定義
より一般的なガウス関数
クロススペクトルと相互相関関数
コヒーレンスと位相
関数を並行移動させた程度では、コヒーレンスは(最大値のまま)変わらないため、本質的な関数の形の議論に役立ちそう。
以下で追加考察。
koideforest.hatenadiary.com
一方で、位相は関数の並行移動をもろに受ける。波数に対して位相が線形ならば、関数が単に並行移動していると言えそう。
ガウス関数のフーリエ変換 (scipy.fftpack)
ガウス関数はフーリエ変換しても、数式上は(幅は変わるが)またガウス関数になる。
高速フーリエ変換(FFT)とフーリエ変換との関係を以前まとめたため、FFTを用いてガウス関数をフーリエ変換する。
koideforest.hatenadiary.com
ここではscipy.fftpackを用いてFFTを行うが、(少なくともこの場合)出力される配列データの並びが「(正の振動数のデータ)、(負の振動数のデータ)」となっており、そのままプロットするとガウス関数に見えない。
試しにやってみる。
import numpy as np import matplotlib.pyplot as plt from scipy.fftpack import fft a = 1 f = lambda x: np.exp( - a * x**2 ) # Gauss function max_x = 10. min_x = -10. N = 2**7 x = np.linspace( min_x, max_x, N ) data = f( x ) delta_x = x[1] - x[0] # ( max_x - min_x ) / ( N - 1 ) # Fast Fourier Transformation tilde_f = fft( data ) * delta_x delta_k = 1. / delta_x / ( N - 1 ) k = np.arange( 0, 1./delta_x + delta_k, delta_k ) * 2. * np.pi plt.plot( k, abs( tilde_f ), "o" ) plt.show()
データの並びに関しては、以下のサイトが参考になる。
fft - memoring
SciPy FFTpack
Fourier Transforms (scipy.fftpack) — SciPy v1.1.0 Reference Guide
もう少し説明すると、FFTにおけるアルゴリズムの性質上、fft で生成した配列は、直感的な[ -N/2, ..., N/2 ]という順番では無く、
[ 0, 1, ,2, ... (N/2 - 1), N/2, -N/2, (-N/2 + 1), ..., -1 ]
という風になっている。
つまり、
[ (正の振動数:昇順)、(負の振動:降順)]
である。
fftpackに内蔵されている fftfreq を使うと、fft と同様な並びの横軸を作成してくれる。
from scipy.fftpack import fftfreq freq = fftfreq( N, delta_x ) * 2. * np.pi plt.plot( freq, abs( tilde_f ), "-o" ) plt.show()
点でプロットするならば問題無いが、線を追加すると N/2 → -N/2 に移るときの線が現れてしまう。
そこで、配列の並びを直感的に戻してくれる関数 fftshift が用意されている。
from scipy.fftpack import fftshift tilde_f_shift = fftshift( tilde_f ) freq_shift = fftshift( freq ) gauss_k = lambda k: np.sqrt( np.pi / a ) * np.exp( - k**2 / ( 4. * a ) ) plt.plot( freq_shift, abs( tilde_f_shift ), "-o" ) put.plot( freq_shift, gauss_k( freq_shift ), "x" ) plt.show()
解析的に求めた解も一緒にプロットした。
解析解の導出は以下のサイトが分かり易い。ただし、指数の係数は省略されている。
ガウス関数のフーリエ変換 | ポップラーン
しかし、FFTを用いたガウス関数のフーリエ変換で最も問題なのは、実は数値解が振動している点である。
虚部は実部に比べて小さいため、実部のみプロットすると、
plt.plot( freq_shift, tilde_f_shift.real, "-o" ) plt.plot( freq_shift, abs( tilde_f_shift ), "-" ) plt.show()
何故、振動が現れるのかは、今の所パッと思い付かないが、フーリエ変換はなかなか侮れないことを教えてくれていると思う。
高速フーリエ変換(FFT)を使ってフーリエ変換する方法
高速フーリエ変換(FFT)は結局は離散フーリエ変換であり、どちらかというとフーリエ級数展開に近い。
高速フーリエ変換 - Wikipedia
FFTは結局は級数展開、つまり「和」なので、フーリエ変換、つまり「積分」にするためには、変換が必要である。
http://hooktail.sub.jp/mathInPhys/riemannIntegra/
多くの場合には、フーリエ変換は信号解析で用いられるため、時間と振動数を変数に用いることが多い。
ここでは信号データを、そのフーリエ成分をとして考える。
データ数がの時のFFTの式:
少なくともpythonのscipy.fftpackではこの定義に従っていると思われるため、横軸は振動数で出力される。
これをフーリエ変換を見比べる。以下ではを使い、に変数変換した後のフーリエ成分をとすると。
「和」を「積分」に変換するには、積分における「」に対応するものが必要である。それがなければ、いくらデータ間隔を細かくしたところで積分にはならない。
データ間隔を以下のように定義する。
これを用いて、和を積分に変換する。
したがって、十分大きいを用いたとき、近似的に以下の関係が成り立つ。
数式上では、を用いることが(少なくとも自分の分野では)多いため、かに揃えておかないと混乱する。
古典的な相関関数のまとめ
「相関とは内積の拡張のようなものである」という説明で納得した。
相関関数 [物理のかぎしっぽ]
離散的な数直線上の点における関数値をベクトルと思えば、内積は
と表せる。これがゼロのとき、関数は互いに直交していると言える。
ここで、を連続量に置き換え、それに伴って和を積分に置き換えれば、内積を拡張したものが得られる。
これを相関と呼ぶ。
もしくはのときをそれぞれ自己相関、相互相関と呼ぶ。
自己相関は、ベクトルで言えば「ノルムの二乗」に相当する。
自己相関がゼロで無いときに、相互相関がゼロになる場合、互いに直交すると言える。
フーリエ変換も、ある波数の平面波との相関を取ったものである。
しかし、自己相関関数と言うと、意味が変わる。
http://www.jspf.or.jp/Journal/PDF_JSPF/jspf2009_09/jspf2009_09-620.pdf
自己相関関数は、自分自身をズラしたものとの相関であり、ズラす量を引数に持つ。
を、つまり、自己相関で規格化したものを自己相関係数と呼ぶ。
関数が周期的な場合、自己相関関数(係数)も周期的に振動し、が最初にを跨ぐは変動の速さを特徴付ける量として使うことが出来る。
自己相関関数の被積分関数をフーリエ変換を用いて表すと、
(ただし、は実関数で積分区間外でゼロ())
つまり、真面目に関数をちょっとずつズラして積分しなくても、フーリエ変換の(強度の)逆フーリエ変換から自己相関関数を求めることが出来る。
今まで、「ふーん」としか思っていなかったが、よく考えると、
であり、フーリエ成分を(絶対値の)二乗して逆フーリエしたら全然違う量になっているのは、なかなか面白いのではないか?
この「逆フーリエ変換すると自己相関関数を与える」関数をスペクトルと呼ぶ。
特に、は実関数の全強度に対応し、から、「スペクトルは(周)波数ごとの強度」と解釈出来ることがわかる。また、計算したスペクトルが正しいかどうか、直接求めた全強度と比較することでチェックすることも出来る。
余談だが、平均値 からのズレであるを扱う方が便利なこともある。
自己相関が分散に対応するように自己相関関数を以下のように定義すると、
この場合、が全分散を与えるため、スペクトルは(周)波数ごとの分散を与えるものとして解釈出来る。これは揺らぎの解析に利用出来ると考えられる。
自己相関関数があるならば、相互相関関数も定義出来る。
はクロススペクトルと呼ばれる。
クロススペクトルは一般に複素数であるから、その位相は次のように求められる。
種々の相関関数から相互相関係数を作ったように、種々のスペクトルからコヒーレンスが定義出来る。
これらが実際の関数を考えた時にどうなっているのかは、今後少しずつ遊んでみたいと思う。