Technology Topics by Brains

ブレインズテクノロジーの研究開発機関「未来工場」で働くエンジニアが、先端オープン技術、機械学習×データ分析(異常検知、予兆検知)に関する取組みをご紹介します。

WindowsのDockerで慣れないこと(実践編)

Neuron開発チームの木村です。

Neuronは、ほとんどの場合Windowsにインストールされるため、テスト環境をWindows上に構築しています。 このテスト環境構築を効率化するため、docker上でのWindows利用を模索しています。

前回に引き続き、 Windows上のdockerでイメージをビルドする際・コンテナを動かす際に慣れないことを紹介します。

前回検討した導入方法から、dockerは「Hyper-Vコンテナ + LCOW」の形態で利用しています。 (そのため、Docker for Windowsを利用する場合とは異なる状況であることに注意ください。)

Windowsのバージョンは、
Windows 10 Pro バージョン 1803 (Version 10.0.17134.165)、

dockerのバージョンは、

Client:
 Version:           master-dockerproject-2018-07-18
 API version:       1.38
 Go version:        go1.10.3
 Git commit:        48fbb12b
 Built:             Wed Jul 18 23:51:23 2018
 OS/Arch:           windows/amd64
 Experimental:      false

Server:
 Engine:
  Version:          master-dockerproject-2018-07-18
  API version:      1.38 (minimum version 1.24)
  Go version:       go1.10.3
  Git commit:       7f91801
  Built:            Thu Jul 19 00:00:43 2018
  OS/Arch:          windows/amd64
  Experimental:     true

LCOWのバージョンは、4.14.29-0aea33bcです。

以下、OSがWindowsであるコンテナを「WindowsOSコンテナ」、OSがLinuxであるコンテナを「LinuxOSコンテナ」と呼びます。

目次

イメージビルド編

multi-stage buildsで、名前に半角スペースを含むファイルをCOPYできない

multi-stage buildsでイメージビルドする際、名前に半角スペースを含むファイル(または、それを含むディレクトリ)を、WindowsのイメージにCOPYしようとすると、ビルドが以下のエラーで止まります。

COPY failed: file does not exist

よくある半角スペースの問題と考え、"で括るために
COPY --from=builder ["src/sample file", "."]のJSON形式で指定しても、今度は別のエラーで止まります。

COPY failed: Forbidden path outside the build context

昔からあるファイル・ディレクトリ名の半角スペース問題に、久しぶりに遭遇しました。(厳密には別の問題のようですが…)


RUN/CMDで使われるシェルを意識する

RUN/CMDにシェル形式でコマンドを指定する際、コマンド実行に使われるシェルがcmd(コマンドプロンプト)なのかpowershellなのかを意識することが重要です。powershell特有のコマンドを使う際は自然と意識すると思いますが、このふたつは環境変数の取得方法がだいぶ異なるため、javaのような自分でパスを通すコマンドを使う時に注意が必要です。

cmd powershell
環境変数VARの取得方法 %VAR% $env:VAR
存在しない環境変数NONEから取得される値 %NONE%という文字列 空文字列

シェル形式で書かれたRUN/CMDのコマンド実行に使われるシェルは、デフォルトがcmd(cmd /S /C)ですが、SHELLで変更することができます。
使用するイメージのDockerfileを読むか、docker inspectでイメージのCmdプロパティからシェルを特定するとよいでしょう。たとえば、openjdk:8u171-jdk-nanoserverpowershellをシェルに指定しています。


コンテナ実行編

[LinuxOSコンテナ] メモリ割り当てが変更できない

LinuxOSコンテナに対して、docker run時に--memoryオプションから割当メモリを指定しても無視され、 以下のように、約1GBのメモリが割り当てられます。

PS C:\> docker run --rm --memory 2g alpine:3.8 cat /proc/meminfo
MemTotal:         985568 kB
MemFree:          870932 kB
MemAvailable:     805808 kB
...

WindowsOSコンテナに対しては、以下のように、指定通り約2GBのメモリが割当たります。

PS C:\> docker run --rm --memory 2g microsoft/windowsservercore:1803 systeminfo
...
Total Physical Memory:     2,559 MB
Available Physical Memory: 2,194 MB
...

mobyのissueを読む感じでは、 LCOWを使用しているために起きる問題のようです…


[LinuxOSコンテナ] user指定が効かない

LinuxOSコンテナに対して、docker run時に--userオプションを指定しても無視され、 以下のように、rootユーザで実行されます。

PS C:\> docker run --rm --user guest alpine:3.8 whoami
root

LCOW使用時の既知のバグのようです。 イメージビルド時のUSER指定も効かないため、LinuxOSコンテナは常にrootで動くことになります。

そして、この影響なのか、コンテナ内のディレクトリ・ファイルの所有者とグループが大体root:rootになってしまいます。


[LinuxOSコンテナ] mountされたvolumeに対してchmod/chownが効かない

「イメージビルド時のVOLUME」や「docker run時の--volumeオプション」でvolumeに指定されたディレクトリ以下に対して、chmod/chownが効きません。
以下の例のように、chmodを行っても、終了コードが0であるにも関わらずパーミッションが変更されません(例ではjenkins/jenkins:2.134イメージを使用)。

root@fcdd7ae24179:/var/jenkins_home# ls -la
total 0
-r-xr-xr-x 1 root root  102 Jul 30 05:02 copy_reference_file.log
drwxrwxrwx 1 root root 4096 Jul 30 05:02 init.groovy.d
root@fcdd7ae24179:/var/jenkins_home# chmod 755 copy_reference_file.log
root@fcdd7ae24179:/var/jenkins_home# echo $?
0
root@fcdd7ae24179:/var/jenkins_home# ls -la
total 0
-r-xr-xr-x 1 root root  102 Jul 30 05:02 copy_reference_file.log
drwxrwxrwx 1 root root 4096 Jul 30 05:02 init.groovy.d

既知の問題ではあるようですが、回避策は特にない状況です。
このことが起因する問題には色々当たりましたが、大きかったのはjenkinsでリポジトリをダウンロードする際にエラーが起きることでした。 volume指定を使わないことでエラーを回避しましたが、docker cpでjenkins設定の永続化を行うことになってしまいました…


稼働中のコンテナに対してdocker cpできない

稼働中のコンテナに対してdocker cpを実行すると、以下のエラーが発生してコピーができません。

PS C:\Users\user\Documents> docker cp 7961ee45d606:/var .
Error response from daemon: filesystem operations against a running Hyper-V container are not supported

なお、停止中のコンテナに対しては正常にコピーできます。


おわりに

Windows上のdockerにて、イメージビルドやコンテナ実行の際に慣れないことを紹介しました。 LinuxOSコンテナについて、「volume周りの問題」と「メモリが最大1GB」がかなり厄介な問題です。

上記の他にも以下のようなことがありました。

  • コンテナ自体を作成できないイメージがある
    • 例えば、openjdk:8u171-jdk-nanoserver
    • しかし、openjdk:8u171-jdk-windowsservercore はコンテナ実行可能
  • LinuxOSコンテナ用のDockerfileを、WindowsOS用に書き換えるには細かいところも気にしないとならない
    • javaのclasspath区切りは、Windowsでは:ではなく;、など
  • dockerdを再起動すると、元々稼働していたコンテナがstartできなくなる

やはり、experimentalなLCOWを使ってLinuxOSコンテナを扱うのは、まだ早いのでしょうか。 WindowsOSコンテナのみに絞ればよいかもしれませんが、そうするとLCOWを使う理由もないですし…

docker上でのWindows利用について、LCOWは引き続き試してゆきますが、 現状では、Docker for Windowsが可能な範囲で模索してゆくほうがよさそうです。


ブレインズテクノロジーではエンジニアを募集中です
採用ページはこちら

AngularでPlotlyのグラフを描画してみた

ブレインズテクノロジーの加藤です。

例年以上に暑い日々が続いていますが、如何お過ごしでしょうか。
ブレインズテクノロジーでは、そんな暑い夏にも負けないくらいの空前のボルダリング熱波が巻き起こっています。
週に1回、有志でジムに通っていますが、いい運動になりますし、課題をクリアしたときの達成感はひとしおです。

さて今回は、AngularでPlotlyを動かしてグラフ描画する方法についてまとめます。

はじめに

Angularを触ったことのない人には、下記のチュートリアルがおすすめです。
ヒーローがどうのこうのという謎アプリを作らされますが、何となくのイメージが掴めると思います。

Angular 日本語ドキュメンテーション

新しいアプリケーションを作成

下記のコマンドで、新しいアプリケーションを作ります。
名前は何でもいいですが、今回はとりあえず「angular-plotly」とします。

$ ng new angular-plotly

Plotlyファイルのimport

先ほど作った「angular-plotly」下に移動し、npm経由でplotlyの型定義ファイルを取得します。

$ cd angular-plotly
$ npm install plotly.js --save
$ npm install @types/plotly.js --save

Plotlyによるグラフ描画の実装

