Technology Topics by Brains

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

速報 re:invent 2018

こんにちは、ブレインズテクノロジーの佐々木です。現在AWSのイベントであるre:invent 2018に参加しています。

こちらの時間で11月28日(水)の午前中にCEOのAndy Jassyによる基調講演が行われました。 基調講演は例年新しいサービス・機能が発表される場であるため、毎年非常に注目されています。

今年もたくさんの新サービスが発表されましたので、簡単に紹介したいと思います。 なお、この記事は速報のため、内容に誤りがある場合がございます。ご了承ください。

f:id:ryotasasaki:20181129042326j:plain
会場の様子

続きを読む

JavaScriptでデータフレーム操作 - dataframe-js -

はじめに

こんにちは、ブレインズテクノロジーの佐々木です。４月に新卒で入社しました。学生時代は熱帯林を這いつくばって植物と戯れていました。

CSVからデータを読み取って、ちょっとした分析をしたい！というシーンはよくありますよね。そうした手軽なデータ操作のツールとしては、PythonのpandasやRのdplyrなどが機能も充実していて使いやすいため、人気があるようです。また、大量のデータを製品環境で扱うような場面では、hadoopやsparkなどの分散処理基盤を活用することも多いようです。ただ、これらの選択肢が便利すぎるがゆえに、フロントエンドでデータ操作に迫られたときにストレスを感じる人も多いのではないでしょうか。

ということで、今回は、ブラウザ上でも便利にデータ操作ができるJavaScriptライブラリ、dataframe-jsに入門します。

インストールと準備

といいつつ実際にブラウザで動かすのは面倒なので、今回はnode上で動かします。基本的にはブラウザでも同様に動く（と信じて）います。またもっと手軽に試したい方はJupyter notebookにIJavaScriptカーネルを入れると良いかもしれません（ここまでやると、pandas使えよ！ってツッコミが入りそうなのでやめました）。

npm install dataframe-jsでインストールできます。 前提として、以下のようなディレクトリ構成で作業しています。

.
├── iris.js
├── main.js
├── package.json

csvからデータを読み込む場合、データそのものではなくてPromiseが返ってきます。今回は以下のようなIrisクラスにメソッドを書き加えて、main.jsから呼び出す形で作業を進めて行きます。

  • iris.js
const DataFrame = require('dataframe-js').DataFrame;

module.exports = class Iris {
  constructor(params) {
    this.url = params.url;
    this.dfPromise = DataFrame.fromCSV(this.url);
  }
  
  showDF() {
    this.dfPromise.then(df => {
      df.show();
    });
  }
}
  • main.js
const Iris = require('./dataframe');

const iris = new Iris({
    url: 'https://raw.githubusercontent.com/uiuc-cse/data-fa14/gh-pages/data/iris.csv'
});

iris.showDF();

node main.jsで実行します。

ちょっとアヤメの話

それでは本題に入ります。この記事のゴールは、iris（アヤメ）データセットを用いて種ごとに花弁(petal)と萼片(sepal)のサイズの関係を調べ（て、I. versicolor, I. virginica, I. setosa各種の形態的特性とその適応的意義について考察す）ることとします。

USDA Forest Serviceによると、この三種はいずれも湖畔沿などの湿地に生息しsetosa（アラスカ）, versicolor（五大湖付近）, virginica（ニューヨークからフロリダまで）の順に寒冷地に分布しています。一般に寒い方が大気飽差が小さく（相対湿度が大きく）蒸散で水を失いにくいことが知られています。また、組織が幅広だと乾燥重量あたりの蒸散量が増えるので、virginica, versicolor, setosa の順に花弁および萼片が細長であるという仮説をたてました。

この仮説を検証するために、各３種ごとに花弁と萼片の縦横の比の平均と標準偏差を計算し、比較します。

f:id:ryotasasaki:20181009144400p:plain http://suruchifialoke.com/2016-10-13-machine-learning-tutorial-iris-classification/

データ分析

まずはデータの基本的な整形をしてみます。行のフィルターにはfilter, 列の選択にはselectが使えます。使い方はpandasというよりも、Rのdplyrやsparkのデータフレームに似ています。

  • iris.js
filterAndSelectDF() {
    this.dfPromise.then(df => {
      df
        .filter(row => row.get("species") === "versicolor")
        .select("sepal_length", "sepal_width", "species")
        .show(3);
    });
}
  • 結果
sepal_length sepal_width species
7 3.2 versicolor
6.4 3.2 versicolor
6.9 3.1 versicolor

