CUBE SUGAR CONTAINER

技術系のこと書きます。

ピクセラ PIX-MT100 を iPad から使ってみる

外出先でパソコンからインターネットを使いたいときがある。 そんなときのために、普段はピクセラの PIX-MT100 という LTE 対応 USB ドングルに MVNO の SIM カードを入れて持ち歩いている。

ピクセラ LTE対応USBドングル ホワイト  PIX-MT100

ピクセラ LTE対応USBドングル ホワイト PIX-MT100

  • 発売日: 2016/06/17
  • メディア: Personal Computers

この製品は Mac などのパソコンに USB 端子でつなぐと、有線の Ethernet デバイスとして認識する。 つまり、USB Ethernet アダプタを挿して LAN ケーブルをつないだような状態になる。 また、使う上でとくに専用のドライバを入れる必要もない。 テザリングに比べると無線 LAN チャネルの混雑状況に関係なく安定した通信ができるので意外と重宝している。

で、今回は PIX-MT100 が iPad から使えるのかが気になった。 製品のサポートページを見ると、一応 iPhone / iPad は iOS 8 以降でサポートされているようだ。

www.pixela.co.jp

とはいえ、事前に Web を軽く調べても、実際に検証している人が見当たらなかったので自分で試してみることにした。

確認に使った iPad の環境は次のとおり。

f:id:momijiame:20200501212844p:plain
検証に使った iPad の情報

結論としては、ちゃんと使えることが分かった。 USB-A to C アダプタ経由で挿すと、しばらくして "4G Modem" という名前で Ethernet デバイスを認識する。 以下の写真を見ると分かるとおり、PIX-MT100 のステータスランプも正常を示す青色が点灯する。

f:id:momijiame:20200501193752j:plain
PIX-MT100 を iPad OS が認識した状態

設定情報を確認すると、ちゃんと IPv4 / IPv6 の設定が配布されていることがわかる。

f:id:momijiame:20200501193802j:plain
ネットワークの設定情報

もちろん、無線 LAN などを切った状態でちゃんとインターネットが使えることも確認できた。

実際にこの状態で使うかは別として、iPad でも PIX-MT100 を使うことはできるようだ。 いじょう。

Word2Vec 形式のファイルフォーマットについて

Word2Vec では、Skip-gram や CBOW といったタスクを学習させたニューラルネットワークの隠れ層の重みを使って単語を特徴ベクトルにエンコードする。 つまり、Word2Vec で成果物として得られるのは、コーパスの各単語に対応する特徴ベクトルになる。 今回は、単語の特徴ベクトルを永続化するために使われる、Word2Vec 形式とか呼ばれているファイルフォーマットについて調べたので書いてみる。

使った環境は次のとおり。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G4032
$ python -V
Python 3.7.7

下準備

はじめに、学習済みモデルを取得するのに使うコマンドをインストールしておく。

$ brew install gzip wget

また、動作確認に使うための Python パッケージとして gensim をインストールしておく。

$ pip install gensim

Word2Vec 形式 (バイナリ)

Word2Vec 形式と呼ばれているフォーマットには、テキスト形式とバイナリ形式の 2 種類がある。 はじめに、バイナリ形式から紹介する。 たとえば、本家 Google が配布している学習済み Word2Vec モデルはバイナリ形式になっている。

code.google.com

上記から「GoogleNews-vectors-negative300.bin.gz」というファイルをダウンロードしてこよう。 ただし、圧縮された状態で 1.5GB あるので結構時間がかかる。

ダウンロードできたら Python の REPL を起動する。

$ python

Python の gzip モジュール経由でダウンロードしたファイルを開く。

>>> import gzip
>>> f = gzip.open('GoogleNews-vectors-negative300.bin.gz', mode='rb')

バイナリ形式といっても、最初の一行は ASCII 文字列で記述されている。 ここには、スペースで区切られたコーパスの単語 (語彙) 数と、単語の特徴ベクトルの次元数が入る。 今回の学習済みモデルでは語彙が 300 万語、特徴ベクトルが 300 次元とわかる。

>>> header = f.readline()
>>> header
b'3000000 300\n'

この数字はバイナリの読み込みに使うので数値に変換して変数に入れておこう。

>>> words, veclen = [int(x) for x in header.strip().split()]
>>> words, veclen
(3000000, 300)

以降のデータには、単語と特徴ベクトルがスペース区切りで記録される。 はじめに、スペースまでを単語として読み込むための関数を用意しよう。

>>> def readword(f):
...    """単語を読み込むための関数"""
...     word = []
...     while True:
...         c = f.read(1)
...         if c == b' ':
...             break
...         if c == b'':
...             raise EOFError()
...         if c != b'\n':
...             word.append(c)
...     return b''.join(word)
... 

ファイルオブジェクトを読むと、単語として '</s>' が得られた。 いや、これ HTML のタグだから、本来は除外すべき単語だと思うんだけどね。

>>> readword(f)
b'</s>'

単語を取り出せたので、次は対応する特徴ベクトルを読み出す。 特徴ベクトルでは、それぞれの次元が倍精度浮動小数点 (32 ビット) として記録されている。 そこで、まずは読み出すべき特徴ベクトルのバイト数を計算しておこう。

>>> import numpy as np
>>> np.dtype(np.float32).itemsize
4
>>> bin_weights_len = np.dtype(np.float32).itemsize * veclen
>>> bin_weights_len
1200

今回であれば 300 次元 x 4 バイト = 1,200 バイトとなった。

計算した特徴ベクトルのバイト数を、ファイルオブジェクトから読み出す。

>>> bin_weights = f.read(bin_weights_len)

読みだしたバイト列を、NumPy を使って配列に変換しよう。 これが '</s>' という単語に対応する特徴ベクトルということ。

>>> word_weights = np.frombuffer(bin_weights, dtype=np.float32)
>>> word_weights
array([ 1.1291504e-03, -8.9645386e-04,  3.1852722e-04,  1.5335083e-03,
        1.1062622e-03, -1.4038086e-03, -3.0517578e-05, -4.1961670e-04,
       -5.7601929e-04,  1.0757446e-03, -1.0223389e-03, -6.1798096e-04,
...(snip)...
       -9.8419189e-04, -5.4931641e-04, -1.5487671e-03,  1.3732910e-03,
       -6.0796738e-05, -8.2397461e-04,  1.3275146e-03,  1.1596680e-03,
        5.6838989e-04, -1.5640259e-03, -1.2302399e-04, -8.6307526e-05],
      dtype=float32)

以降は同様に単語と特徴ベクトルが交互に記録されている。 たとえば、もう一単語読み出すと 'in' という単語が得られた。 あとは、語彙の数だけ繰り返せば良いことがわかる。

>>> readword(f)
b'in'

念のため gensim を使って答え合わせをしよう。 学習済みモデルを読み込む。 この処理には数分単位で時間がかかる。

>>> import gensim
>>> model = gensim.models.KeyedVectors.load_word2vec_format('GoogleNews-vectors-negative300.bin.gz', binary=True)

読み込んだモデルを使って '</s>' に対応する特徴ベクトルを得る。

>>> model.wv['</s>']
array([ 1.1291504e-03, -8.9645386e-04,  3.1852722e-04,  1.5335083e-03,
        1.1062622e-03, -1.4038086e-03, -3.0517578e-05, -4.1961670e-04,
       -5.7601929e-04,  1.0757446e-03, -1.0223389e-03, -6.1798096e-04,
...(snip)...
       -9.8419189e-04, -5.4931641e-04, -1.5487671e-03,  1.3732910e-03,
       -6.0796738e-05, -8.2397461e-04,  1.3275146e-03,  1.1596680e-03,
        5.6838989e-04, -1.5640259e-03, -1.2302399e-04, -8.6307526e-05],
      dtype=float32)