シンプルに下記サンプルの棒グラフ「Basic Bar Chart」を描画します。
グラフ出力時のデータなどは、サイトの情報をそのまま用いることにします。
Javascript Graphing Library D3.js-based Bar Charts | Examples | Plotly

ここからは、実際のtypescriptのコードとhtmlファイルを記載します。

  • src/app/app.component.ts
import { Component, OnInit } from '@angular/core';
import * as Plotly from 'plotly.js';

@Component({
  selector: 'app-root',
  templateUrl: './app.component.html',
  styleUrls: ['./app.component.css']
})

export class AppComponent implements OnInit {
  ngOnInit() {
    const data = [
      {
        x: ['giraffes', 'orangutans', 'monkeys'],
        y: [20, 14, 23],
        type: 'bar'
      } as Plotly.ScatterData
    ];
    Plotly.newPlot('myDiv', data);
  }
}
  • src/app/app.component.html
<div id='myDiv'></div>

解説

簡単にコードを解説します。

import * as Plotly from 'plotly.js';

で、先程npm経由で取得したplotlyを呼べるようにしています。

また、

import { Component, OnInit } from '@angular/core';

で、起動時の動作を定義するインターフェース「OnInit」をimportしています。
「OnInit」をimplementsした「ngOnInit」関数にてPlotly描画の定義をすることにより、Angular起動と同時にグラフ描画処理が実行されます。

export class AppComponent implements OnInit {
  ngOnInit() {
    const data = [
      {
        x: ['giraffes', 'orangutans', 'monkeys'],
        y: [20, 14, 23],
        type: 'bar'
      } as Plotly.ScatterData
    ];
    Plotly.newPlot('myDiv', data);
  }
}

「ngOnInit」関数で定義している内容は、前述の「Basic Bar Chart」のJSコードを一部改変したものです。
Plotlyで描画する際に、dataの型を「Plotly.ScatterData」に変換しています。

Angular起動

下記のコマンドを実行することにより、Angularが起動します。

$ ng serve --open

・・・が、画面に何も描写されないと思います。

デベロッパーツールで確認すると、下記のようなエラーが出ています。

Uncaught ReferenceError: global is not defined

回避策として、src/polyfills.tsの文末に以下を記載します。

 // "global is not defined"の対応
(window as any).global = window;

この状態で再起動すると、無事にグラフが描画されました。

f:id:kato_takayuki:20180723171618p:plain

おわりに

今回は最もシンプルなグラフを生成しましたが、同じ要領でいろいろなグラフが生成出来ると思います。

おわり。


ブレインズテクノロジーでは「共に成長できる仲間」を募集中です。
採用ページはこちら

Pythonの機械学習ライブラリtslearnを使った時系列データのクラスタリング

こんにちは、ブレインズテクノロジーの柏木です。

今回はPythonで扱える機械学習ライブラリのtslearnを使って、時系列データをクラスタリングしていきたいと思います。

github.com

tslearnとは

時系列分析のための機械学習ツールを提供するPythonパッケージで、scikit-learnをベースとして作られているみたいです。 主な機能として、クラスタリング、教師ありの分類、複数の時系列を重ねた際の重心の計算ができたりします。 今回使用するに至った一番のモチベーションは、波形や振動などの時系列データに対してクラスタリングできるというところです。

tslearnインストール

pipコマンドでインストールできます。

pip install tslearn

Kshapeというクラスタリング手法

今回tslearnで使用するモジュールとして、Kshapeというクラスタリング手法を時系列データに適用していきたいと思います。 Kshapeは2015年に下記の論文で提唱された方法で、以下の流れで実行されるアルゴリズムになります。

  1. 相互相関測定に基づいた距離尺度を使う(Shape-based distance: SBD)
  2. 1.に基づいた時系列クラスタの重心を計算(時系列の形状を保持する新しい重心ベースのクラスタリングアルゴリズム)
  3. 各クラスタに分割する方法は一般的なkmeansと同じだが、距離尺度や重心を計算する時に上記1と2を使う

細かい数式や理論については論文を確認下さい。 個人的に、内容がしっかり記載されていてわかりやすい論文だと思います。

Paparrizos, John and Luis Gravano. “k-Shape: Efficient and Accurate Clustering of Time Series.” SIGMOD Record 45 (2015): 68.

今回使用するデータセットについて

今回使うデータセットはthe UEA & UCR Time Series Classification RepositoryにあるDataset: SonyAIBORobotSurface1を加工したものになります。 このサイトには様々な時系列データがありますので、色々見てみるのも楽しいかもしれないです。

使用したサンプルデータはこちらに置いておきます。 あるワーク単位でファイルが10個あり、1秒単位のセンサーデータを想定したものになります。

時系列データのクラスタリング

それでは、時系列データのクラスタリングをしていきたいと思います。 実際のPythonコードを記載します。 まず、ライブラリのインポートです。

#Python 3.6.4
import pandas as pd
import numpy as np
import glob
from tslearn.clustering import KShape
from tslearn.preprocessing import TimeSeriesScalerMeanVariance

import matplotlib.pyplot as plt
%matplotlib inline

次に、いくつか関数を定義しておきます。

def read_file_to_dataframe(filenames):
    #ファイルデータを読み込み、データフレームを返す
    dfs = []
    for filename in filenames:
        original_df = pd.read_csv(filename, index_col=None, header=0)
        dfs.append(original_df)
    return dfs

def time_series_data_to_array(dataframes, target_col=''):
    #データフレームを読み込み、それらを時系列の配列にする
    tsdata = []
    for i, df in enumerate(dataframes):
        tsdata.append(df[target_col].values.tolist()[:])
        #それぞれの時系列データの最大の長さを確認
        len_max = 0
        for ts in tsdata:
            if len(ts) > len_max:
                len_max = len(ts)
        #時系列データの長さを揃えるために、最後のデータを付け加える
        for i, ts in enumerate(tsdata):
            len_add = len_max - len(ts)
            tsdata[i] = ts + [ts[-1]] * len_add
    
    tsdata = np.array(tsdata)
    return tsdata

def transform_vector(time_series_array):
    #ベクトルに変換
    stack_list = []
    for j in range(len(time_series_array)):
        data = np.array(time_series_array[j])
        data = data.reshape((1, len(data))).T
        stack_list.append(data)
    #一次元配列にする
    stack_data = np.stack(stack_list, axis=0)
    return stack_data


filenames = sorted(glob.glob('sample_data/sample_data*.csv'))
df = read_file_to_dataframe(filenames=filenames)
tsdata = time_series_data_to_array(dataframes=df, target_col='data')
stack_data = transform_vector(time_series_array=tsdata)

ここでは、使用するファイルを読み込んで配列に変換します。 この時に、仮に時系列のデータの長さがファイル毎に異なる場合、長さを揃える必要があります。 そのため、今回は時系列上一番最後の値で埋め合わせを行なっています。 長さを揃えた後は、一次元のベクトルに変換しています。

それでは、加工したデータを学習して実際にクラスタリングしていきます。

seed = 0
np.random.seed(seed)
#相互相関を計算するために、正規化する必要があります。
#TimeSeriesScalerMeanVarianceがデータを正規化してくれるクラスになります。
stack_data = TimeSeriesScalerMeanVariance(mu=0.0, std=1.0).fit_transform(stack_data)

#KShapeクラスのインスタンス化
ks = KShape(n_clusters=2, n_init=10, verbose=True, random_state=seed)
y_pred = ks.fit_predict(stack_data)

#クラスタリングして可視化
plt.figure(figsize=(16,9))
for yi in range(2):
    plt.subplot(2, 1, 1 + yi)
    for xx in stack_data[y_pred == yi]:
        plt.plot(xx.ravel(), "k-", alpha=.2)
    #plt.plot(ks.cluster_centers_[yi].ravel(), "r-")
    plt.title("Cluster %d" % (yi + 1))

plt.tight_layout()
plt.show()

以下、結果をプロットしたものになります。

f:id:brains_mkashiwagi:20180607143544p:plain

KShapeクラスの引数n_clustersでクラスタ数を指定することができ、今回はn_clusters=2としました。 実際のデータセットのページを見てみると、AIBOがカーペットを歩いている場合とコンクリートを歩いている場合の二種類のデータとなっています。 それを踏まえて、結果を見てみると良い感じに分離できているように見えます。 今回はクラスタ数を2としましたが、クラスタリングする際によくあるクラスタ数をどのように決めたら良いかという問題について、エルボー法で確認した結果を次に示します。

エルボー法によるクラスターの数の計算

エルボー法とは... 各点から割り当てられたクラスタ中心との距離の二乗の合計をクラスタ内誤差平方和(SSE)として計算します。 クラスタ数を変えて、それぞれのSSE値をプロットし、「肘」のように曲がった点を最適なクラスタ数とする方法です。 また、SSEが小さいほど歪みのない良いモデルであり、うまくクラスタリングできていることになります。