次に、花弁縦横比（縦/幅）を計算してみます。

  • iris.js
mutateWLratio() {
    this.dfPromise.then(df => {
      df
        .map(row => row.set('sepal_wlratio', row.get('sepal_length') / row.get('sepal_width')))
        .map(row => row.set('petal_wlratio', row.get('petal_length') / row.get('petal_width')))
        .select("sepal_wlratio", "petal_wlratio", "species")
        .show(3);
    });
}
  • 結果
sepal_wlratio petal_wlratio species
1.4571 6.9999 setosa
1.6333 6.9999 setosa
1.46875 6.5 setosa

最後に、種ごとにグループ演算をして各種の形質の平均・標準偏差を求めます。groupByで指定した列に対しGroupedDataFrameオブジェクトが返されます。GroupedDataFrameオブジェクトはaggregateメソッドをもち、グループ（種）ごとに関数を適用できます。

今回は花弁・萼片それぞれに平均と分散を種ごとに計算し、一つのデータフレームにまとめています。

  • iris.js
calcSppStats() {
    this.dfPromise.then(df => {
      const groupedDF = df
      .chain(
        row => row.set('sepal_wlratio', row.get('sepal_length') / row.get('sepal_width')),
        row => row.set('petal_wlratio', row.get('petal_length') / row.get('petal_width'))
      )
      .select('sepal_wlratio', 'petal_wlratio', 'species')
      .groupBy('species');
      groupedDF
        .aggregate(group => group.stat.mean('sepal_wlratio'))
        .rename('aggregation', 'sepal_wlratio_sp_mean')
        .join(
          groupedDF
            .aggregate(group => group.stat.sd('sepal_wlratio'))
            .rename('aggregation', 'sepal_wlratio_sp_sd')
        , 'species', 'inner')
        .join(
          groupedDF
            .aggregate(group => group.stat.mean('petal_wlratio'))
            .rename('aggregation', 'petal_wlratio_sp_mean')
        , 'species', 'inner')
        .join(
          groupedDF
            .aggregate(group => group.stat.sd('petal_wlratio'))
            .rename('aggregation', 'petal_wlratio_sp_sd')
        , 'species', 'inner')
        .show(3);
    });
}
  • 結果
species sepal_wlratio_sp_mean sepal_wlratio_sp_sd petal_wlratio_sp_mean petal_wlratio_sp_sd
setosa 1.4745 0.1186 7.0779 3.1237
versic 2.1604 0.2286 3.2428 0.3124
virginica 2.2304 0.2469 2.7806 0.4073

また、group.stat.statsというメソッドを使うと各種統計量がまとめて計算できます。

分析まとめ（アヤメ）

以上の結果をテーブルにまとめました。

species 花弁縦横比（縦/幅）・平均 ± 標準偏差 萼片縦横比（縦/幅）・平均 ± 標準偏差
setosa 1.48 ± 0.12 7.10 ± 3.12
versicolor 2.16 ± 0.23 3.24 ± 0.31
virginica 2.23 ± 0.25 2.78 ± 0.41

花弁は仮説の通り、virginica, versicolor, setosa の順に細長でした。しかし、萼片については反対にsetosa, versicolor, virginicaの順に細長でした。考察として、花弁や萼片の形状（縦横比）だけでなく面積そのものや生理学的特性も考慮する必要があると考えられます。

まとめ（dataframe-js）

流石にpandasには劣りますが、基本的なデータ操作はストレスなくできた気がします。クライアント側でデータ操作がスムーズにできることで、plotly.jsなどと組み合わせて（df.toCollection()とすると簡単にオブジェクトに変換できます）1. リッチな可視化が比較的簡単に実装できるようになる、2. pandasなどを使ってデータ処理をするためだけにサーバーサイドに処理を投げる必要がなくなりデータフローがすっきりさせられる、などのメリットがあるのではないでしょうか。

追記

はてなブログでJavaScriptが実行できることを知ったので、一応試しておきました。表示ボタンを押すとhttps://raw.githubusercontent.com/uiuc-cse/data-fa14/gh-pages/data/iris.csv のデータがロードされ、一連のこれでブラウザでも同様に動く（と信じて）ことも確認できました。

結果はここに出ます。
<script src="https://cdn.rawgit.com/Gmousse/dataframe-js/master/dist/dataframe-min.js"></script>
<script src="https://ajax.googleapis.com/ajax/libs/jquery/3.3.1/jquery.min.js"></script>