すると、先ほど手動で読みだした数値と一致していることがわかる。

Word2Vec 形式 (テキスト)

続いてはテキスト形式について紹介する。 ただし、バイナリ形式に比べると非常に単純なので説明は一瞬でおわる。 なお、テキスト形式の Word2Vec モデルとして、あまり良いものが見当たらなかったので Facebook の fastText を使うことにした。

次のコマンドを使って配布されているモデルの先頭 3 行を読み出そう。

$ wget -qO - https://dl.fbaipublicfiles.com/fasttext/vectors-crawl/cc.en.300.vec.gz | gzip -d | head -n 3
2000000 300
, 0.1250 -0.1079 0.0245 -0.2529 0.1057 -0.0184 0.1177 -0.0701 -0.0401 -0.0080 0.0772 -0.0226 0.0893 -0.0487 -0.0897 -0.0835 0.0200 0.0273 -0.0194 0.0964 0.0875 0.0098 0.0453 0.0155 0.1462 0.0225 0.0448 0.0137 0.0570 0.1764 -0.1072 -0.0826 0.0173 0.1090 0.0207 -0.1271 0.2445 0.0375 -0.0209 -0.0445 0.0540 0.1282 0.0437 0.0588 0.0984 0.0539 0.0004 0.1290 0.0242 -0.0120 -0.0480 0.0346 -0.0664 -0.0330 -0.0625 -0.0708 -0.0579 0.1738 0.4448 0.0370 -0.1001 -0.0032 0.0359 -0.0685 -0.0361 0.0070 0.1316 -0.0945 -0.0610 0.0178 -0.0763 -0.0192 0.0033 0.0056 0.1878 -0.0754 -0.0095 0.0446 -0.0588 0.0244 -0.0251 -0.0493 0.0308 -0.0359 -0.1884 -0.0988 0.1887 0.0459 -0.0816 -0.1524 -0.0375 -0.0692 0.0427 -0.0471 -0.0086 -0.2190 -0.0064 0.0877 -0.0074 -0.1400 -0.0156 0.0161 0.1040 -0.1445 -0.0719 -0.0144 -0.0293 -0.0126 0.0619 -0.0373 -0.1471 -0.2552 -0.0685 0.2892 -0.0275 0.0436 0.0311 0.0249 0.0142 0.0403 0.1729 0.0023 -0.0255 -0.0212 0.0701 -0.0727 0.0279 0.1151 -0.0394 -0.0962 -0.0598 -0.0459 -0.0326 -0.2317 0.0945 0.0110 -0.2511 0.1087 -0.0699 0.0359 0.0208 -0.0536 0.0478 0.0178 0.0095 0.0354 -0.7726 -0.0790 0.0472 0.0584 -0.1013 0.0448 -0.1202 0.0376 0.0510 -0.0616 0.4321 0.0179 0.0263 0.0271 0.0473 -0.0951 -0.2261 0.0261 0.0262 -0.0235 -0.0369 -0.1655 -0.0697 0.0122 -0.0303 0.0427 0.0787 -0.0360 0.0206 -0.0068 0.1257 0.0447 -0.0776 -0.1122 -0.0291 0.4654 0.1010 0.4440 0.0095 0.1312 0.0766 0.0873 -0.0878 -0.0296 0.0046 0.0416 -0.0134 0.0571 -0.0109 -0.0655 0.0082 -0.0563 -0.0830 -0.0550 -0.0688 0.0091 -0.0677 -0.1001 0.0200 -0.0979 0.1134 -0.0188 0.0136 0.0782 0.0207 0.0133 -0.0492 -0.0139 0.0123 0.0360 0.1249 0.0503 0.0015 0.1246 -0.0897 -0.0121 -0.0182 0.2245 -0.0313 -0.1596 0.0073 -0.0772 -0.0830 -0.0716 0.0112 0.0218 0.1245 -0.0361 0.0312 0.0652 0.0560 0.0670 0.0709 -0.0248 0.0449 -0.1296 0.1408 -0.0359 1.1585 0.0027 0.0185 0.0549 0.1367 0.2742 -0.0472 -0.0414 -0.0548 -0.0911 -0.0015 -0.1613 0.0179 -0.0040 0.0610 0.0559 0.1138 0.2978 -0.1511 -0.0079 -0.0838 0.0296 -0.1041 0.1627 -0.0670 0.0504 -0.0420 -0.0020 0.1840 0.0596 0.0448 0.0989 -0.2157 -0.0117 0.2142 -0.1672 -0.0444 0.2045 -0.4620 -0.0482 0.0688 -0.0304 0.0478 0.1583 0.0920 0.0949 0.0650 -0.0398 -0.1376 -0.0436 0.0578 0.0188 0.0148 0.2305 -0.0696 -0.0215
the -0.0517 0.0740 -0.0131 0.0447 -0.0343 0.0212 0.0069 -0.0163 -0.0181 -0.0020 -0.1021 0.0059 0.0257 -0.0026 -0.0586 -0.0378 0.0163 0.0146 -0.0088 -0.0176 -0.0085 -0.0078 -0.0183 0.0088 0.0013 -0.0938 0.0139 0.0149 -0.0394 -0.0294 0.0094 -0.0252 -0.0104 -0.2213 -0.0229 -0.0089 -0.0322 0.0822 0.0021 0.0282 0.0072 -0.0091 -0.0352 -0.0178 -0.0706 0.0630 -0.0092 -0.0223 -0.0056 0.0515 -0.0307 0.0436 -0.0110 -0.0555 0.0089 -0.0673 0.0105 0.0574 0.0099 -0.0283 0.0470 0.0053 0.0030 0.0007 0.0443 0.0069 -0.0334 0.0091 -0.0076 0.0066 0.0917 0.0311 0.0543 0.0282 -0.0200 -0.0334 0.0053 0.0364 0.2249 0.0928 -0.0123 0.0086 -0.0599 0.0676 0.0402 0.0012 0.0464 -0.0437 0.0059 0.0917 -0.0412 -0.0151 -0.0231 0.0095 0.0588 0.0279 0.0647 -0.0568 -0.0130 0.0474 0.0354 -0.0121 -0.0077 -0.1306 0.0134 -0.0506 0.0111 0.0119 -0.0221 0.0394 0.0221 0.0245 0.0039 0.1146 0.0228 -0.0468 -0.0459 -0.0189 0.0076 -0.0302 -0.0348 -0.0289 -0.0399 0.0245 -0.0102 0.0578 -0.0388 -0.0117 -0.0304 0.2466 -0.0111 0.0356 0.0046 0.2094 -0.1022 0.0336 0.0688 -0.0708 0.0269 -0.0423 0.0077 -0.0267 0.0072 0.0035 0.0351 -0.0063 -0.4462 0.0105 -0.0120 -0.0449 -0.1696 0.0504 0.0930 -0.0044 -0.0041 0.0322 0.2026 0.0613 -0.0295 0.0228 -0.0190 0.0173 0.1480 -0.0175 -0.0125 0.0687 0.0333 -0.0303 0.0428 0.0051 0.0228 0.0104 0.0731 0.0079 -0.0051 0.0543 -0.0325 0.0512 0.0288 -0.0586 -0.0000 0.0493 0.0166 -0.0143 0.0359 0.0543 -0.0005 -0.0589 0.0162 -0.0222 -0.0199 0.0235 -0.0678 0.0179 0.0033 0.0114 0.0473 -0.0443 0.0323 0.0195 -0.0647 0.3388 0.0699 -0.0215 -0.0244 -0.0034 -0.0034 -0.0622 0.0123 0.0375 -0.0197 0.0241 -0.0877 0.0201 -0.0061 -0.0256 -0.0191 -0.0264 0.0190 -0.0422 0.0251 0.0825 -0.0096 0.1288 0.0621 0.0538 0.0189 0.0422 0.1803 -0.0010 -0.0328 -0.0559 -0.0157 0.0490 0.0352 -0.0417 0.0159 -0.0766 -0.0657 0.0497 0.0102 0.1471 -0.0711 -0.1469 0.4737 -0.0169 -0.0051 0.0159 0.0550 -0.0634 -0.0210 0.0122 0.0269 0.0060 0.0664 0.0106 -0.0705 -0.0207 -0.0784 -0.0291 -0.0283 -0.1568 -0.0393 0.0050 0.0204 -0.0026 0.0436 0.0279 -0.0393 0.0367 -0.0042 -0.0156 -0.0733 -0.1637 0.0652 -0.0062 -0.0650 -0.1984 -0.0411 -0.1534 0.0020 0.0133 -0.2364 -0.0528 -0.0042 -0.0447 0.0112 -0.0332 -0.0550 0.0013 0.0169 -0.0439 -0.0578 0.0223 -0.0777 -0.0432 -0.0251 0.2370 0.0004 -0.0042