エルボー法について、wikipediaに記載されています。
Elbow method (clustering) - Wikipedia

distortions = []

#1~10クラスタまで計算 
for i in range(1,11):
    ks = KShape(n_clusters=i, n_init=10, verbose=True, random_state=seed)
    #クラスタリングの計算を実行
    ks.fit(stack_data)
    #ks.fitするとks.inertia_が得られる
    #inertia_でSSEを取得できる
    distortions.append(ks.inertia_)

plt.plot(range(1,11), distortions, marker='o')
plt.xlabel('Number of clusters')
plt.ylabel('Distortion')
plt.show()

f:id:brains_mkashiwagi:20180607152931p:plain

n_clusters=2で肘のように曲がっているので、今回はクラスタ数2で正解のようです。 ただし、一般的にエルボー法ではこのようにはっきりと曲がりがある場合は少なく、あまり精度はよくないようです。 他にもクラスタ数を決める方法は色々ありシルエット分析、X-means、GAP法などがあります。 興味がある方はX-means、GAP法は論文がありますので、参考文献をご覧下さい。

クラスタリングの決め方について、wikipediaに記載されています。
Determining the number of clusters in a data set - Wikipedia

まとめ

今回はtslearnを使って、時系列データのクラスタリングをしてみました。 振動や波形データをうまくクラスタリングできるので、非常に重宝すると思います。 tslearnにはKshape以外にもTimeSeriesKMeansという手法もあり、距離尺度としてDTW(動的時間伸縮法)やSoft-DTWを使ってクラスタリングをすることができます。さらに時系列に対して重心線を引いてくれたりもします。
今回はデータ点数も少なく振動も緩やかであるため、クラスタ数を決めることが容易でしたが、一般的に時系列データのクラスタ数を決めることは難しい問題だったりするので、その辺りを決める良い方法を模索していきたいと思います。


ブレインズテクノロジーでは「共に成長できる仲間」を募集中です。
採用ページはこちら

参考文献および参考ウェブサイト

PythonでのARIMAモデルを使った時系列データの予測の基礎[後編]

こんにちは。ブレインズテクノロジーの岩城です。

こちらの記事は、2部構成でお送りしている、Pythonを使用してのARIMAモデルの作成・予測の流れの整理の後半です。

前半の記事では、ARIMAモデルの基礎を簡単におさらいしました。
後半の記事では、Pythonを使って実際にコードを確認しながら、ARIMAモデルを作成し、予測するまでの流れを説明します。

ARIMAモデルの復習が必要な方は、是非こちらの記事も参考にしてください。

blog.brains-tech.co.jp

今回も、少し理論の話が出てきます。 理論の部分は前回に引き続き、こちらの教科書を参考にしています。

沖本竜義(2010), "経済・ファイナンスデータの計量時系列分析", 朝倉書店http://www.asakura.co.jp/G_12.php?isbn=ISBN978-4-254-12792-8

www.asakura.co.jp

きちんと理論を整理したい場合は、是非読んでみてください。

記事執筆時の実行環境は下記の通りです。

実行環境

  • MacBook Pro(13-inch, 2017, Four Thunderbolt 3 Ports)
  • プロセッサ: 3.1 GHz Intel Core i5
  • メモリ: 16 GB 2133 MHz LPDDR3
  • OS: macOS High Sierra (10.13.3)
  • Python 3.6.4
  • statsmodels 0.8.0
  • numpy 1.14.0
  • pandas 0.22.0

PythonでのARIMAモデル作成

説明の流れの整理

それでは、早速モデル作成に入っていきましょう。 先に全体の概要を整理しておきます。

  1. 使用するデータについて
  2. データの特徴の確認
  3. 単位根検定の実施
  4. 総当たり法でのモデル候補探索
  5. 良さそうなモデルの確認と評価
  6. 選択したモデルを用いた予測

はじめにどういうモデルになりそうか確認します。その後、いくつかのパラメータを組み合わせた総当たりでのモデル作成を行い、良さそうな候補を選択した上で、結果を確認します。最後に、今回の選択手法でもっとも良いと思われたモデルを使用して、予測を行います。

なお、一連の流れはJupyter notebook上で行いました。

1. 使用するデータについて

今回は元になるデータのモデルが分かっている上で、総当たりの手法で選択される結果について整理してみることにしました。
データはARIMA(0, 1, 1)として設計しています。差分を取る前の、元のデータの導出式は下記の通りとしました。

   y_t = 0.1 + y_{t-1} + \varepsilon_t - \varepsilon_{t-1}

この式は変形すると、

   y_t - y_{t-1} = \Delta y_t = 0.1 + \varepsilon_t - \varepsilon_{t-1}

となり、1階差分の系列が c = 0.1 \theta_1 = -1のMA(1)過程になるようになっています。誤差項は \mu = 0 \sigma^2 = 25  (\sigma = 5)の正規分布としました。

また、最後に日付を利用した予測をさせるので、ここでは架空の日付を作成してデータに付与しておきました。

ここまでのコードは下記の通りです。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
import datetime
%matplotlib inline

def create_data(data_length, delta=0, ar1 = 1, ma1 = -1, set_std=1, y0 = 0, random_seed=0):
    np.random.seed(random_seed)
    cur_y = y0
    e_m1 = np.random.normal(loc=0, scale=set_std)
    val_list = []

    for i in range(data_length):
        val_list.append(cur_y)
        e_0 = np.random.normal(loc=0, scale=set_std)
        cur_y = delta + ar1 * cur_y + e_0 + ma1 * e_m1
        e_m1 = e_0
        
    return val_list

# 架空の時間データを作成する。与えた年、月、日から1時間ごとの情報を作成する。
def create_dummy_date(start_info, data_num):
    """
    Args:
        start_info: List of integers, [year, month, day]
        data_num: int, number of data.
    Returns:
        List of datetime.
    """
        
    date_list = []
    time_info = datetime.datetime(year=start_info[0], month=start_info[1], day=start_info[2], hour=0)
    for i in range(data_num):
        date_list.append(time_info)
        time_info = time_info + datetime.timedelta(hours=1)
    return date_list

# データの作成
data_num = 2000
delta = 0.1
coef_a = 1
ma1 = -1
set_std = 5
random_seed=0

base_data = create_data(data_num, delta=delta,coef_a=coef_a, ma1=ma1, set_std=set_std,random_seed=random_seed)
linear_data = [i * delta for i in range(data_num)]

index_info = create_dummy_date([2000, 5,1], data_num)

df_all = pd.DataFrame(index=index_info, columns=["raw", "diff1", "linear"])
df_all["raw"] = base_data
df_all["linear"] = linear_data
df_all["diff1"] = df_all.raw - df_all.raw.shift()

fig, ax = plt.subplots( figsize=(6, 4))
df_all.plot(y=["raw", "linear"], ax=ax)

作成したグラフはこちらです。オレンジのラインが比較のための  y = 0.1t のグラフです。作成したデータはこのグラフに近い位置に存在しています。

f:id:brains_iwaki:20180604175326p:plain

記事の最後で、作成したモデルを用いた予測を行うので、ここでモデル作成用のデータを分けておきます。

df_base = df_all[:1000][["raw", "diff1"]]
raw_data = df_base.raw
diff1_data = df_base.diff1.dropna()

モデル作成用データの先頭3行と末尾3行を確認しましょう。

# データの始めの3行と末尾3行の確認
df_base.head(3).append(df_base.tail(3))
raw diff1
2000-05-01 00:00:00 0.000000 NaN
2000-05-01 01:00:00 -6.719476 -6.719476
2000-05-01 02:00:00 -3.726572 2.992904
2000-06-11 13:00:00 91.350700 1.562956
2000-06-11 14:00:00 85.241684 -6.109016
2000-06-11 15:00:00 89.289168 4.047484

このようなデータが、今回のモデル作成の対象になります。

2. データの特徴の確認

次に、作成したデータの時系列的な特徴を確認します。
以下に、元のデータと1階差分をとったデータそれぞれについて、生データ、自己相関関数(ACF)、偏自己相関関数(PCAF) のグラフを示します。

plt.rcParams["font.size"] = 10
fig, ax = plt.subplots(2, 3, figsize=(20, 8))
axes = ax.flatten()
diff_color = "#800000"

axes[0].plot(raw_data)
axes[0].set_title("raw_data")
fig = sm.graphics.tsa.plot_acf(raw_data, lags=40, ax=axes[1])
fig = sm.graphics.tsa.plot_pacf(raw_data, lags=40, ax=axes[2])

axes[3].plot(diff1_data.dropna(), color=diff_color, alpha=1)
axes[3].set_title("diff1")
fig = sm.graphics.tsa.plot_acf(diff1_data, lags=40, ax=axes[4], color=diff_color)
fig = sm.graphics.tsa.plot_pacf(diff1_data, lags=40, ax=axes[5], color=diff_color)
plt.savefig("../output/df_base_info.png", bbox_inches='tight', dpi=120)

f:id:brains_iwaki:20180604175615p:plain

