今回もメモがてらです。わかっている人にとっては、ほぼほぼどうでもいい内容ですので、あまり期待しないでください。
ＡＲ過程のスペクトラル・デンスィティー
さて、表題の通り、今回は（も）あまり内容のない記事です。実は、最近「時系列解析」的なものに凝っていて、自分でAR時系列などを作って遊んだりしています。ただ、それだけだと、タダのノイズを眺めているだけの詰まらない遊びになってしまうので、ちょっと格好をつけて（パワー）スペクトラル・デンシティなるものを計測してみました（計測、といってもシミュレーションしたデータをFFTしただけです）。AR時系列を作ること自体は、まったくもって何の難しさもないので、コードは晒しません。とりあえず、PSDは下記のように計算しました。
def psd(y, fs=1, onesided=False, scale=None):
T = len(y)
yf = fft.fftshift(fft.fft(y))
spec = (np.abs(yf)**2)/(T*fs)
f = fft.fftshift(fft.fftfreq(T, 1/fs))
spec = spec[T//2:] if onesided else spec
spec = np.log(spec) if scale == 'log' else spec
f = f[T//2:] if onesided else f
return f, spec
一方で、AR過程には理論的なスペクトルというものが存在します。ただのノイズではないので、周波数に特性が生じるのですが、それはAR係数と呼ばれるパラメータから推計できてしまいます。計算方法は、下記のサイトから「拝借」しました。
一応、導出なども載っているので、数式とお友達になれる人は見てみましょう。例によって、にょろにょろ及びでんでんむしはありません。シグマとexpだけです。恐れるに足らず。それにしても、脳波形（EEG)とかの解析をしている人たちは凄いですねー。ほぼほぼノイズの中から、カオスだの周期だのをちゃんと見つけるというね。しかも、「ストレスと自律神経の科学」なんて、私のためにあるような学問です（患者としてお世話になるという意味で）。
class AR:
def __init__(self, a, omega):
F, G, Q = MatrixModel.AR(a, omega)
self.F = F
self.G = G
self.Q = Q
self.omega = omega
# Algorithm introduced in the document below.
# https://www3.nd.edu/~esims1/arp_companion.pdf
def response(self, lags):
h, lam = [1.0], self.F
for i in range(lags):
h.append(lam[0][0])
lam = lam @ self.F
return np.array(h)
def psd(self, n, fs=1.0, onesided=False, scale=None):
def peak(f):
p = self.F.shape[0]
num = self.omega**2
exps = np.exp(-2*np.pi*1j*np.linspace(1,p,p)*f)
denom = np.abs(1 - np.sum(self.F[0]*exps)) ** 2
return num / denom
freqz = fftshift(fftfreq(n, 1.0))
f = fs*freqz
spec = np.array([peak(f)/fs for f in freqz])
spec = np.log(spec) if scale == 'log' else spec
f = f[n//2:] if onesided else f
spec = spec[n//2:] if onesided else spec
return f, spec
※このARというクラスは状態空間モデルに組み込むことを想定しているのでわかりにくいですが、self.F[0]というやつが、ar係数をベクトルとしたものに相当します。念のため補足。ちなみにF自体は「状態空間モデル、ARモデル」などと検索すれば出てくる代物です。状態のシフトと、AR過程の最新値を行列演算だけで計算できるような仕組みになっています。行列って便利ですね。
細かい説明は割愛しますが、ＰＳＤに関しては、まぁ、数式をそのまま実装しただけです。まだ、「しゅうはすうりょういき」という考え方に慣れない私は、ふとこの式に渡すfってなんやねん、と思ってしまうのですが、まぁ、周波数ですね。周波数（まだ、この辺がうすぼんやり）。一般に、教科書に載っているような数式は、離散型のＦＦＴを想定しているわけではないので、サンプリングレートなどという概念については触れられていません。たぶん、あっていると思うのですが、その分の補正は私が勝手に、多分こうだろうと思われる方法で行っています（間違ってたら教えてくれるとありがたす）。
あと、パワースペクトルとそれなりに深い仲のインパルス応答の計算式も実装してみました。これに関しては、ソースに参考文献へのリンクを載せていますが、一読です。当初、インパルス応答に関しては、fibみたいな再帰関数でも書かなければ実装できそうになかったので（面倒だから）やめておこうと思ったのですが、当該の文献に行列の累乗でサクッと係数を計算できるアルゴリズムが紹介されていたので実装してみました。ライフハックｗとか、何でもかんでも、ハックｗハックｗというご時世ですが、こういうのを真のハックと言うのだよ、君。
動かしてみる
適当なＡＲ時系列を生成して、理論上のスペクトルと、シミュレーションデータから得られたスペクトルを比較してみました。
見ての通り、オレンジ色が理論上のやつで、青いギザギザが実測値です。これを見たら、誰だって間違ってると思う。俺もそー思う。そして、終わりのないデバッグが始まったのですが、どこをどう見ても、ぜってー間違ってる気がしねぇぜと思ったのでこれが正しいです。
さて、何が悪いのかというと、図が悪いのでした。Y軸を対数軸にすると：
このとおり、理論通りになっているように見えます。
いやね、「かがくしゃ」の人たちが、対数軸を使うのにはこういう理由があるんだね。単に見やすいからそうしていると思っていたのですが、ちゃんと図を書かないとあってるんだか間違ってるんだか分からないほどの差が出るようです。こういうシミュレーションをするときは、そもそもどうやって確認したらいいのかくらいは知っておいた方がいいようです。
教訓
よく知らない分野のプログラミングに首を突っ込むときは、ちゃんと下調べをすること。プロットしてみて、よさそうか、だめそうかで、バグっているか？動いているか？を判断するのは禁物。
なんか、どーでもいい感じになってしまいましたが、以上です。