ITの隊長のブログ

ITの隊長のブログです。Rubyを使って仕事しています。最近も色々やっているお（^ω^ = ^ω^）

時系列の勉強雑メモ

Python

ちゃんと理解したいので何度も読み直している。

  • 時系列分析の難しさ
    • 1990年1月1日の気温を予測する場合、無数にある1990年1月1日が母集団。母平均を推定しようとした場合、データは手元に1つ（1990年1月1日は現実に1つしかない）しか無いはずなのでむずい
  • なので、推測統計学の性質（？）を使って、DGP（データ生成過程）をモデリングできるならいけるんじゃないか？ってことを考える
    • 例えばサイコロで考えるとき、{1,2,...,6}は確率変数で{1/6,1/6,...,1/6}は母集団の確率分布
      • これを時系列に置き換える。1990年1月1日の気温（確率変数）は{1,2,...,6}で、確率分布は{1/6,1/6,...,1/6}とする
        • こうすると期待値も計算できるし、分散も計算できる
      • また時間が経過すると確率分布は変わるが同様に計算可能
    • 考え方は、「実際にある1つのデータ（1990年1月1日）を本来ならばあり得た確率変数の1つの実現値とみなし、もしも今日（1990年1月1日）という日が複数あれば次は違う値が得られたかもしれない。ということを想定する
  • データの表記方法（というか計算を理解したかった）
    • 期待値: \mu_t = E(y_t) 例えば1990年1月1日がもしも無数にあったら、平均的に \mu_t となる。ということを表す
    • 自己共分散: \gamma_{1t} = Conv(y_{t}, y_{t-1}) = E[(y_t - \mu_t )(y_{t-1} - \mu_{t-1})]
    • 自己相関: p_{kt} = Corr(y_t, y_{t-k}) = \frac{Conv(y_t, y_{t-k})}{\sqrt{Var(y_t)Var(y_{t-k})}}

自己相関がまでがでてきたので、前の記事に書いたコードと突き合わせたい。

www.aipacommander.com

今頃気づいたけど、参考にした記事と数式が違うな

f:id:aipacommander:20200414001359p:plain

引用元: https://momonoki2017.blogspot.com/2018/03/python7.html

引用元は全体から算出（？）、今見ている本はt時点とk時点の自己相関を算出しているという認識だから違う？？？

qiita.com

naruhodo. kはラグなのか（記事はhだけど）、nはデータ数。

あれ？じゃあt時点の期待値ってなんで全体のデータを使っているんだろ？？？

うーん、、、謎

statsmodelsのソースみてたらなんか引数で切り替えられるっぽい.

def autocorrelation2(y, lag=40, unbiased=False):
    # 自己相関係数を試しに計算してみた（Numpy利用）
    y_mean = np.mean(y) # 乗客数の平均値
    if unbiased:
        under = len(y) - np.arange(lag + 1)
    else:
        under = np.sum((y - y_mean) ** 2) 

    # ラグ1〜40の自己相関係数：-1.0〜1.0
    if unbiased:
        rk = [np.sum((y - y_mean)**2)]
        rk += [np.sum((y[k:] - y_mean) * (y[:-k] - y_mean)) for k in range(1, lag + 1)]
        rk /= under
    else:
        # ラグ0の自己相関係数：1.0
        rk = [1]
        rk += [np.sum((y[k:] - y_mean) * (y[:-k] - y_mean)) / under for k in range(1, lag + 1)]

    return np.array(rk)

全く同じ値にならなかったので、多分どっか間違えていると思うんだけど、1. 全データの数、2. lag=k時点の数で除算するという方法があった。

しかし参考にした記事の場合、1は「1. 共通の値（全データで算出した分散）で除算」になるので？になっている。

うーーーー、わからないを纏める（数式〜〜〜〜〜〜）

  • 書籍が指す分母とネット記事が指す分母が違うのはよくわからない（これはよく読んでいないから不明なだけかもしれないが、statsmodelsの実装も気になっている）

その次の偏自己相関も理解。

  • 端折るが、t時点とt-2時点前を比較したいのに、t-1時点のデータがあるので正確な自己相関を計算できない。から、t-1時点のデータを除去したい。という算出方法
  • \hat{y_t} = αy_{t-1} と書き、係数αは E[ (y_t\ - \hat{y_t})^2 ] を最小となるように選ぶ。こうすることで、t-1時点のデータから表現することができなかった残りが計算できる

実装したら理解できると思ったんだけど先が長い。。。

時系列の勉強雑メモ

Python

これまでゴリ推してきた時系列をなんとか数式含め脳内補完できないかなといことで再勉強中です。

早速ですが、ホワイトノイズは自己相関が0らしいので、本当か試してみる。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