自己相関関数(ACF, Autocorrelation function)と偏自己相関関数(PCAF, Parcial Autocorrelation function)は、任意の時刻のデータ y_tとそのk個前のデータ y_{t-k}の関係性を示す-1 から1の間の値です。

ACFは y_k y_{t-k}との間の線形的な相関の強さを、PCAFは y_k y_{t-k}に対するy_{t-k}以外の項( y_{t-k+1}など)の影響を除いた相関の強さを示しています。

細かな説明は参考文献をご参照いただくことにして、グラフから分かることを整理します。
その前に、元データの式と1階差分データの式は下記の通りでした。

   y_t = 0.1 + y_{t-1} + \varepsilon_t - \varepsilon_{t-1}
   y_t - y_{t-1} = \Delta y_t = 0.1 + \varepsilon_t - \varepsilon_{t-1}

まずは自己相関係数(ACF)の確認をしましょう。

ARのACFの特徴は、一つの項の影響が長期的に減衰しながら残ることです。 例えば、

  y_t= ay_{t-1} = a(ay_{t-2}) =  a^2 (ay_{t-3}) = ...

のように考えてみると分かりやすいかと思います。このように、ある時刻のデータは、ずっと離れた時刻のデータから表すことができます。

MAのACFの特徴は、MAの次数以降の項が0に近くなることです。これは例えばMA(1)なら、

   y_t =  \varepsilon_t - \theta_{1} \varepsilon_{t-1}
   y_{t-1} =  \varepsilon_{t-1} -  \theta_{1} \varepsilon_{t-2}
   y_{t-2} =  \varepsilon_{t-2} -  \theta_{1} \varepsilon_{t-3}

のように式を並べてみると、  y_t y_{t-2}は直接同じ項を共有していないことが確認できます。このため、直接的な相関関係はない、と考えることができます。
今回は、元データにARの項があり、差分データはMAのみになっているので、ともに期待通りのグラフになっています。

次に、偏自己相関係数(PCAF)です。

AR項の偏自己相関係数は、もともとARの式が次数pの過去の項から構成されていたので、次数pより後ろの項とのPCAFは0になります。 一方で、MA過程は反転可能であればAR(∞)と書き直すことができることもあり、偏自己相関係数も無限に残ることになります。

今回のグラフでは、両方ともMAの項を含んでいることもあり、急激に小さくなるのではなく、減衰していく様子が確認できます。

3. 単位根検定の実施

続いて、単位根検定を行ってみましょう。

その前に、ここで初めて単位根という言葉が出てきたので、まずは定義を確認します。

定義 5.1 (単位根過程)
原系列 y_tが非定常過程であり,差分系列 \Delta y_t = y_t - y_{t-1}が定常過程であるとき,過程は単位根過程(unit root process)といわれる.
(沖本竜義(2010), "経済・ファイナンスデータの計量時系列分析", P. 105, 朝倉書店より引用)

要は、元のデータが非定常過程であり、差分系列が定常過程であるような系列が単位根過程です。今回は1階差分によって定常の式が得られるようにデータを作成したので、検定の結果、単位根過程である可能性を示せることを期待します。

単位根検定を実施する背景には、データが単位根過程であるか、別の過程かによって、作成するモデルの特性や得られる予測が大きく異なってくることがあります。こちらも非常に重要なので、教科書をご参照ください。

単位根検定にもADF検定、PP検定、KPSS検定など、いくつか種類があります。 こちらも詳細は参考文献等をご参照ください。

今回はADF検定を用いることにしました。

ADF検定では、過程が単位根AR(p)であることを帰無仮説とした検定を行います。そのため、棄却されなければ単位根がないとはいえない、という結論になります。なお、単位根検定の際は、仮説のモデルが定数項とトレンドを含むか考慮する必要があります。今回はデータを与えるモデルをこちらで決めていたので、モデル通りの、定数を含む場合での結果を示します。

statsmodelsのadfullerはこちらを参考にしてください。
statsmodels.tsa.stattools.adfuller — statsmodels 0.9.0 documentation

それでは、元データへのADF検定を行いましょう。

# 単位根検定
adf_result = sm.tsa.stattools.adfuller(raw_data)
adf_result

出力結果は、

    (0.10136608604360935,
     0.9661542394395962,
     18,
     981,
     {'1%': -3.4370334797663844,
      '10%': -2.568341114581742,
      '5%': -2.8644907213150725},
     5938.189071549921)

となりました。一番初めの値が検定値、2番目がP値です。
期待通り、P値はおよそ0.966で、0.05よりも大きいので、単位根AR(p)過程でないとはいえないという検定結果が確認できました。

4. 総当たり法でのモデル候補探索

今回は総当たり法でモデルの候補を探します。

ちなみに、総当たりを用いずに、データからモデルを考える場合には、先ほどのACF、PCA、単位根検定などを利用して候補を考える方法があります。

その場合は、次数が減衰する位置の確認や検定を組み合わせながら、モデルを考えていきます。
今回の例では1階差分のACFが2次以降で0近くになったので、ARIMA(p,1,1) (p≥1)が候補になりそうですね。

総当たり法では単純に、様々なパラメータでモデルを作成し、AIC等を用いて有力だと思われるモデルを探します。

今回は、p={0, 1, 2}、q={0, 1, 2}, d={0, 1}の各パラメータの組み合わせを試すことにします。

結果の確認のため、作成したモデルのAICを一つのデータフレームに集約させます。
コードは下記の通りです。

# シンプルな総当たり
p_num = 3
d_num = 2
q_num = 3
trial_num = p_num * d_num * q_num
use_columns = ["p", "d", "q", "AIC"]
step_cnt = 0
train_data = raw_data

info_df = pd.DataFrame(index=range(trial_num-2), columns=use_columns)

for p in range(p_num):
    for d in range(d_num):
        for q in range(q_num):
            if p == 0 and q == 0:
                continue
            model = sm.tsa.statespace.SARIMAX(train_data, order=(p, d, q), enforce_invertibility=False, trend='t')
            param_list = [p, d, q]
            for name, val in zip(use_columns[:-1], param_list):
                info_df.iloc[step_cnt][name] = val
            try:
                result = model.fit(disp=False)
                # print(result.aic)
                info_df.iloc[step_cnt]["AIC"] = result.aic
            except:
                pass
            step_cnt += 1
            print("Finish trial No. {}, Param: {}".format(step_cnt, param_list))

AICが低い順番に並べ変えた結果が下記の表です。

info_df.sort_values(by="AIC")
p d q AIC
12 2 0 2 6070.66
3 0 1 2 6085.99
14 2 1 1 6086.7
8 1 1 1 6092.45
5 1 0 1 6126.84
6 1 0 2 6134.79
2 0 1 1 6149.16
11 2 0 1 6218.62
13 2 1 0 6295.7
10 2 0 0 6305.9
4 1 0 0 6423.68
7 1 1 0 6433.19
1 0 0 2 6439.94
0 0 0 1 6511.52
9 1 1 2 NaN
15 2 1 2 NaN

AIC最小のモデルは(2, 0, 2)、2番目が(0, 1, 2)、3番目が(2, 1, 1)、元のモデルのパラメータである(0, 1, 1)は7番目でした。

それではこれらの4つの組み合わせに対して改めてモデルを作成し、結果を確認してみましょう。

5. 良さそうなモデルの確認と評価

先ほどのパラメータを使用したモデルを再作成し、結果を確認します。

良いモデルかどうかは、作成したモデルから出力される結果と実際の値との残差が、平均0、自己相関係数が全て0のホワイトノイズになっているかどうかで判断します。

全ての自己相関係数が0かどうかを判断するには、Ljung-Box検定を使用します。
この検定では、自己相関係数に対して \rho_1 = \rho_2 = ... = \rho_m= 0 (1 \leq m \leq k, mは整数 )という帰無仮説を立てます。得られたP値から、

  • P値が0.05より小さく、帰無仮説が棄却されれば、kまでの間の少なくとも一つは自己相関係数が0ではない。
  • P値が0.05より大きく、帰無仮説が棄却されなかった場合は、kまでの間の全ての自己相関係数が0ではないとはいえない。

という検定結果を得ることができます。

モデル作成結果の中にもLjung-Box検定の結果は記載されていますが、今回は明示的に検定を行うことにしました。検定の結果を見やすくするために、下記の関数を使用することにします。

# 残差に対するLjung-Box検定
def make_lbtest_df(test_data, lags=10):
    test_result = sm.stats.diagnostic.acorr_ljungbox(test_data, lags=lags)
    lbtest_df = pd.DataFrame({
        "Q": test_result[0],"p-value": test_result[1]})
    lbtest_df = lbtest_df[["Q", "p-value"]]
    return lbtest_df

(1) (2, 0, 2)でのモデル作成

それでは、まずは総当たりでAIC最小だった(2, 0, 2)のモデルを再度作成し、結果を確認します。