<input id='show_button' type=button value='表示'>
<script type="text/javascript">
$(function() {
    $('#show_button').click(function() {
        const DataFrame = dfjs.DataFrame;
        DataFrame.fromCSV('https://raw.githubusercontent.com/uiuc-cse/data-fa14/gh-pages/data/iris.csv')
            .then(df => {
                const groupedDF = df
                    .chain(
                    row => row.set('sepal_wlratio', row.get('sepal_length') / row.get('sepal_width')),
                    row => row.set('petal_wlratio', row.get('petal_length') / row.get('petal_width'))
                    )
                    .select('sepal_wlratio', 'petal_wlratio', 'species')
                    .groupBy('species');
                const resDF = groupedDF
                    .aggregate(group => group.stat.mean('sepal_wlratio'))
                    .rename('aggregation', 'sepal_wlratio_sp_mean')
                    .join(
                        groupedDF
                        .aggregate(group => group.stat.sd('sepal_wlratio'))
                        .rename('aggregation', 'sepal_wlratio_sp_sd')
                    , 'species', 'inner')
                    .join(
                        groupedDF
                        .aggregate(group => group.stat.mean('petal_wlratio'))
                        .rename('aggregation', 'petal_wlratio_sp_mean')
                    , 'species', 'inner')
                    .join(
                        groupedDF
                        .aggregate(group => group.stat.sd('petal_wlratio'))
                        .rename('aggregation', 'petal_wlratio_sp_sd')
                    , 'species', 'inner');
                const resJSON = JSON.stringify(resDF.toCollection(), null, 2);
                $('#result').text(resJSON);
            })
    });
});
</script>


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

参考

まだ日本語対応していないAmazon LexのBotをAmazon Translateで翻訳させてみた

AWS Machine Learning Bot

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

ブレインズテクノロジーは、ImpulseやNeuronといった、企業の生産性を劇的に向上させるサービス開発に従事していますが、私達自身の生産性も向上させるため、利用するシステムの探求と刷新を繰り返しています。

最近では、アマゾンウェブサービスが提供する機械学習サービスのAmazon Lexを利用して、営業記録やスケジュールの登録を簡素化できないかと模索をしていましたが、Amazon Lexとの会話は全て英語でした。

この度、Amazon Translateが日本語に対応したことを記念して、Amazon Lexとの英会話をリアルタイムに日本語に翻訳するBotを作ってみたいと思います。

f:id:brains-tech:20180830010820p:plain:w200
アマゾンウェブサービスが提供する機械学習サービス

Slack Botがやること

Slackでスケジュールを登録したいと呟くと、Amazon Lexが何の予定なのかや日時を聞いてくれるBotです。本記事ではAmazon Lexと会話するところまでを解説しますが、これをGoogle Calendar等のカレンダーと連携させて実際に予定を入れることもできたりします。

なお、ZapierがなくてもSlackとAWS Lambdaをつなぐことは出来ますが、今回は極力プログラムを書かないことを目指すためZapierを採用しています。

f:id:brains-tech:20180826220337p:plain

Amazon Translateとは？

深層学習モデルを使用して、従来の統計ベースやルールベースの翻訳アルゴリズムよりも正確で自然な翻訳を提供するニューラル機械翻訳サービスです。

2018年4月にサービスが提供され、当初は下記6言語が英語間翻訳に対応していましたが、

  • アラビア語
  • 中国語 (簡体字)
  • フランス語
  • ドイツ語
  • ポルトガル語
  • スペイン語

AWS Summit 2018 New Yorkにて、2018年7月から日本語を含めた下記6言語にも対応することが発表されました。

  • 日本語
  • ロシア語
  • イタリア語
  • 中国語 (繁体字)
  • トルコ語
  • チェコ語

Amazon Translateの詳細はこちら

Amazon Lexとは？

音声やテキストを使用して、任意のアプリケーションに対話型インターフェイスを構築するサービスです。

Amazon Alexaと同様、音声のテキスト変換には自動音声認識 (ASR)、テキストの意図認識には自然言語理解 (NLU) という高度な深層学習機能が使われているため、チャットBotのようなアプリケーションでリアルな会話を実現することができます。しかし、残念ながら現時点の対応言語は英語のみとなっており、日本語には対応していません (2018年8月時点) 。

Amazon Lexの詳細はこちら

Slack Botを作る

それでは、Slack Botを作っていきます。AWS Lambdaのところで簡単なプログラムを書きますが、後はポチポチするだけで構築完了です。大きくハマることがなければ1〜2時間程度で終わるはずです。

構築ステップ

Amazon LexとSlackの融合