dlen = 1024 #ノイズデータのデータ長
mean = 0.0  #ノイズの平均値
std  = 1.0  #ノイズの分散

y = np.random.normal(mean, std, dlen)

_, ax = plt.subplots(figsize=(30, 5))
ax.plot(y)
plt.show()

f:id:aipacommander:20200410021846p:plain

自己相関係数を計算する関数をパクって用意。

def autocorrelation(y, lag=40):
    # 自己相関係数を試しに計算してみた（Numpy利用）
    y_mean = np.mean(y) # 乗客数の平均値
    under = np.sum((y - y_mean) ** 2) 

    # ラグ0の自己相関係数：1.0
    rk = [1]

    # ラグ1〜40の自己相関係数：-1.0〜1.0
    rk += [np.sum((y[k:] - y_mean) * (y[:-k] - y_mean)) / under for k in range(1, lag + 1)]

    return np.array(rk)

momonoki2017.blogspot.com

a = autocorrelation(y, lag=40)
_, ax = plt.subplots(figsize=(100, 20))
ax.stem(np.arange(len(a)), a, use_line_collection=True)
plt.show()

f:id:aipacommander:20200410021956p:plain

さてはて、ここからが言語化です。雑メモ。

lagをもう少しわかりやすくしたさ

  • 今回のデータは1024日分のデータを（ランダムで生成）用意
  • コレログラムの1個目はlag=1のときのデータ　→ k=1のズレで算出したもの
    • 2日 ~ 1024日目と1日 ~ 1023日目で算出
    • 30日 ~ 1024日目と1日 ~ 994日目で算出
    • 35日 ~ 1024日目と1日 ~ 989日目で算出
  • 求めたいものって何だっけ？ → 周期とかズレでどの程度相関があるかそれぞれ確認する
  • いろいろ情報漁っている中で、6個目は「このコレログラムを素直に受け取れば現在と6日前までは自己相関を表している」という説明があったけど、ほんとう？？?
    • lag=6のときのデータは、6日 ~ 1024日目と1日 ~ 1018日目で算出
    • やっぱりギャップがある. 現在のデータ is 何？ と仮に正しいとしても、計算は点ではなく範囲で計算している（つまり6日目の値ではなく6日目〜のデータを使っている）
    • ん？ 6日前 か → 1日 ~ 1018日目 これは確かに6日前かな？（6日前までのデータ）
    • 現在とさしているのは、現在（1024日目）　〜 6日目まで
    • わかったかもしれん、1024日目は1024 - 6 = 1018日目なので、現在と6日前のデータで計算している（計算途中）
    • → これをズラした分だけ行う。ということか

これ、一年前もやった気がする。。。忘れたイメージを取り戻すの大変orz

seaborn.pairplotで左と下にでてくるlabelのrotationを変更したい

Python

とてもめんどかった（探索が）。ドキュメントを探してたらmatplotlibのようなクラスを扱えることを知りできた。

_x_columns =['width', 'height']
g = sns.pairplot(df[_x_columns], height=8, corner=True, hue='sex')

for ax in g.axes.flat:
    if ax is None:
        continue

    ax.set_xlabel(ax.get_xlabel(), rotation=90)  # 下側のラベル
    ax.set_ylabel(ax.get_ylabel(), rotation=0)  # 左側のラベル

fadeInとfadeOutを実装して死んだ

JavaScript

ウンコード

export function fade(node, duration, _type) {
    // style属性にdisplay: noneが設定されていたとき
    console.log(node.style.display);
    if (node.style.display === 'none') {
      node.style.display = 'none';
    } else {
      node.style.display = 'block';
    }
    node.style.opacity = 0;
    if (_type === 'out') {
      node.style.opacity = 1;
    }
  
    var start = performance.now();
    let tick = function(timestamp) {
      // イージング計算式（linear）
      var easing = (timestamp - start) / duration;
      if (_type === 'out') {
        easing = 1 - easing;
      }
  
      // opacityが1を超えないように
      console.log(node.style.opacity);
      if (_type === 'in') {
        node.style.opacity = Math.min(easing, 1);
      } else {
        node.style.opacity = Math.max(easing, 0);
      }
  
      if (_type === 'in') {
        // opacityが1より小さいとき
        if (easing < 1) {
          requestAnimationFrame(tick);
          node.style.display = 'block';
        } else {
          node.style.opacity = '';
        }
      } else {
        if (easing > 0) {
          requestAnimationFrame(tick);
        } else {
          node.style.opacity = '';
          node.style.display = 'none';
        }
      }
    }
    requestAnimationFrame(tick);
  }

jQueryなら4行ぐらいでおわりなのに....