model202 =  sm.tsa.statespace.SARIMAX(train_data, order=(2, 0, 2), trend='c')
result202 = model202.fit(disp=False)
print(result202.summary())

出力結果は下記のようになりました。

                               Statespace Model Results                           
    ==============================================================================
    Dep. Variable:                    raw   No. Observations:                 1000
    Model:               SARIMAX(2, 0, 2)   Log Likelihood               -3055.083
    Date:                Mon, 04 Jun 2018   AIC                           6122.166
    Time:                        17:07:09   BIC                           6151.612
    Sample:                    05-01-2000   HQIC                          6133.357
                             - 06-11-2000                                         
    Covariance Type:                  opg                                         
    ==============================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
    ------------------------------------------------------------------------------
    intercept      0.0084      0.042      0.201      0.841      -0.074       0.090
    ar.L1          0.0136      0.017      0.802      0.423      -0.020       0.047
    ar.L2          0.9862      0.017     58.070      0.000       0.953       1.019
    ma.L1          0.1036      0.020      5.301      0.000       0.065       0.142
    ma.L2         -0.8833      0.018    -49.170      0.000      -0.919      -0.848
    sigma2        26.2190      1.187     22.090      0.000      23.893      28.545
    ===================================================================================
    Ljung-Box (Q):                       43.34   Jarque-Bera (JB):                 0.26
    Prob(Q):                              0.33   Prob(JB):                         0.88
    Heteroskedasticity (H):               1.01   Skew:                             0.03
    Prob(H) (two-sided):                  0.90   Kurtosis:                         2.95
    ===================================================================================
    
    Warnings:
    [1] Covariance matrix calculated using the outer product of gradients (complex-step).

平均値と残差を用いたLjung-Box検定の結果を確認しましょう。

resid202 = result202.resid
print(resid202.mean())
lbtest202 = make_lbtest_df(resid202)
lbtest202.index = range(1, 11)
lbtest202
Mean of residual error is 0.7901548893874051
Q p-value
1 5.794679 0.016075
2 6.866166 0.032287
3 7.270120 0.063769
4 13.793465 0.007984
5 14.626877 0.012081
6 16.375478 0.011874
7 20.259926 0.005035
8 20.272314 0.009353
9 20.306049 0.016115
10 20.725448 0.023091

残差の平均は0.8程度でした。また、k=10までのLjung-Box検定結果のP値はいずれも0.05より小さいので、自己相関係数のうち少なくとも一つは0ではないという検定結果となります。

以上の結果から、このモデルは今のところ必ずしも良いモデルとはいえなさそうです。
(p, d, q) = (2, 0, 2)というパラメータの通り、差分を一度も取らないモデルである、ということも関係しているのでしょうか。

(2) (0, 1, 2)でのモデル作成

次は、AICが2番目に小さかった、(p, d, q) = (0, 1, 2)の結果を確認します。

model012 =  sm.tsa.statespace.SARIMAX(train_data, order=(0, 1, 2), trend='c')
result012 = model012.fit(disp=False)
print(result012.summary())

resid012 = result012.resid
print("\nMean of residual error is {}".format(resid012.mean()))
lbtest012 = make_lbtest_df(resid012)
lbtest012.index = range(1, 11)
lbtest012
    /Users/brainsiwaki/.pyenv/versions/anaconda3-5.1.0/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
      "Check mle_retvals", ConvergenceWarning)


                               Statespace Model Results                           
    ==============================================================================
    Dep. Variable:                    raw   No. Observations:                 1000
    Model:               SARIMAX(0, 1, 2)   Log Likelihood               -3015.781
    Date:                Mon, 04 Jun 2018   AIC                           6039.561
    Time:                        18:10:22   BIC                           6059.192
    Sample:                    05-01-2000   HQIC                          6047.022
                             - 06-11-2000                                         
    Covariance Type:                  opg                                         
    ==============================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
    ------------------------------------------------------------------------------
    intercept      0.0996      0.002     64.571      0.000       0.097       0.103
    ma.L1         -0.9864      0.031    -31.372      0.000      -1.048      -0.925
    ma.L2         -0.0052      0.032     -0.162      0.871      -0.068       0.057
    sigma2        24.3259      1.105     22.019      0.000      22.161      26.491
    ===================================================================================
    Ljung-Box (Q):                       34.17   Jarque-Bera (JB):                 0.22
    Prob(Q):                              0.73   Prob(JB):                         0.90
    Heteroskedasticity (H):               0.98   Skew:                             0.03
    Prob(H) (two-sided):                  0.86   Kurtosis:                         2.96
    ===================================================================================
    
    Warnings:
    [1] Covariance matrix calculated using the outer product of gradients (complex-step).
    
    Mean of residual error is -0.08916347484994207
Q p-value
1 1.313842 0.251700
2 1.396998 0.497331
3 2.765295 0.429245
4 4.687174 0.320927
5 4.750206 0.447120
6 4.932049 0.552557
7 13.140888 0.068749
8 13.630489 0.091920
9 14.180136 0.116060
10 14.187540 0.164608

残差の平均はおよそ-0.09、Ljung-Box検定ではいずれもP値が0.05よりも大きく、帰無仮説は棄却されていません。そのため、残差の自己相関係数がいずれも0ではないとはいえない、という検定結果になります。

モデルとしては悪くなさそうですが、モデル作成時にConvergenceWarningが発生しています。mle_retvalsを確認するようにメッセージが出ているので確認しておきます。

print(result012.mle_retvals)

出力結果はこちらです。

{'fopt': 3.0157805026031124, 'gopt': array([-0.14644907, -0.0001639 , -0.00720694, -0.00082886]), 'fcalls': 355, 'warnflag': 1, 'converged': False}

最後の'converged'がFalseになっています。最尤推定の最適化の際に収束性しないような問題がありそうです。

(3) (2, 1, 1)でのモデル作成

次は、AICが3番目に小さかった、(p, d, q) = (2, 1, 1)の結果を確認します。

model211 =  sm.tsa.statespace.SARIMAX(train_data, order=(2, 1, 1), trend='c')
result211 = model211.fit(disp=False)
print(result211.summary())

resid211 = result211.resid
print("\nMean of residual error is {}".format(resid211.mean()))
lbtest211 = make_lbtest_df(resid211)
lbtest211.index = range(1, 11)
lbtest211
                               Statespace Model Results                           
    ==============================================================================
    Dep. Variable:                    raw   No. Observations:                 1000
    Model:               SARIMAX(2, 1, 1)   Log Likelihood               -3015.019
    Date:                Mon, 04 Jun 2018   AIC                           6040.039
    Time:                        18:09:27   BIC                           6064.577
    Sample:                    05-01-2000   HQIC                          6049.365
                             - 06-11-2000                                         
    Covariance Type:                  opg                                         
    ==============================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
    ------------------------------------------------------------------------------
    intercept      0.1022      0.005     20.028      0.000       0.092       0.112
    ar.L1         -0.0315      0.032     -0.977      0.329      -0.095       0.032
    ar.L2          0.0089      0.033      0.270      0.787      -0.056       0.074
    ma.L1         -0.9917      0.006   -163.987      0.000      -1.004      -0.980
    sigma2        24.3877      1.111     21.959      0.000      22.211      26.564
    ===================================================================================
    Ljung-Box (Q):                       33.84   Jarque-Bera (JB):                 0.26
    Prob(Q):                              0.74   Prob(JB):                         0.88
    Heteroskedasticity (H):               0.98   Skew:                             0.03
    Prob(H) (two-sided):                  0.85   Kurtosis:                         2.95
    ===================================================================================
    
    Warnings:
    [1] Covariance matrix calculated using the outer product of gradients (complex-step).
    
    Mean of residual error is -0.12554654527010825

Q p-value
1 0.000015 0.996934
2 0.000184 0.999908
3 1.281983 0.733416
4 3.151668 0.532772
5 3.285055 0.656130
6 3.398996 0.757356
7 11.631083 0.113362
8 12.315876 0.137657
9 12.860859 0.169004
10 12.871481 0.230946

平均はおよそ-0.12、Ljung-Box検定ではいずれもP値が0.05よりも大きく、帰無仮説は棄却されていません。そのため、自己相関係数がいずれも0ではないとはいえない、という検定結果になります。

このモデルは、それほど悪くはなさそうです。

(4) (0, 1, 1)でのモデル作成

最後に、元データを作成するのに使用した、(p, d, q) = (0, 1, 1)でのモデル作成を行います。

model011 =  sm.tsa.statespace.SARIMAX(train_data, order=(0, 1, 1), trend='c')
result011 = model011.fit(disp=False)
print(result011.summary())