公式サイトのマニュアルに詳しい手順が記載されているので、これを参考にしながら設定をしていきます。本記事では、少しややこしい連携設定を補足します。

  • ステップ 1: Amazon Lexボットを作成する
  • ステップ 2: Slackにサインアップして Slackチームを作成する
  • ステップ 3: Slackアプリケーションを作成する
  • ステップ 4: SlackアプリケーションとAmazon Lexボットを統合する
  • ステップ 5: Slack統合を完了する
  • ステップ 6: 中間テストをする

Amazon TranslateとSlackの融合

Amazon LexがSlackに返答した英語をZapierがひろってAWS Lambdaに投げる仕組みです。こちらもGUIを使ってお手軽に設定をしていきます。

  • ステップ 7: Slackに着信Webフックを追加する
  • ステップ 8: AWS Lambda関数に付与するロールを作成する
  • ステップ 9: AWS Lambda関数を作成する
  • ステップ 10: ZapierでSlackとAWS Lambdaを統合する
  • ステップ 11: 統合テストをする

ステップ 1: Amazon Lexボットを作成する

AWSマネージドコンソールからAmazon Lexを選択します。いくつかサンプルが準備されていますので、今回はその中から[ScheduleAppointment]を選択します。後は、好きなようにサンプルをいじり、[Build]と[Publish]まで実行します。

f:id:brains-tech:20180830124454p:plain

なお、[Build]が完了すると画面右側の[Test Chatbot]でテストチャットを行うことができます。こちらが入力する言葉を変えても、それに応じた返答をし、アポイントを予約するところまで導いてくれることが確認できます。

f:id:brains-tech:20180826234133p:plain:w300

ステップ 2: SlackにサインアップしてSlackチームを作成する

既存のSlackチームを使いたい場合は、このステップをとばしてOKです。

ステップ 3: Slackアプリケーションを作成する

Slack APIコンソールを開き、[Start Buildings]を押下してアプリを作成します。 [Bot Users]と[Interactive Components]の設定をした後、 [Basic Information]のClient ID、Client Secret、Verification Tokenをメモしておきます。

f:id:brains-tech:20180827112047p:plain

ステップ 4: SlackアプリケーションとAmazon Lexボットを統合する

AWSマネージドコンソールからAmazon Lexを選択します。先ほど作成したAmazon Lexボットを選択し、[Channels]タブからSlackとの連携設定をします。Client ID、Client Secret、Verification Tokenにはステップ 3のSlack APIコンソールでメモしておいた値を入力します。

f:id:brains-tech:20180827210833p:plain

[Activate]が完了するとPostback URLとOAuth URLが表示されるので、こちらをメモしておきます。

f:id:brains-tech:20180827212342p:plain

ステップ 5: Slack統合を完了する

Slack APIコンソールを開き、[OAuth & Permissions]で権限設定をしたあと、Redirect URLsにステップ 4のAWSマネージドコンソールでメモしておいたOAuth URLを入力します。

f:id:brains-tech:20180827213800p:plain

[Interactive Components]と[Event Subscriptions]のRequest URLに、ステップ 4のAWSマネージドコンソールでメモしておいたPostback URLを入力します。

f:id:brains-tech:20180827215339p:plain
f:id:brains-tech:20180827215414p:plain

これで、連携設定は完了です！

ステップ 6: 中間テストをする

Slackにログインし、テストをしたいチャンネルにSlack APIで作ったBotユーザを招待します。 ここでは、realbot子さんが#test-channelにlexbot子さんを招待しています。

f:id:brains-tech:20180827223216p:plain

少し分かりづらいのですが、realbot子さんがこちらが入力したメッセージ。lexbot子さんがAmazon Lexが返してくれているメッセージです。

f:id:brains-tech:20180830114539p:plain

このように、lexbot子さんは英語しか話せないので、次のステップ以降でAmazon TranslateとSlackを融合させて、日本語を話してもらうようにします。

ステップ 7: Slackに着信Webフックを追加する

Slack App ディレクトリを開き、着信 Webフックを追加します。

# lexbot子さんの着信 Webフックを利用してもOKですが、今回は表示名を変えたかったため、別の着信 Webフックを追加しています。

f:id:brains-tech:20180829225619p:plain

名前やアイコンを指定したあと、Webhook URLをメモしておきます。

f:id:brains-tech:20180830125305p:plain

ステップ 8: AWS Lambda関数に付与するロールを作成する

AWSマネージドコンソールからIAMを選択します。[ロール]->[ロールの作成]を押下し、AWSサービスのLambdaを選択します。

f:id:brains-tech:20180829214638p:plain