先頭の行は、バイナリ形式と同じでコーパスの語彙数と特徴ベクトルの次元数になっている。 そして、以降は各行が単語と特徴ベクトルを表している。 最初の列が単語で、以降はスペース区切りで特徴ベクトルが浮動小数点の文字列として記録されている。

テキスト形式はとても単純なので、答え合わせをする必要もないかな。 ちなみに、単語の特徴ベクトルはモデルや配布元によってフォーマットがかなりバラバラ。 ここで紹介した Word2Vec 形式は、たくさんあるフォーマットのひとつに過ぎない。

いじょう。

ゼロから作るDeep Learning ❷ ―自然言語処理編

ゼロから作るDeep Learning ❷ ―自然言語処理編

  • 作者:斎藤 康毅
  • 発売日: 2018/07/21
  • メディア: 単行本(ソフトカバー)

NAS を買ったら両親に孫の動画を見せやすくなった話

子どもが生まれると、必然的に動画を撮影する機会が増える。 今回は、子どもを撮影した動画を保存するために NAS (Network Attached Storage) を買ったら、副次的な効果として遠隔にいる両親に孫の動画を見せやすくなって良かったという話について。

TL; DR

  • 最近の売ってる NAS は使いやすい
  • 家の外からコンテンツを閲覧するのも簡単
  • NAS の選び方について

背景

子どもが生まれる前後に、きっとこれから動画をたくさん撮るだろうと思ってビデオカメラを購入した。

そして、ビデオカメラを購入すると、今度は撮影した動画を保存しておく場所として NAS がほしくなった。 なぜなら、撮影した動画を microSD カードのような外部記録媒体に入れたままだと、次のような問題点がある。

  • 撮影した動画を鑑賞しにくい
  • 外部記録媒体が故障するリスクは比較的高い

かといって、保存する先としてクラウドストレージを選ぶと、中長期的に見て月額が大きな出費になる。 同時に、最近の製品としての NAS というものを体験してみたい気持ちもあった *1

NAS の提供する機能

はじめに、NAS の提供する主だった機能について整理しておく。

ネットワークファイルシステム

これは、NAS の提供する最も基本的な機能のひとつ。 たとえば SMB (Server Message Block) や NFS (Network File System) といったプロトコルのサーバになる。 家の中でパソコンなどから、NAS のボリュームをネットワークごしにマウントしてファイルを操作する。

リモートアクセス

この機能も、最近の NAS では標準的に提供されている。 ようするに、家の外からでも NAS に保存したファイルにアクセスするための機能。

一昔前なら Dynamic DNS やルータの Destination NAT の設定など、自分で色々と頑張る必要があったけど今は昔。 今どきのものは、そういった手間が不要になっている。 たとえば専用のアプリや Web ページに、管理用 UI から自分で設定したアカウントを入力するだけで使えます、という感じ。

今回のいわば本題である遠隔地にいる両親にも、iPad にアプリを入れてもらってコンテンツを閲覧できるようにした。

ストレージの冗長化

ほとんどの NAS は、HDD (Hard Disk Drive) や SSD (Solid State Drive) といったストレージを複数搭載できる。 そのため、RAID (Redundant Arrays of Inexpensive Disks) を組むことでストレージの物理的な故障に対する冗長性を確保できる。 RAID の方式に関しては、とりあえずストレージが 2 本積めるものなら RAID1 さえ選んでおけば良い。

その他

その他にも、メーカーやモデルによっては NAS という枠にとらわれず多種多様な機能を提供している。 たとえば、アドオンで機能を追加することで仮想マシンや Docker コンテナが動いたりするようなものまである。 もはや、管理のしやすい Linux サーバという感じ。

NAS の選び方について

ここからは、NAS を買うために見るべきポイントについて書いていく。

スマホ・タブレット向けアプリの有無

今どきの NAS は、たいていリモートアクセスのできるスマホ・タブレット向けのアプリが用意されている。 ただし、以下の点については、あらかじめ確認しておいた方が良い。

  • 自分が求めているユースケースに合致するか
  • 使うまでに必要なステップが複雑でないか

ストレージに関して

NAS には、あらかじめストレージが内蔵されているものと、自分で買ってきて入れるものがある。 日本国内のメーカーが作っている製品や、ライトユーザー向けの製品ほどストレージを内蔵している傾向がある。 ストレージを内蔵しているモデルは手間がかからない一方で、ストレージの種類が自分で選べなかったり自分で同じモデルを買うよりも割高な傾向がある。

なお、ストレージには NAS で使うことを想定して作られた専用のシリーズがある。 専用のモデルは、一般的なモデルよりも信頼性が高かったり、ファームウェアにチューニングが施されている。 たとえば Western Digital であれば Red シリーズ、Seagate なら IronWolf シリーズがそれに当たる。

国内テレビ・レコーダーとの連携

これが一番やっかいで分かりにくいポイント。 せっかく大きなストレージがあるならと、国内のテレビ・レコーダーで録画した番組を NAS に保存して再生したいというニーズはあるだろう。 そのためには、NAS が DLNA (Digital Living Network Alliance) の DTCP-IP (Digital Transmission Content Protection over Internet Protocol) という規格に対応している必要がある。 これは、著作権保護の機能が働いた状態でコンテンツをやり取りするための規格だけど、注意点が山盛りにある。 結論から先に述べると、このニーズを満たしたいときはあきらめて「テレビ・レコーダーに NAS の機能が付与されたもの」を買った方が良いと思われる。 例えば、次のようなもの。

相性問題

DTCP-IP はサーバとクライアントの相性問題が強くある。 たとえば特定のレコーダーと NAS ではうまく記録できない、NAS に記録した内容を特定のアプリでは再生できない、というもの。 そのため、実際に使ってみたりレビューで確認できる環境でないと記録・再生できるかはわからないというレベルにある。

一度入れたら取り出せない?!

ほとんどの製品は、NAS に記録したら最後、別の製品には録画した番組を移動できない。 そのため、その NAS が故障したら記録した番組はあきらめる他ない、という状態になる。 テレビ・レコーダーに NAS の機能がついたものであれば、ブルーレイに焼くような選択肢もある。 また、一部の国内製品は同じシリーズ間でのみ番組を転送できるものがある。

国外製品は基本的にアドオン機能で対応する

DTCP-IP は、基本的に日本国内の製品でのみ使われている規格となる。 そのため、国外の製品には基本的に標準では載っていないと考えた方が良い。 対応したい場合には、アドオン機能を使って有料のサードパーティ製のメディアサーバーの機能を追加することになる。

もちろん、機能を追加しても前述した問題は存在する。 また、対応しているのが特定の機種に限られるような場合もあって、なんというか地雷だらけ。 あきらめてテレビ・レコーダーを買った方が良い、という気持ちになる。

アドオン機能の有無