resid011 = result011.resid
print("\nMean of residual error is {}".format(resid011.mean()))
lbtest011 = make_lbtest_df(resid011)
lbtest011.index = range(1, 11)
lbtest011
                               Statespace Model Results                           
    ==============================================================================
    Dep. Variable:                    raw   No. Observations:                 1000
    Model:               SARIMAX(0, 1, 1)   Log Likelihood               -3015.563
    Date:                Mon, 04 Jun 2018   AIC                           6037.126
    Time:                        18:09:38   BIC                           6051.850
    Sample:                    05-01-2000   HQIC                          6042.722
                             - 06-11-2000                                         
    Covariance Type:                  opg                                         
    ==============================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
    ------------------------------------------------------------------------------
    intercept      0.1000      0.002     66.342      0.000       0.097       0.103
    ma.L1         -0.9923      0.006   -169.478      0.000      -1.004      -0.981
    sigma2        24.4122      1.107     22.045      0.000      22.242      26.583
    ===================================================================================
    Ljung-Box (Q):                       34.11   Jarque-Bera (JB):                 0.21
    Prob(Q):                              0.73   Prob(JB):                         0.90
    Heteroskedasticity (H):               0.98   Skew:                             0.03
    Prob(H) (two-sided):                  0.84   Kurtosis:                         2.96
    ===================================================================================
    
    Warnings:
    [1] Covariance matrix calculated using the outer product of gradients (complex-step).
    Mean of residual error is -0.12792704758113668
Q p-value
1 0.950670 0.329549
2 1.037009 0.595410
3 2.404111 0.492870
4 4.300639 0.366845
5 4.366195 0.497983
6 4.532642 0.604989
7 12.781170 0.077623
8 13.305484 0.101762
9 13.867792 0.127106
10 13.874346 0.178801

こちらも平均はおよそ-0.12、Ljung-Box検定でのP値はいずれも0.05より大きく、帰無仮説は棄却されません。よってこのモデルの残差はホワイトノイズに近く、ある程度良いモデルであることが期待できます。

以上、4つのパラメータでモデルの作成と評価を行いました。

最初に試したモデルのように、AICが最小の場合でも、必ずしも理想的なモデルになるとは限らないことも確認できました。また、パラメータによっては収束性に問題が発生する場合があることも確認できました。

6. 選択したモデルを用いた予測

最後に、statsmodelsを使用した場合の予測を行ます。 モデルのパラメータには、総当たり法で3番目にAICが低く、モデル確認の結果も悪くなかった(p, d, q) = (2, 1, 1)を使用します。

statsmodelsのpredictionのうち、in-sample、つまりデータがある部分の予測には、one-step-predictionとdynamic predictionの2種類の方法があります。

one-step-predictionでは、すでにデータがある期間に対しては、データの予測に必要な過去のデータとしては実際のデータを使用します。一方、dynamic predictionでは、モデルによって作成されたデータを利用して次のデータを予測します。

詳しくは下記のドキュメントやチュートリアル等をご参照ください。

それでは、グラフを確認してみましょう。

pred_start = '2000-06-01'
pred_end = '2000-06-20'
focus_start = '2000-06-01'
dynamic_start = '2000-06-05'
focus_end = pred_end
focus_time = [pred_start, focus_end]

# one-step-ahead forecast
predict211 = result211.get_prediction(start=pred_start, end=pred_end)
predict211_ci = predict211.conf_int()
predict211_mean = predict211.predicted_mean

# dynamic prediction
predict211_dy = result211.get_prediction(start=pred_start, end=pred_end, dynamic=dynamic_start)
predict211_dy_ci = predict211_dy.conf_int()
predict211_dy_mean = predict211_dy.predicted_mean


fig, ax = plt.subplots(figsize=(12, 6))

df_base.raw.plot(color='black', ax=ax)
df_all.iloc[1001:].raw.plot(color='b', ax=ax, alpha=0.5, style="--")

# one-step-ahead forecast
predict211_mean.plot( color='r',ax=ax, label='One-step-ahead forecast')
ax.fill_between(predict211_ci.index,  predict211_ci.iloc[:, 0], predict211_ci.iloc[:, 1], alpha=0.1)

# dynamic prediction
predict211_dy_mean.plot( color='g',ax=ax, label='dynamic forecast')
ax.fill_between(predict211_dy_ci.index,  predict211_dy_ci.iloc[:, 0], predict211_dy_ci.iloc[:, 1], alpha=0.1)

ax.set_xlim(focus_time)
ax.set_ylim([50, 150])

f:id:brains_iwaki:20180604175648p:plain

黒の実線が、今回モデル作成に使用したデータを示し、青の破線がはじめに分けておいた、モデル作成に使用していないデータを示します。

予測は赤の実線でone-step predictionの平均値、緑の実線でdynamic predictionの平均値を示し、
95%信頼区間をそれぞれ薄い赤、薄い緑で示しました。今回はこの 2種類の信頼区間がほとんど重なっています。

一部信頼区間から外れていることを含め、この予想結果をどう評価するかは、状況によりけりといったところでしょうか。

まとめ

Pythonのstatsmodelsを用いてのモデル作りから予測までの一通りの流れを行いました。

今回は単純なデータを対象にしましたが、扱った範囲だけでも、次数を試す範囲やモデル選択、モデルの解釈と、検討事項はまだまだ多い印象を受けました。現実のデータでは、季節項やトレンド項、あるいは外的な要因が入ってきたりと、更に難易度が上がっていきそうです。

目的との兼ね合いも考慮しながら、取り入れていけるかどうか検討していきたいと思います。

ブレインズテクノロジーでは「共に成長できる仲間」を募集中です。
採用ページはこちら

参考文献および参考ウェブサイト

PythonでのARIMAモデルを使った時系列データの予測の基礎[前編]

こんにちは。ブレインズテクノロジーの岩城です。

今回は、2部構成でPythonを使用してのARIMAモデルの作成・予測の流れを整理したいと思います。

前半の本記事では、時系列データの予測でよく利用される、ARIMAモデルの基礎を簡単におさらいします。

後半の記事では実際にコードを追いながら一連の流れを確認していきますので、先にコードを読みたい場合はこちらをご参照ください。

blog.brains-tech.co.jp

主な参考文献

本記事で紹介するARIMAモデルの理論の部分は、主にこちらの教科書を参考にしています。

沖本竜義(2010), "経済・ファイナンスデータの計量時系列分析", 朝倉書店http://www.asakura.co.jp/G_12.php?isbn=ISBN978-4-254-12792-8

www.asakura.co.jp

きちんと理論を把握したい場合は、是非読んでみてください。

ARIMAとは

まずは、ARIMA(Autoregressive integrated moving average、自己回帰和分移動平均)モデルについて、大まかな内容を確認します。細かな部分は扱いませんので、きちんとした理解が必要な方はご紹介した教科書等をご参照ください。

用語は長いですが、いくつかの用語を組み合わせただけなので、一つずつみていけばそこまで難しくないと思います。

先に、ARIMA過程の定義を先程の教科書からお借りします。

定義 5.3 (ARIMA過程)
d階差分を取った系列が定常かつ反転可能なARMA(p, q)過程に従う過程は次数(p, d, q)の自己回帰和分移動平均過程もしくはARIMA(p, d, q)過程と呼ばれる.
(沖本竜義(2010), "経済・ファイナンスデータの計量時系列分析", P. 106, 朝倉書店より引用)

要はある時系列データに対して、

  • (1)d階差分をとったデータが、
  • (2)定常かつ
  • (3)反転可能な
  • (4)ARMA過程

で表すことができるような過程がARIMA過程で、これを数式化したのがARIMAモデルです。

順に一つずつ確認していきます。

(1)d階差分について

d階差分については問題ないかと思いますが、念のため簡単に説明します。

元の系列が y_0, y_1, y_2, ... , y_tだったとします。

1階差分を取った系列は \Delta y_k = y_k - y_{k-1}  (1\leq k \leq t)で構成される系列です。

2階差分は1階差分の系列から同様の差分を取った系列で、以下d階差分についても同様です。

(4)ARMA過程について ((2)と(3)は後回し)

次は、(2)と(3)は後回しにして、AR過程、MA過程、ARMA過程について説明することにします。

AR過程

AR過程(autoregressive process, 自己回帰過程)は、ある時刻のデータが、過去の時刻のデータから表現することができる過程です。p次のAR過程であるAR(p)は、

  y_t = c + \psi _1 y_{t-1} + \psi _2 y_{t-2} + ... + \psi _p y_{t-p} + \varepsilon _{t}

と表現されます。ここで、 c は定数、 \varepsilon_tは誤差項で、平均0、標準偏差σの正規分布に従います。

MA過程

MA過程(moving average process, 移動平均過程)は、ある時刻のデータが、過去の時刻での誤差項を用いて表現することができる過程で、q次のMA過程MA(q)は、

   y_t = c + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + ... + \theta_q \varepsilon_{t-q} + \varepsilon_t

と表現されます。 \varepsilon_tは先程と同様に誤差項で、各時間の誤差項が平均0、標準偏差σの正規分布に従います。ある時刻と別の時刻の式の中に、共通して存在する誤差項を持たせることにより、相関関係を表現します。

AR過程とMA過程の特徴

AR過程、MA過程ともに、注目しているデータと過去のデータとの関係を表す方法です。