TranslateReadOnlyポリシーを選択します。

f:id:brains-tech:20180829215132p:plain

任意のロール名を入力しロールを作成します。

f:id:brains-tech:20180829215301p:plain

ステップ 9: AWS Lambda関数を作成する

AWSマネージドコンソールからAWS Lambdaを選択します。[関数の作成]を選択して、関数を一から作成します。 Zapierを利用する場合、関数の名前は必ずzapier_evalとしてください。ステップ 8で作成したロールを選択して関数を作成します。

f:id:brains-tech:20180830130047p:plain

関数コードを入れて保存します。

import json
import boto3
import unicodedata
from urllib.parse import parse_qs
import urllib.request

def chk_lang(string):
    
    for char in string:
        name = unicodedata.name(char)
        if 'CJK UNIFIED' in name or 'HIRAGANA' in name or 'KATAKANA' in name:
            return True
    return False
    
def lambda_handler(event, context):
    
    plain_text = event['plain']
    lang = chk_lang(plain_text)

    if lang: return {'statusCode': 101}

    translate = boto3.client(service_name = 'translate', region_name = 'us-east-1', use_ssl = True)
        
    response = translate.translate_text(Text = plain_text, 
                SourceLanguageCode = 'en', TargetLanguageCode = 'ja')

    translated_text = response['TranslatedText']
    
    hook = 'https://hooks.slack.com/services/xxxxx/xxxxx' #ステップ 7で確認したWebhook URLを入力
    headers = {'Content-Type' : 'application/json'}
    method = 'POST'
    
    json_data = {'text' : '`(ja) ' + translated_text + '`'}
    json_data = json.dumps(json_data).encode('utf-8')

    request = urllib.request.Request(hook,
                data = json_data,
                method = method, headers=headers)

    urllib.request.urlopen(request)
    
    return translated_text

ステップ 10: ZapierでSlackとAWS Lambdaを統合する

Zapier Exploreを開き、[Make a Zap!]を選択します。

# Zapierは様々なサービスを無料で連携させることができますが、AWS Lambdaを連携させる場合は有償となります。

Triggerは、SlackにNew Message Posted to Channelがあった場合とします。

f:id:brains-tech:20180829234023p:plain

Channelを指定し、Trigger for Bot Messages?をnoからyesに変更します。

f:id:brains-tech:20180829235032p:plain

今回は、Invoke Functionを選択します。返り値を受け取ってあれこれしたい場合はRun Codeがいいでしょう。

f:id:brains-tech:20180829235959p:plain

Functionにはステップ 9で作成したzapier_evalを選択、Argumentsでeventオブジェクトに渡すデータとしてSlackのTextを指定します。

f:id:brains-tech:20180830001744p:plain

ZapのステータスをONにしたら、全ての設定が終了です。ふぅ

ステップ 11: 統合テストをする

Slackにログインし、realbot子さんがlexbot子さんに話しかけます。すると、tranbot子さんが英語を日本語に翻訳してくれました！

f:id:brains-tech:20180830131538p:plain

おわりに

トリガーとなる英語を適切に入力しないとAmazon Lexは反応してくれませんが、Hi (やあ) やHello (こんにちは) といった簡単な言葉をトリガーにしておけば、

  1. 英語で入力
  2. 英語で返答 (by Amazon Lex)
  3. 日本語に翻訳 (by Amazon Translate)

から

  1. 日本語で入力
  2. 英語に翻訳 (by Amazon Translate)
  3. 英語で返答 (by Amazon Lex)
  4. 日本語に翻訳 (by Amazon Translate)

とすることもできそうです。

それが実現できると、InputとOutputの言語はAmazon Translateが対応する12言語を自由に組み合わせることができるようになります。 Amazon Lexが英語にしか対応していなくても、多言語コミュニケーションが可能になる訳ですね。

余談

無茶なキーワードを入れたら返ってきた返答です。普段は敬語のtranbot子さんですが、たまに不機嫌になるのでご注意を (笑)

f:id:brains-tech:20180830004242p:plain


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

WindowsのDockerで慣れないこと（実践編）

Docker Windows

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-nanoserverはpowershellをシェルに指定しています。

コンテナ実行編

[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のグラフを描画してみた

Angular Plotly

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

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

さて今回は、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を使った時系列データのクラスタリング

Machine Learning Python 時系列

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

今回は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モデルを使った時系列データの予測の基礎[後編]

Python 時系列

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

こちらの記事は、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を用いてのモデル作りから予測までの一通りの流れを行いました。

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

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

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

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