前述したとおり、国外の製品では特に、アドオン機能を使って機能を NAS に追加できるものが多い。 サードパーティーのメーカーがアドオンを配布するストアなども整備されていて、エコシステムが構築されている。 アドオンで追加できる機能が豊富にあると、NAS という用途に限らない使い方が楽しめる。

管理用 UI の使いやすさ

たいていの NAS は、管理用の Web UI で全て設定が完結する。 この点は自前のサーバマシンに Linux を入れてネットワークファイルシステムを構築するのとは比べ物にならないほど楽。 最近は台湾の NAS 専業メーカー (QNAP や Synology) が、この管理用 UI の分かりやすさという点で一歩リードしているようだ。

パフォーマンス

一般的なパソコンと比べると、お世辞にも NAS に積んである CPU やメモリのパフォーマンスは優れているとは言えない。 たとえば、エントリーグレードの CPU はシングルコアでメモリは 512MB とかになっている。 意外とここらへんは実際に使っているときのパフォーマンスに影響してくる。 基本的に上位のグレードを使うほどスループットなどのパフォーマンスに悩まされることは少ないと考えて良い。

主要メーカーと特徴など

  • Synology

    • 台湾の NAS 専業メーカー
    • 主な製品はストレージ非内蔵タイプ
    • DTCP-IP に対応するときはアドオン機能で対応する
      • あまり期待しない方が良い
      • DiXiM Media Server が使えるのは DS218j だけ
  • QNAP

    • 台湾の NAS 専業メーカー
    • 主な製品はストレージ非内蔵タイプ
    • DTCP-IP に対応するときはアドオン機能で対応する
      • あまり期待しない方が良い
  • NETGEAR

    • ネットワーク機器を手広く手掛けるアメリカのメーカー
    • 主な製品はストレージ非内蔵タイプ
    • DTCP-IP には対応していない
  • I-O DATA

    • IT 機器を手広く手掛ける国内メーカー
    • 主な製品はストレージ内蔵タイプ
    • DTCP-IP には標準で対応しているモデルがある
  • バッファロー

    • IT 機器を手広く手掛ける国内メーカー
    • 主な製品はストレージ内蔵タイプ
    • DTCP-IP には標準で対応しているモデルがある
  • (メーカーではないけど) テレビ・レコーダーに NAS 機能がついたもの (Panasonic ほか)

    • 基本的にストレージ内蔵タイプ
    • DTCP-IP 対応を重視するならこれ

購入した NAS について

我が家では Synology の DS218+ という製品を購入した。 なお、DTCP-IP については、前述したとおり対応できても中途半端なのであきらめている。

Synology の NAS では、リモートアクセスに QuickConnect という機能が用意されている。 これは、専用のアプリ (DS File や DS Video など) や Web ページから NAS にアクセスできるというもの。 もちろん、ルータなどに特に設定をしなくても NAT 越えできる。 設定も NAS の管理用 Web UI で接続に使う識別子を確保して、あとはユーザのアカウントを用意するだけで使える。

この製品はストレージが非内蔵タイプとなっている。 うちでは Amazon のセールで購入した Western Digital の外付けハードディスクを殻割りして取り出した HDD を使っている。 この HDD については以下を参照のこと。

blog.amedama.jp

まとめ

うちでは、両親の iPad にアプリを入れて、上記の環境に接続してもらっている。 NAS のアクセス状況を見ると、最近は毎日のように実家で孫の動画を楽しんでいるようだ。 両親いわく COVID-19 もあって物理的に会うことが難しいので余計に嬉しいとのことだった。

いじょう。

*1:以前は自宅のサーバに Linux を入れて運用していた

Python: Keras で Convolutional AutoEncoder を書いてみる

以前に Keras で AutoEncoder を実装するエントリを書いた。 このときは AutoEncoder を構成する Neural Network のアーキテクチャとして単純な全結合層から成る MLP (Multi Layer Perceptron) を使っている。

blog.amedama.jp

一方で、データとして画像を扱う場合にはアーキテクチャとして CNN (Convolutional Neural Network) が使われることも多い。 そこで、今回は CNN をアーキテクチャとして採用した Convolutional AutoEncoder を書いてみた。

使った環境は次のとおり。 CNN は MLP に比べると計算量が大きいので GPU もしくは TPU が使える環境を用意した方が良い。

$ python -V
Python 3.6.9
$ uname -a
Linux b5244776fd7d 4.19.104+ #1 SMP Wed Feb 19 05:26:34 PST 2020 x86_64 x86_64 x86_64 GNU/Linux
$ pip list | egrep -i "(keras |tensorflow-gpu )"
Keras                    2.3.1          
tensorflow-gpu           2.1.0
$ nvidia-smi
Thu Apr 16 09:25:07 2020       
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 440.64.00    Driver Version: 418.67       CUDA Version: 10.1     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|===============================+======================+======================|
|   0  Tesla P100-PCIE...  Off  | 00000000:00:04.0 Off |                    0 |
| N/A   54C    P0    39W / 250W |   4685MiB / 16280MiB |      0%      Default |
+-------------------------------+----------------------+----------------------+
                                                                               
+-----------------------------------------------------------------------------+
| Processes:                                                       GPU Memory |
|  GPU       PID   Type   Process name                             Usage      |
|=============================================================================|
+-----------------------------------------------------------------------------+

下準備

はじめに、必要なパッケージをインストールしておく。

$ pip install keras tensorflow-gpu matplotlib

Convolutional AutoEncoder を Keras の Sequential API で実装する

以下のサンプルコードでは、Keras の Sequential API を使って Convolutional AutoEncoder を実装している。 ポイントは Conv2D 層のパディングで "same" を指定しないと次元をうまく合わせることが難しい。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
from keras import layers
from keras import models
from keras import callbacks
from keras import backend as K
from keras.datasets import mnist
from matplotlib import pyplot as plt
from matplotlib import cm