AR過程では過去のデータに係数をかけたものをいくつか組み合わせて注目するデータを表現します。

MA過程では、注目しているデータと過去のデータに共通する項を持たせることで関係性を表現します。

ARMA過程

ARMA(autoregressive moving average, 自己回帰移動平均)過程は、AR過程とMA過程を組み合わせて表現される過程です。ARMA(p, q)は、

   y_t = c + \psi_1 y_{t-1} + \psi_2 y_{t-2} + ... + \psi_py_{t-p} + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + ... + \theta_q \varepsilon_{t-q} + \varepsilon_t

と表現されます。

(2)定常性

定常性は、モデルの特性の一つです。こちらについても、教科書から定義をお借りします。

定義 1.1(弱定常性) 任意のtとkに対して

   E(y_t) = \mu
   Cov(y_t, y_{t-k}) = E[(y_t - \mu)(y_{t-k} - \mu)] = \gamma_{k}

が成立する場合,過程は弱定常(weak stationarity)といわれる.

(沖本竜義(2010), "経済・ファイナンスデータの計量時系列分析", P. 8, 朝倉書店より引用)

もしよく分からない場合は、厳密な表現ではありませんが、

  • 各時刻のデータ y_tの平均値(期待値)が一定で、
  • 任意の時刻kだけ離れた y_tとy_{t-k}の相関係数がkにのみ依存し、tに依存しない。

の二つの条件を満たすデータ系列だと考えてください。
詳細は省きますが、MA過程は常に定常なので、AR過程の部分が定常性を満たせばこの条件を満たすことができます。
また、AR過程が定常である時、AR(p)はMA(∞)で書き直すことができます。

(3)反転可能性

反転可能性は、MA過程に関する条件です。反転可能である時、MA(q)をAR(∞)と書き直すことが可能です。こちらも詳細な説明はできませんが、少しだけ注釈を加えておくと、一般にMA過程は複数存在するのに対し、反転可能性を満たすMA過程は1種類しかないという特徴があります。なお、ARMA過程が反転可能である条件は、ARMAを構成するMA過程が反転可能であることを意味しています。

ARIMAモデルのまとめ

ARIMAモデルの定義と、求められる性質について整理しました。不正確な表現であることを覚悟して、ARIMAモデルを簡潔に表現しようとすると、

「ある時系列データに対して何度か差分を取ったデータがARMAモデルで表現できて、そのARMAモデルは、定常で反転可能であり、 無限に項を使えばAR項だけでもMA項だけでも表現できる」

といったところでしょうか。正確な理解が必要な場合は、やはり参考書をご確認ください。

まとめ

今回の記事では、ARIMAモデルについて、最低限必要な情報を整理しました。大まかな内容だけでも把握していただき、必要に応じて読み直したり、詳しく調べる際の足がかりにしていただれば幸いです。詳細な内容を確認したい場合は、今回何度も引用させていただいている教科書をおすすめします。

後半の記事では実際にPythonを使用して、一連のモデル作成の流れを確認していきます。こちらの記事も是非ご一読ください。 blog.brains-tech.co.jp
ブレインズテクノロジーでは「共に成長できる仲間」を募集中です。
採用ページはこちら

参考文献および参考ウェブサイト

WindowsのDockerで慣れないこと

Neuron開発チームの木村です。

Neuronは、ほとんどの場合Windowsにインストールされるため、テスト環境をWindows上に構築しています。 このテスト環境構築を効率化するため、Docker上でのWindows利用を模索し始めています。

今回は「Docker for Macを使っていた私が、WindowsでDockerを使う際に慣れないこと」を紹介します。

なお、利用しているWindowsのバージョンは
「Windows 10 Pro バージョン 1803 (Version 10.0.17134.48)」です。

以下、OSがWindowsであるコンテナを「Windows OSコンテナ」、OSがLinuxであるコンテナを「Linux OSコンテナ」と呼びます。

目次

どのDockerを使うべきか分からない

(この節は、Windows 10 Pro または Windows Server 2016 のバージョン 1709以上を前提とします。)

macでは、とりあえずDocker for Macをインストールすればよかったのですが、 WindowsでDockerを使用するには、以下の複数の方法があります。

  • Docker Toolbox
  • Docker for Windows
  • Hyper-Vコンテナ + LCOW
  • Windowsコンテナ

Docker Toolbox以外は、Windows OSコンテナもLinux OSコンテナも利用できますが、一体どれを使うべきなのでしょうか。

それぞれの特徴をまとめました。

Docker Toolbox

インストール方法

長所

  • (強いて言えば、)Hyper-Vなしで利用できる
    • 代わりに、VirtualBoxが使われる

短所

  • Windows OSコンテナを利用できない
  • 「Legacy desktop solution.」と説明されるなど、deprecated感がある

Docker for Windows

インストール方法

長所

  • docker composeなどのツール群が一緒にインストールされる

短所

  • Windows OSコンテナの実行にも、Hyper-Vによる仮想化が行われる
  • Linux OSコンテナ用のdocker daemonがMobyLinuxVMというHyper-Vインスタンス上で動作するため、二段の仮想化となり、コンテナ実行時にファイルのpermissionに引っかかるなど面倒なことが起きる可能性がある
    • (experimental機能(LCOW利用)を有効にすれば、MobyLinuxVMは介さないかもしれません)
  • Windows OSコンテナとLinux OSコンテナを同時に操作できない
    • 操作したいコンテナに合わせて、「モードを切り替える」または「experimental機能の--platformオプションを指定する」必要がある

Hyper-Vコンテナ + LCOW

インストール方法

長所

  • Windows OSコンテナとLinux OSコンテナを同時に操作できる
    • 現段階では、docker pullには--platformオプションが必要

短所

  • Windows OSコンテナの実行にも、Hyper-Vによる仮想化が行われる
  • LCOWは、現在experimental
  • 上記のインストール方法では、docker clientにパスを通したり、docker daemonのサービス化は自身で行うこととなる

Windowsコンテナ

インストール方法

長所

  • Windows OSコンテナについては、ホストOSのカーネルを使って実行されるため、Hyper-Vによる仮想化が行われない