def main():
    # MNIST データセットを読み込む
    (x_train, train), (x_test, y_test) = mnist.load_data()
    image_height, image_width = 28, 28

    # バックエンドに依存したチャネルの位置を調整する
    if K.image_data_format() == 'channels_last':
        x_train = x_train.reshape(x_train.shape[0],
                                  image_height, image_width, 1)
        x_test = x_test.reshape(x_test.shape[0],
                                image_height, image_width, 1)
        input_shape = (image_height, image_width, 1)
    else:
        x_train = x_train.reshape(x_train.shape[0],
                                  1, image_height, image_width)
        x_test = x_test.reshape(x_test.shape[0],
                                1, image_height, image_width)
        input_shape = (1, image_height, image_width)

    # Min-Max Normalization (0. ~ 1. の範囲に値を収める)
    x_train = x_train.astype('float32')
    x_test = x_test.astype('float32')
    x_train = (x_train - x_train.min()) / (x_train.max() - x_train.min())
    x_test = (x_test - x_test.min()) / (x_test.max() - x_test.min())

    # 畳み込み演算を用いた AutoEncoder のネットワーク (Sequential API)
    model = models.Sequential()
    # 28 x 28 x 1
    model.add(layers.Conv2D(16, kernel_size=(3, 3),
                            activation='relu', padding='same',
                            input_shape=input_shape))
    # 28 x 28 x 16
    model.add(layers.MaxPooling2D(pool_size=(2, 2), padding='same'))
    # 14 x 14 x 16
    model.add(layers.Conv2D(8, kernel_size=(3, 3),
                            activation='relu', padding='same'))
    # 14 x 14 x 8
    model.add(layers.MaxPooling2D(pool_size=(2, 2), padding='same'))
    # 7 x 7 x 8
    model.add(layers.Conv2D(8, kernel_size=(3, 3),
                            activation='relu', padding='same'))
    # 7 x 7 x 8
    model.add(layers.UpSampling2D(size=(2, 2)))
    # 14 x 14 x 8
    model.add(layers.Conv2D(16, kernel_size=(3, 3),
                            activation='relu', padding='same'))
    # 14 x 14 x 16
    model.add(layers.UpSampling2D(size=(2, 2)))
    # 28 x 28 x 16
    model.add(layers.Conv2D(1, kernel_size=(3, 3),
                            activation='sigmoid', padding='same'))
    # 28 x 28 x 1
    model.compile(optimizer='adam',
                  loss='binary_crossentropy')

    # モデルの構造を確認する
    print(model.summary())

    fit_callbacks = [
        callbacks.EarlyStopping(monitor='val_loss',
                                patience=5,
                                mode='min')
    ]

    # モデルを学習させる
    model.fit(x_train, x_train,
              epochs=1000,
              batch_size=4096,
              shuffle=True,
              validation_data=(x_test, x_test),
              callbacks=fit_callbacks,
              )

    # テストデータの損失を確認しておく
    score = model.evaluate(x_test, x_test, verbose=0)
    print('test xentropy:', score)

    # 学習済みのモデルを元に、次元圧縮だけするモデルを用意する
    encoder = models.clone_model(model)
    encoder.compile(optimizer='adam',
                    loss='binary_crossentropy')
    encoder.set_weights(model.get_weights())

    # 中間層までのレイヤーを取り除く
    encoder.pop()
    encoder.pop()
    encoder.pop()
    encoder.pop()

    # テストデータからランダムに 10 点を選び出す
    p = np.random.randint(0, len(x_test), 10)
    x_test_sampled = x_test[p]
    # 選びだしたサンプルを AutoEncoder にかける
    x_test_sampled_pred = model.predict_proba(x_test_sampled,
                                              verbose=0)
    # 次元圧縮だけする場合
    x_test_sampled_enc = encoder.predict_proba(x_test_sampled,
                                               verbose=0)

    # 処理結果を可視化する
    fig, axes = plt.subplots(3, 10, figsize=(12, 12))
    for i, label in enumerate(y_test[p]):
        # 元画像を上段に表示する
        img = x_test_sampled[i].reshape(image_height, image_width)
        axes[0][i].imshow(img, cmap=cm.gray_r)
        axes[0][i].axis('off')
        axes[0][i].set_title(label, color='red')
        # AutoEncoder で次元圧縮した画像を中段に表示する
        enc_img = x_test_sampled_enc[i].reshape(7, 7 * 8).T
        axes[1][i].imshow(enc_img, cmap=cm.gray_r)
        axes[1][i].axis('off')
        # AutoEncoder で復元した画像を下段に表示する
        pred_img = x_test_sampled_pred[i].reshape(image_height, image_width)
        axes[2][i].imshow(pred_img, cmap=cm.gray_r)
        axes[2][i].axis('off')

    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python cae.py
Model: "sequential_1"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv2d_11 (Conv2D)           (None, 28, 28, 16)        160       
_________________________________________________________________
max_pooling2d_5 (MaxPooling2 (None, 14, 14, 16)        0         
_________________________________________________________________
conv2d_12 (Conv2D)           (None, 14, 14, 8)         1160      
_________________________________________________________________
max_pooling2d_6 (MaxPooling2 (None, 7, 7, 8)           0         
_________________________________________________________________
conv2d_13 (Conv2D)           (None, 7, 7, 8)           584       
_________________________________________________________________
up_sampling2d_5 (UpSampling2 (None, 14, 14, 8)         0         
_________________________________________________________________
conv2d_14 (Conv2D)           (None, 14, 14, 16)        1168      
_________________________________________________________________
up_sampling2d_6 (UpSampling2 (None, 28, 28, 16)        0         
_________________________________________________________________
conv2d_15 (Conv2D)           (None, 28, 28, 1)         145       
=================================================================
Total params: 3,217
Trainable params: 3,217
Non-trainable params: 0
_________________________________________________________________
None
Train on 60000 samples, validate on 10000 samples
Epoch 1/1000
60000/60000 [==============================] - 1s 18us/step - loss: 0.6363 - val_loss: 0.5522
Epoch 2/1000
60000/60000 [==============================] - 1s 14us/step - loss: 0.5042 - val_loss: 0.4346
Epoch 3/1000
60000/60000 [==============================] - 1s 14us/step - loss: 0.3655 - val_loss: 0.2909

...(省略)...

Epoch 373/1000
60000/60000 [==============================] - 1s 14us/step - loss: 0.0728 - val_loss: 0.0721
Epoch 374/1000
60000/60000 [==============================] - 1s 14us/step - loss: 0.0727 - val_loss: 0.0721
Epoch 375/1000
60000/60000 [==============================] - 1s 14us/step - loss: 0.0727 - val_loss: 0.0721
test xentropy: 0.07211270518302917

検証用データに対するクロスエントロピーは約 0.072 と、前述した MLP の AutoEncoder よりも小さくなっていることがわかる。

以下は、上段が検証用データの画像、中段が AutoEncoder が次元圧縮した特徴の画像表現、下段が復元した画像となっている。

f:id:momijiame:20200416183423p:plain
Convolutional AutoEncoder の入出力の画像表現

入力に比べれば少しかすれているものの、ちゃんと復元できている。

Functional API を使った場合

おまけとして Functional API を使った例も以下に示す。 学習済みモデルから中間層の出力を取り出すところだけは Functional API を使う方法が分からなかったので Sequential API を使った。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
from keras import layers
from keras import models
from keras import callbacks
from keras import backend as K
from keras.datasets import mnist
from matplotlib import pyplot as plt
from matplotlib import cm


def main():
    (x_train, train), (x_test, y_test) = mnist.load_data()
    image_height, image_width = 28, 28

    if K.image_data_format() == 'channels_last':
        x_train = x_train.reshape(x_train.shape[0],
                                  image_height, image_width, 1)
        x_test = x_test.reshape(x_test.shape[0],
                                image_height, image_width, 1)
        input_shape = (image_height, image_width, 1)
    else:
        x_train = x_train.reshape(x_train.shape[0],
                                  1, image_height, image_width)
        x_test = x_test.reshape(x_test.shape[0],
                                1, image_height, image_width)
        input_shape = (1, image_height, image_width)

    x_train = x_train.astype('float32')
    x_test = x_test.astype('float32')
    x_train = (x_train - x_train.min()) / (x_train.max() - x_train.min())
    x_test = (x_test - x_test.min()) / (x_test.max() - x_test.min())

    # Functional API を使う場合
    input_ = layers.Input(shape=input_shape)
    # 28 x 28 x 1
    x = layers.Conv2D(16, kernel_size=(3, 3),
                      activation='relu', padding='same')(input_)
    # 28 x 28 x 16
    x = layers.MaxPooling2D(pool_size=(2, 2), padding='same')(x)
    # 14 x 14 x 16
    x = layers.Conv2D(8, kernel_size=(3, 3),
                      activation='relu', padding='same')(x)
    # 14 x 14 x 8
    x = layers.MaxPooling2D(pool_size=(2, 2), padding='same')(x)
    # 7 x 7 x 8
    x = layers.Conv2D(8, kernel_size=(3, 3),
                      activation='relu', padding='same')(x)
    # 7 x 7 x 8
    x = layers.UpSampling2D(size=(2, 2))(x)
    # 14 x 14 x 8
    x = layers.Conv2D(16, kernel_size=(3, 3),
                      activation='relu', padding='same')(x)
    # 14 x 14 x 16
    x = layers.UpSampling2D(size=(2, 2))(x)
    # 28 x 28 x 16
    output_ = layers.Conv2D(1, kernel_size=(3, 3),
                            activation='sigmoid', padding='same')(x)
    # 28 x 28 x 1
    model = models.Model(inputs=input_, output=output_)
    model.compile(optimizer='adam',
                  loss='binary_crossentropy')

    print(model.summary())

    fit_callbacks = [
        callbacks.EarlyStopping(monitor='val_loss',
                                patience=5,
                                mode='min')
    ]

    model.fit(x_train, x_train,
              epochs=1000,
              batch_size=4096,
              shuffle=True,
              validation_data=(x_test, x_test),
              callbacks=fit_callbacks,
              )

    score = model.evaluate(x_test, x_test, verbose=0)
    print('test xentropy:', score)

    encoder = models.clone_model(model)
    encoder.compile(optimizer='adam',
                    loss='binary_crossentropy')
    encoder.set_weights(model.get_weights())

    # Sequential API を使ってモデルを構築し直す
    encoder = models.Sequential()
    # 真ん中の層までを取り出す
    for layer in model.layers[:-4]:
        encoder.add(layer)

    p = np.random.randint(0, len(x_test), 10)
    x_test_sampled = x_test[p]
    x_test_sampled_pred = model.predict(x_test_sampled,
                                        verbose=0)
    x_test_sampled_enc = encoder.predict_proba(x_test_sampled,
                                               verbose=0)

    fig, axes = plt.subplots(3, 10, figsize=(12, 12))
    for i, label in enumerate(y_test[p]):
        img = x_test_sampled[i].reshape(image_height, image_width)
        axes[0][i].imshow(img, cmap=cm.gray_r)
        axes[0][i].axis('off')
        axes[0][i].set_title(label, color='red')

        enc_img = x_test_sampled_enc[i].reshape(7, 7 * 8).T
        axes[1][i].imshow(enc_img, cmap=cm.gray_r)
        axes[1][i].axis('off')

        pred_img = x_test_sampled_pred[i].reshape(image_height, image_width)
        axes[2][i].imshow(pred_img, cmap=cm.gray_r)
        axes[2][i].axis('off')

    plt.show()


if __name__ == '__main__':
    main()

結果は同じなので省略する。

いじょう。

Python: gensim で学習済み単語ベクトル表現を扱ってみる

Python で自然言語処理を扱うためのパッケージのひとつに gensim がある。 今回は、gensim で学習済み単語ベクトル表現 (Pre-trained Word Vectors) を使った Word Embedding を試してみた。 Word Embedding というのは単語 (Word) をベクトル表現の特徴量にする手法のこと。 ごく単純な One-Hot Encoding から、ニューラルネットワークの学習時の重みを元にしたものまで様々な手法が知られている。

今回は、その中でも Facebook の公開している fastText と呼ばれる学習済み単語ベクトル表現を使うことにした。 学習に使われているコーパスは Web クローラのデータと Wikipedia があるようだけど、とりあえずクローラの方にしてみる。 詳細については以下を参照のこと。

fasttext.cc

使った環境は次のとおり。

$ sw_vers           
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G3020
$ python -V             
Python 3.7.7

下準備

下準備として、まずは公開されている学習済み単語ベクトル表現をダウンロードしておく。 次の URL からは Web クローラで収集した日本語コーパスを使って、単語を 300 次元のベクトルに埋め込むためのモデルが得られる。 このファイルは 1GB ほどあるので割と時間がかかる。

$ wget https://dl.fbaipublicfiles.com/fasttext/vectors-crawl/cc.ja.300.vec.gz

なお、fastText の Python パッケージを使って、次のようにしても良い。 この機能は今のところリポジトリの HEAD にしか入っていないけど、次のリリースでは含まれるはず。

$ pip install -U git+https://github.com/facebookresearch/fastText.git
$ python -c "from fasttext import util; util.download_model('ja')"

そして、gensim パッケージをインストールしておく。

$ pip install gensim

インストールできたら、Python の REPL を起動しておく。

$ python

gensim から fastText の学習済み単語ベクトル表現を使う

はじめに gensim をインポートする。

>>> import gensim

先ほどダウンロードしたファイルを gensim で読み込む。 これで単語をベクトル表現に変換するためのモデルができる。 なお、この操作にはかなり時間がかかるので気長に待つ必要がある。 だいたい 10 分くらい 1

>>> model = gensim.models.KeyedVectors.load_word2vec_format('cc.ja.300.vec.gz', binary=False)
>>> model
<gensim.models.keyedvectors.Word2VecKeyedVectors object at 0x127797310>

このモデルのざっくりした使い方は以下のドキュメントに記載されている。 ただ、ソースコードと見比べていくとイマイチ更新が追いついていない雰囲気も伺える。

radimrehurek.com

モデルの読み込みが終わったら、早速単語をベクトル表現にしてみよう。 試しに「猫」という単語をベクトル表現にすると、次のようになる。 モデルには辞書ライクなインターフェースがあるので、ブラケットを使ってアクセスする。

>>> model['猫']
array([-2.618e-01, -7.520e-02, -1.930e-02,  2.088e-01, -3.005e-01,
        1.936e-01, -1.561e-01, -3.540e-02,  1.220e-01,  2.718e-01,
        7.460e-02,  1.356e-01,  2.299e-01,  1.851e-01, -2.684e-01,
...(snip)...
       -9.260e-02, -8.890e-02,  1.143e-01, -3.381e-01, -1.913e-01,
        8.160e-02, -4.420e-02, -2.405e-01, -2.170e-02, -1.062e-01,
        3.230e-02, -2.380e-02,  1.860e-02, -2.750e-02, -1.900e-01],
      dtype=float32)

確認すると、この単語ベクトルはちゃんと 300 次元ある。

>>> model['猫'].shape
(300,)

学習した単語のボキャブラリーには vocab という名前でアクセスできる。 このモデルでは 200 万の単語を学習しているようだ。

>>> list(model.vocab.keys())[:10]
['の', '、', '。', 'に', 'は', 'が', 'を', 'て', 'た', 'で']
>>> len(model.vocab.keys())
2000000

most_similar() メソッドを使うと、特定の単語に似ている単語が得られる。

>>> from pprint import pprint
>>> pprint(model.most_similar('猫', topn=10))
[('ネコ', 0.8059155941009521),
 ('ねこ', 0.7272598147392273),
 ('子猫', 0.720253586769104),
 ('仔猫', 0.7062687873840332),
 ('ニャンコ', 0.7058036923408508),
 ('野良猫', 0.7030349969863892),
 ('犬', 0.6505385041236877),
 ('ミケ', 0.6356303691864014),
 ('野良ねこ', 0.6340526342391968),
 ('飼猫', 0.6265145540237427)]

単語の類似度は similarity() メソッドで得られる。 たとえば、「猫」と「犬」の方が「猫」と「人」よりも単語としては似ているようだ。

>>> model.similarity('猫', '犬')
0.65053856
>>> model.similarity('猫', '人')
0.23371725

ところで、Word2vec なんかだと、まるで単語の意味を理解しているかのように単語ベクトル間の足し引きができるとかよく言われている。 このモデルではどんなもんだろうか。

ありがちな例として「王様」から「男」を抜いて「女」を足してみよう。 女王になるみたいな説明、よく見る。

>>> new_vec = model['王様'] - model['男'] + model['女']

similar_by_vector() メソッドを使って得られたベクトルに似ている単語を確認してみる。

>>> pprint(model.similar_by_vector(new_vec))
[('王様', 0.8916897773742676),
 ('女王', 0.527921199798584),
 ('ラジオキッズ', 0.5255386829376221),
 ('王さま', 0.5226017236709595),
 ('王妃', 0.5000214576721191),
 ('裸', 0.487439900636673),
 ('タプチム', 0.4832267761230469),
 ('アンナ・レオノーウェンズ', 0.4807651937007904),
 ('ゲムケン', 0.48058977723121643),
 ('お姫様', 0.4792743921279907)]

最も似ている単語は変換前の単語ということになってしまった。 変換前の単語は取り除くことにするのかな? だとしたら、次に「女王」が得られているけど。

そもそも元の単語ベクトルに類似した単語は次のとおり。

>>> pprint(model.similar_by_word('王様'))
[('王さま', 0.5732064247131348),
 ('ラジオキッズ', 0.5712020993232727),
 ('女王', 0.5381238460540771),
 ('ゲムケン', 0.5354433059692383),
 ('裸', 0.5288593769073486),
 ('ブランチ', 0.5283758640289307),
 ('王', 0.5277552604675293),
 ('タプチム', 0.5203500986099243),
 ('乃女', 0.5120265483856201),
 ('王妃', 0.5055921077728271)]

学習済み単語ベクトル表現の読み込みが遅い問題について

一般的に、公開されている学習済み単語ベクトル表現は、特定の言語や環境に依存しないフォーマットになっているようだ。 そのため、gensim に限らず読み込む際には結構重たいパース処理が必要になっている。 この点は、一般的な機械学習において CSV などのデータを読み込むときのテクニックを応用すれば高速化が見込めそうだ。

具体的には Python であれば Pickle を使って直列化・非直列化する。 まずは pickle パッケージを読み込む。

>>> import pickle

ファイル名を指定してモデルをストレージに直列化する。

>>> with open('gensim-kvecs.cc.ja.300.vec.pkl', mode='wb') as fp:
...     pickle.dump(model, fp)
... 

直列化したら、あとは使いたいタイミングで非直列化するだけ。 ストレージのスループット次第だけど、それでも 10 分待たされるようなことはないはず。

>>> with open('gensim-kvecs.cc.ja.300.vec.pkl', mode='rb') as fp:
...     model = pickle.load(fp)
... 

ただし、直列化されない内部的なキャッシュもあるようで一発目の呼び出しではちょっとモタつく。 とはいえ、もちろん 10 分とか待たされるようなことはない。

いじょう。


  1. 今の実装には明確な処理のボトルネックがあって、それを解消すると 1/3 くらいの時間に短縮できる。

Python: statsmodels で時系列データを基本成分に分解する

時系列データを扱うとき、原系列が傾向変動・季節変動・不規則変動という基本成分の合成で成り立っていると捉えることがある。 傾向変動は中長期的な増加・減少といった変化であり、季節変動は例えば 1 ヶ月や 1 年といった周期的な変化を指している。 不規則変動は、前者 2 つに当てはまらない変化で、誤差変動と特異的変動に分けて考える場合もあるようだ。

原系列が基本成分の合成と考える場合でも、捉え方として加法モデルと乗法モデルにさらに分かれる。 まず、加法モデルでは傾向変動  T(t) と季節変動  S(t)、誤差変動  I(t) の和を原系列  O(t) と捉える。 この考え方では、傾向変動の大きさに関係なく一定の季節変動が加えられる。

 O(t) = T(t) + S(t) + I(t)

一方で、乗法モデルでは積と捉える。 つまり、傾向変動が大きくなれば、それに比例して季節変動も大きくなる、という考え方。

 O(t) = T(t) \times S(t) \times I(t)

これらは扱うデータやタスクによって使い分ける必要があるらしい。 今回は、Python の statsmodels を使って原系列を基本成分に分解してみる。

使った環境は次のとおり。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G3020
$ python -V                      
Python 3.7.7

下準備

まずは、必要なパッケージをインストールしておく。

$ pip install statsmodels seaborn

現系列を基本成分に分解する

現系列を基本成分に分解するには seasonal_decompose() 関数を使う。 次のサンプルコードでは、旅客機の乗客数のデータを分解してみた。 なお、モデルとしては乗法モデルを仮定している。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import pandas as pd
import seaborn as sns
from statsmodels import api as sm
from matplotlib import pyplot as plt
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()


def main():
    # 旅客機の乗客数のデータセットを読み込む
    df = sns.load_dataset('flights')

    df['year-month'] = df.month.astype(str) + ', ' + df.year.astype(str)
    df['year-month'] = pd.to_datetime(df['year-month'], format='%B, %Y')
    df = df.set_index('year-month')

    # 時系列データを傾向変動・季節変動・残差に分解する
    decompose_result = sm.tsa.seasonal_decompose(df['passengers'],
                                                 # 乗法モデルを仮定する
                                                 model='multiplicative')

    # これでもグラフが描ける
    # decompose_result.plot()

    # 描画する領域を用意する
    fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(8, 8), sharex=True)

    # 原系列
    axes[0].set_title('Observed')
    axes[0].plot(decompose_result.observed)

    # 傾向変動
    axes[1].set_title('Trend')
    axes[1].plot(decompose_result.trend)

    # 季節変動
    axes[2].set_title('Seasonal')
    axes[2].plot(decompose_result.seasonal)

    # 残差 (不規則変動 = 誤差変動 + 特異的変動)
    axes[3].set_title('Residual')
    axes[3].plot(decompose_result.resid)

    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python decompose.py

次のようなグラフが得られる。 上から原系列、傾向変動、季節変動、不規則変動 (残差) となっている。

f:id:momijiame:20200406232249p:plain
原系列を基本成分に分解する

注意すべきポイントとして、seasonal_decompose() 関数では傾向変動と不規則変動は端の要素が NaN になってしまうようだ。

原系列と傾向変動で偏自己相関を見比べてみる

時系列データに周期性、つまり季節成分が含まれるか調べる方法のひとつに、偏自己相関を調べるというのがある。 今回は、分解する前の原系列と、分解した後の傾向変動について、自己相関と偏自己相関を調べて見比べてみよう。 自己相関では、ある時点のデータ  O(t) には 1 つ前のデータ  O(t - 1) が関係している可能性がある。 そして、1 つ前のデータには、さらにその 1 つ前のデータが、というような積み重ねの影響を受けている恐れがある。 偏自己相関は、そのような積み重ねの影響を消去することを目的とした統計量になっている。

次のサンプルコードでは、分解する前の原系列と分解した後の傾向変動について自己相関と偏自己相関を計算してグラフにプロットしている。 なお、データに NaN があると自己相関がうまく計算できないようなので、傾向変動からは NaN を除外している。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import pandas as pd
import seaborn as sns
from statsmodels import api as sm
from matplotlib import pyplot as plt
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()


def main():
    # 旅客機の乗客数のデータセットを読み込む
    df = sns.load_dataset('flights')

    df['year-month'] = df.month.astype(str) + ', ' + df.year.astype(str)
    df['year-month'] = pd.to_datetime(df['year-month'], format='%B, %Y')
    df = df.set_index('year-month')

    # 時系列データを傾向変動・季節変動・残差に分解する
    decompose_result = sm.tsa.seasonal_decompose(df['passengers'],
                                                 # 乗法モデルを仮定する
                                                 model='multiplicative')

    _, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 8))

    # 原系列の ACF
    sm.tsa.graphics.plot_acf(df['passengers'], ax=axes[0][0])
    # 原系列の PACF
    sm.tsa.graphics.plot_pacf(df['passengers'], ax=axes[1][0])

    # 傾向変動の ACF
    sm.tsa.graphics.plot_acf(decompose_result.trend.dropna(), ax=axes[0][1])
    # 傾向変動の PACF
    sm.tsa.graphics.plot_pacf(decompose_result.trend.dropna(), ax=axes[1][1])

    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python decompacf.py

次のようなグラフが得られる。 左側が原系列、右側が傾向変動について計算したもの。 青い帯は「相関がない」を帰無仮説とした 95% 信頼区間を表している。 つまり、帯の外にある値は 5% の有意水準で相関があることになる。

f:id:momijiame:20200406233714p:plain
原系列と傾向変動の自己相関 (ACF) と偏自己相関 (PACF)

原系列の偏自己相関では、ところどころに相関があることがわかる。 それに対して、傾向変動の偏自己相関では、ラグ 1 以外に相関がある箇所は見当たらない。 つまり、季節変動を取り除くことができていることになるようだ。

いじょう。

Python: 中心化移動平均 (CMA: Centered Moving Average) について

以前から移動平均 (MA: Moving Average) という手法自体は知っていたけど、中心化移動平均 (CMA: Centered Moving Average) というものがあることは知らなかった。 一般的な移動平均である後方移動平均は、データの対応関係が原系列に対して遅れてしまう。 そこで、中心化移動平均という手法を使うことで遅れをなくすらしい。 この手法は、たとえば次のような用途でひとつのやり方として使われているようだ。

  • 不規則変動の除去
  • 季節変動の除去

使った環境は次のとおり。

$ sw_vers  
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G3020
$ python -V
Python 3.7.7

下準備

下準備として、必要なパッケージをインストールしておく。

$ pip install pandas seaborn

インストールできたら Python のインタプリタを起動する。

$ python

起動できたら、1 ~ 12 までの数字が入った Series のオブジェクトを作る。

>>> import pandas as pd
>>> x = pd.Series(range(1, 12 + 1))
>>> x
0      1
1      2
2      3
3      4
4      5
5      6
6      7
7      8
8      9
9     10
10    11
11    12
dtype: int64

上記のオブジェクトを使って移動平均の計算方法について学んでいく。

(後方) 移動平均 (MA: Moving Average)

はじめに、一般的な移動平均である後方移動平均から試す。 Pandas では rolling() メソッドを使うことで移動平均を計算できる。 以下では 3 点の要素で平均をとる 3 点移動平均を計算している。

>>> backward_3p_ma = x.rolling(window=3).mean()
>>> backward_3p_ma
0      NaN
1      NaN
2      2.0
3      3.0
4      4.0
5      5.0
6      6.0
7      7.0
8      8.0
9      9.0
10    10.0
11    11.0
dtype: float64

この計算では、たとえばインデックスで 2 に入る要素は、次のように計算される。

backward_3p_ma[2] = (x[0] + x[1] + x[2]) / 3

つまり、自分よりも後ろの値を使って平均を計算することになる。

この計算方法では、原系列との対応関係を考えたとき次のように 1 つ分の遅れが出る。

>>> pd.concat([x, backward_3p_ma], axis=1)
     0     1
0    1   NaN
1    2   NaN
2    3   2.0
3    4   3.0
4    5   4.0
5    6   5.0
6    7   6.0
7    8   7.0
8    9   8.0
9   10   9.0
10  11  10.0
11  12  11.0

中心化移動平均 (CMA: Centered Moving Average)

そこで、遅れの影響を取り除いたものが中心化移動平均と呼ばれるらしい。 たとえば、平均をとる要素数が奇数のときは、単に rolling() メソッドの center オプションを有効にするだけで良い。

>>> centered_3p_ma = x.rolling(window=3, center=True).mean()
>>> centered_3p_ma
0      NaN
1      2.0
2      3.0
3      4.0
4      5.0
5      6.0
6      7.0
7      8.0
8      9.0
9     10.0
10    11.0
11     NaN
dtype: float64

このオプションが有効だと、たとえばインデックスで 2 に入る要素は、次のように計算される。 つまり、自分の前後にある値を使って平均を計算することになる。

centered_3p_ma[2] = (x[1] + x[2] + x[3]) / 3

現系列との対応関係を確認すると、このやり方では遅れがなくなっている。

>>> pd.concat([x, centered_3p_ma], axis=1)
     0     1
0    1   NaN
1    2   2.0
2    3   3.0
3    4   4.0
4    5   5.0
5    6   6.0
6    7   7.0
7    8   8.0
8    9   9.0
9   10  10.0
10  11  11.0
11  12   NaN

ただし、このオプションは実装的には単に後方移動平均を計算した上で shift() しているだけに過ぎないらしい。 ようするに、こういうこと。

>>> backward_3p_ma_shifted = x.rolling(window=3).mean().shift(-1)
>>> pd.concat([x, backward_3p_ma_shifted], axis=1)
     0     1
0    1   NaN
1    2   2.0
2    3   3.0
3    4   4.0
4    5   5.0
5    6   6.0
6    7   7.0
7    8   8.0
8    9   9.0
9   10  10.0
10  11  11.0
11  12   NaN

そのため、平均をとる要素数が偶数のときに困ったことが起こる。

平均をとる要素数が偶数のときの問題点について

試しに、平均をとる要素数を 4 点に増やして、そのまま計算してみよう。

>>> centerd_4p_ma = x.rolling(window=4, center=True).mean() 
>>> centerd_4p_ma
0      NaN
1      NaN
2      2.5
3      3.5
4      4.5
5      5.5
6      6.5
7      7.5
8      8.5
9      9.5
10    10.5
11     NaN
dtype: float64

すると、計算結果に端数が出ている。

原系列と比較すると、対応関係に 0.5 の遅れがあることがわかる。

>>> pd.concat([x, centerd_4p_ma], axis=1)
     0     1
0    1   NaN
1    2   NaN
2    3   2.5
3    4   3.5
4    5   4.5
5    6   5.5
6    7   6.5
7    8   7.5
8    9   8.5
9   10   9.5
10  11  10.5
11  12   NaN

つまり、中心化移動平均では平均をとる要素数が偶数と奇数のときで計算方法を変えなければいけない。

要素数が偶数のときの計算方法

要素数が偶数のときの中心化移動平均は、計算が 2 段階に分かれている。 はじめに、後方移動平均をそのまま計算する。

>>> backward_4p_ma =x.rolling(window=4).mean()
>>> backward_4p_ma
0      NaN
1      NaN
2      NaN
3      2.5
4      3.5
5      4.5
6      5.5
7      6.5
8      7.5
9      8.5
10     9.5
11    10.5
dtype: float64

その上で、移動平均に使った要素数の半分だけデータをずらし、もう一度要素数 2 で平均を取り直す。

centerd_4p_ma = backward_4p_ma.shift(-2).rolling(window=2).mean()

原系列との対応関係を確認すると、遅れが解消していることがわかる。

>>> pd.concat([x, centerd_4p_ma], axis=1)
     0     1
0    1   NaN
1    2   NaN
2    3   3.0
3    4   4.0
4    5   5.0
5    6   6.0
6    7   7.0
7    8   8.0
8    9   9.0
9   10  10.0
10  11   NaN
11  12   NaN

別のデータで中心化移動平均を計算してみる

もうちょっとちゃんとしたデータでも計算してみよう。 次のサンプルコードでは、旅客機の乗客数の推移に対して 12 点で中心化移動平均を計算している。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from calendar import month_name

import seaborn as sns
from matplotlib import pyplot as plt
import pandas as pd
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()


def main():
    # 航空機の乗客を記録したデータセットを読み込む
    df = sns.load_dataset('flights')
    # 月の名前を数字に直す
    month_name_mappings = {name: str(n).zfill(2) for n, name in
                           enumerate(month_name)}
    df['month'] = df['month'].apply(lambda x: month_name_mappings[x])
    # 年と月を結合したカラムを用意する
    df['year-month'] = df.year.astype(str) + '-' + df.month.astype(str)
    df['year-month'] = pd.to_datetime(df['year-month'], format='%Y-%m')

    # 原系列
    sns.lineplot(data=df, x='year-month', y='passengers', label='original')
    # 中心化移動平均
    df['passengers-cma'] = df['passengers'].rolling(window=12).mean().shift(-6).rolling(2).mean()
    sns.lineplot(data=df, x='year-month', y='passengers-cma', label='CMA')

    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみる。

$ python flightscma.py

得られるグラフが次のとおり。

f:id:momijiame:20200401184339p:plain
旅客機の乗客数データに対して 12 点中心化移動平均を計算する

ここまでデータ点数が多いと、正直 0.5 の遅れとか言われてもよくわからないけど、これで計算できているはず。