短所

  • Windows Server 2016が必要
  • Windows OSコンテナについて、コンテナのWindowsバージョンとホストのWindowsバージョンが、ビルド番号まで一致している必要がある(バージョンが異なる場合、起動がブロックされる


どれを使えばよいかは要件次第なので、今回は、

  • Windows OSコンテナとLinux OSコンテナの両方を使いたい
  • OSをなるべく意識せずに、dockerコマンドを実行したい
  • テスト環境構築に使用するため、パフォーマンスは多少犠牲にしてもよい
  • (experimentalでも使ってみたい)

といった理由から、「Hyper-Vコンテナ + LCOW」を使うこととしました。


Windowsのdockerイメージはどれなのか

まだDocker上でWindowsが動かなかった頃は、どのdockerイメージをpullしてもLinuxのイメージであり、話は単純でした。
さて、では現在、Windowsのイメージをどう見分けてpullするのでしょうか。

結論から言うと、「これはWindows、これはLinux」といった目印は見つけられませんでした。

Windowsのイメージを探す時は、dockerHubのwinamd64microsoftをあたっていくと良いと思います(ただし、microsoftのほうはLinuxのイメージも多い)。

希望のイメージが見つからない場合は、以下のWindows Server 2016基底イメージを使ってDockerfileを書いてゆきましょう。

なお、Windows Server 2012など旧Windowsのイメージは存在しないと思われます。


DockerfileのRUN/CMDは、コマンドプロンプトで実行される

当然なのですが、DockerfileからWindowsのイメージをビルドする時、RUNCMDの内容は、Windows上で実行されます。 つまり、Linuxのshellコマンドをそのまま使うことはできません。

こればかりは、コマンドプロンプト(または、そこから呼び出すPowerShell)のコマンドに慣れてゆくしかありません…
さらに、パス区切りがややこしく、COPYなど/でよい命令と、WORKDIRなど\でなければならない命令が混在しています。

(bashのコマンドからWindowsのコマンドに変換する何かを、docker buildに挟みたい)


おわりに

今回は、dockerのインストールやdockerイメージの作成について、WindowsでDockerを使う際に慣れないことを紹介しました。
次回は、実際にコンテナを動かしていく際のことを紹介する予定です。 何事もなくすんなり使えればよいですが…

実践編に続きます。


ブレインズテクノロジーでは「共に成長できる仲間」を募集中です
採用ページはこちら

参考

Dockerコンテナで起動するElasticsearchの性能情報を取得した話

ブレインズテクノロジーの加藤です。

今回は、Dockerコンテナ上で起動するjavaアプリケーションの性能情報の取得方法についてまとめます。
例としてElasticsearchを使用しましたが、javaで起動しているアプリケーションであれば何にでも応用が効くはずです。

経緯

あるお客様の環境で、Dockerコンテナ上にElasticsearchを起動してサービスを提供しているのですが、GC系のERRORが頻発するようになりました。

原因の特定と解決策の検討のため、jstatによる性能情報を取得しようとしたのですが、ホストから見たElasticsearchのプロセスIDが異なるからか、ホスト側からは取得ができず。
仕方なくDockerコンテナに入って取得しようとしたのですが、今度はDockerコンテナにjdkが入っていないため、インストールするところから始めました。

以下はその手順です。

Dockerコンテナにログインする

Elasticsearchが起動していること、およびプロセス名を確認します。
今回のプロセス名は"brave_jepsen"のようです。

$ docker ps
CONTAINER ID        IMAGE                  COMMAND                  CREATED                  STATUS              PORTS                NAMES
3596bfd346d4        elasticsearch:latest   "/docker-entrypoint.…"   Less than a second ago   Up 2 seconds        9200/tcp, 9300/tcp   brave_jepsen

docker exec コマンドにより、elasticsearchのコンテナにログインします。

$ docker exec -it brave_jepsen bash
root@3596bfd346d4:/usr/share/elasticsearch#

ログインできました。

Dockerコンテナにjdkをインストールする

試しにDockerコンテナの中で、jstatコマンドを実行してみます。

$ jstat
bash: jstat: command not found

入ってないですね。仕方ないのでインストールします。
インストールにはapt-getコマンドを使用します。

まずはapt-getを最新化しましょう。

$ apt-get update
Ign:1 http://deb.debian.org/debian stretch InRelease
Get:2 http://deb.debian.org/debian stretch-updates InRelease [91.0 kB]
・・・(略)
Fetched 10.3 MB in 3s (3173 kB/s)
Reading package lists... Done

ちなみに環境によっては、apt-getを実行しても失敗することがあります。
その際は、環境変数「http_proxy」にプロキシサーバを設定すれば通るかもしれません。

export http_proxy=${プロキシサーバのホスト}:${プロキシサーバのポート}

続いて、apt-fileをインストールします。

$ apt-get install -y apt-file
Reading package lists... Done
Building dependency tree
Reading state information... Done
・・・(略)
Setting up apt-file (3.1.4) ...
The system-wide cache is empty. You may want to run 'apt-file update'
as root to update the cache.

apt-fileをupdateしろと言っているので、仰せのままに実行します。

$ apt-file update
Ign:1 http://deb.debian.org/debian stretch InRelease
Get:2 http://security.debian.org stretch/updates InRelease [63.0 kB]
・・・(略)
Building dependency tree
Reading state information... Done
3 packages can be upgraded. Run 'apt list --upgradable' to see them.

apt-fileのsearchコマンドにより、jstatが含まれるパッケージを見つけます。

$ apt-file search jstat
libmpj-java: /usr/bin/mpjstatus
openjdk-8-jdk-headless: /usr/lib/jvm/java-8-openjdk-amd64/bin/jstat
openjdk-8-jdk-headless: /usr/lib/jvm/java-8-openjdk-amd64/bin/jstatd
openjdk-8-jdk-headless: /usr/lib/jvm/java-8-openjdk-amd64/man/ja_JP.UTF-8/man1/jstat.1.gz
openjdk-8-jdk-headless: /usr/lib/jvm/java-8-openjdk-amd64/man/ja_JP.UTF-8/man1/jstatd.1.gz
openjdk-8-jdk-headless: /usr/lib/jvm/java-8-openjdk-amd64/man/man1/jstat.1.gz
openjdk-8-jdk-headless: /usr/lib/jvm/java-8-openjdk-amd64/man/man1/jstatd.1.gz

「openjdk-8-jdk-headless」に含まれるようです。
なので対象のパッケージをインストールします。

$ apt-get install -y openjdk-8-jdk-head
Reading package lists... Done
Building dependency tree
・・・(略)
update-alternatives: using /usr/lib/jvm/java-8-openjdk-amd64/bin/wsgen to provide /usr/bin/wsgen (wsgen) in auto mode
update-alternatives: using /usr/lib/jvm/java-8-openjdk-amd64/bin/jcmd to provide /usr/bin/jcmd (jcmd) in auto mode

jpsコマンドにより、Elasticsearchのプロセスが表示されれば成功です。
elasticsearchの場合、プロセスIDは1になるようです。

$ jps -v
1 Elasticsearch -Xms2g -Xmx2g -XX:+UseConcMarkSweepGC -XX:CMSInitiatingOccupancyFraction=75 -XX:+UseCMSInitiatingOccupancyOnly -XX:+AlwaysPreTouch -Xss1m -Djava.awt.headless=true -Dfile.encoding=UTF-8 -Djna.nosys=true -Djdk.io.permissionsUseCanonicalPath=true -Dio.netty.noUnsafe=true -Dio.netty.noKeySetOptimization=true -Dio.netty.recycler.maxCapacityPerThread=0 -Dlog4j.shutdownHookEnabled=false -Dlog4j2.disable.jmx=true -Dlog4j.skipJansi=true -XX:+HeapDumpOnOutOfMemoryError -Des.path.home=/usr/share/elasticsearch
802 Jps -Dapplication.home=/usr/lib/jvm/java-8-openjdk-amd64 -Xms8m

Javaの性能情報取得

jstatコマンドにより性能情報を取得します。
使用可能なオプションの種類や、出力項目の詳細は、下記のオラクルのページに載っています。

docs.oracle.com

今回は gccauseオプションを使用したいと思います。
コマンドの意味としては、"-h1000"で1000行ごとにヘッダーを表示。
途中の1でelasticsearchのプロセスIDを指定。
1000は、1秒ごとに1000回繰り返すことでその間の性能情報を取得することができます。

$ jstat -gccause -h1000 1 1000 > /usr/share/elasticsearch/logs/gccause.txt
$ cat /usr/share/elasticsearch/logs/gccause.txt
  S0     S1     E      O      M     CCS    YGC     YGCT    FGC    FGCT     GCT    LGCC                 GCC
100.00   0.00  98.21   0.53  93.29  83.32      2    0.115     2    0.037    0.152 CMS Final Remark     No GC
100.00   0.00  98.21   0.53  93.29  83.32      2    0.115     2    0.037    0.152 CMS Final Remark     No GC
100.00   0.00  98.21   0.53  93.29  83.32      2    0.115     2    0.037    0.152 CMS Final Remark     No GC
100.00   0.00  98.21   0.53  93.29  83.32      2    0.115     2    0.037    0.152 CMS Final Remark     No GC
100.00   0.00  98.21   0.53  93.29  83.32      2    0.115     2    0.037    0.152 CMS Final Remark     No GC
100.00   0.00  98.21   0.53  93.29  83.32      2    0.115     2    0.037    0.152 CMS Final Remark     No GC
100.00   0.00  98.21   0.53  93.29  83.32      2    0.115     2    0.037    0.152 CMS Final Remark     No GC

GCログはDockerコンテナ内に出力されるので、volume mappingでホスト側と共有しておくと良いです。

ちなみにお客様の環境ではdocker-composeを使っています。
docker-compose.ymlの設定は以下のような感じです。

(抜粋)
  volumes:
        - /home/logs/elasticsearch:/usr/share/elasticsearch/logs

この設定をしておけば、ホスト側の「/home/logs/elasticsearch」ディレクトリにgccause.txtが出力されているはずです。
gccauseだと実行時刻が出力されないため、ホスト側で下記のようなコマンドを実行するといい感じです。

$ tail -f /home/logs/elasticsearch/gccause.txt | awk '{print strftime("%Y/%m/%d %H:%M:%S"), $0; fflush()}' > /home/logs/elasticsearch/gccause_new.txt

(出力結果)
2018/02/23 15:54:50   S0     S1     E      O      M     CCS    YGC     YGCT    FGC    FGCT     GCT    LGCC                 GCC                 
2018/02/23 15:54:50  99.35   0.00  81.41   5.55  98.06  94.74     40    2.183     2    0.108    2.292 Allocation Failure   No GC 
2018/02/23 15:54:50  99.35   0.00  81.43   5.55  98.06  94.74     40    2.183     2    0.108    2.292 Allocation Failure   No GC
2018/02/23 15:54:50  99.35   0.00  81.45   5.55  98.06  94.74     40    2.183     2    0.108    2.292 Allocation Failure   No GC
2018/02/23 15:54:50  99.35   0.00  81.46   5.55  98.06  94.74     40    2.183     2    0.108    2.292 Allocation Failure   No GC
2018/02/23 15:54:50  99.35   0.00  81.48   5.55  98.06  94.74     40    2.183     2    0.108    2.292 Allocation Failure   No GC

終わりに

というわけで、無事にDockerコンテナで起動するElasticsearchの性能情報を取得できました。
冒頭でも書きましたが、javaアプリケーションであれば、Elasticsearch以外でもこの方法で取得できるはずです。

この記事が誰かの役に立てればいいと思うー。そだねー。


ブレインズテクノロジーでは「共に成長できる仲間」を募集中です
採用ページはこちら