分析ノート

2019-04-152019-04-13

記事数が100を超えていました

ついこの間50記事を超えたことを書いたばかりなのですが、
その後も1日1記事のペースを維持でき、めでたく100記事を突破していました。
参考:ブログ記事数が50を超えていたので振り返り

2ヶ月弱しか経っていないのですが、それでも少々変わってきたことがあるのでその辺中心に振り返っておきます。

アクセス数が安定してきた

50記事の頃は毎日2~3人しか訪問者がいらっしゃいませんでしたが、
最近は安定して20〜30人くらいのかたが訪問してくださっているようです。

初回訪問のかたの多くはMeCabのインストール時のエラーや
grahpvizのイストール方法を探しているようです。

2019年第1四半期によく読まれた記事
の１位２位ですね。

一方で時系列解析やpandasなどのプログラミングの話、kerasなどの機械学習記事の人気はイマイチ。
まだチュートリアルレベルのことしか書けてない僕のコンテンツ力不足ですね。

$TeX$の入力に慣れてきた

大学/大学院時代に使って以来すっかり触らなくなっていた$TeX$ですが、
ARMA過程関連の記事を書く中で徐々に勘を取り戻してきました。

最初は複数行の数式の$=$の位置を揃える方法も行列の書き方も、
ギリシャ文字の入力も不等号もアクセント記号もすっかり忘れていたのですが、
徐々に楽に使えるようになってきていい感じです。

正直、最初は数式が多くなりそうってだけで書くのを敬遠したテーマもあったので、
改めてそれらの記事化も検討します。

やっぱり勉強になる(再掲)

50記事の時も書きましたが、新たに知った内容だけでなく、元々理解しているつもりの内容であっても、
記事にするための再検証の中で新たな学びが多くありました。

ARMA過程の一連の記事など、もともとはstatsumodelsの使い方だけ紹介して終わろうと思ってたところ、
最低限の用語だけでも説明せねばと思って、結構前に読んだ沖本本を引っ張り出してきて改めて書いたものです。

沖本本をさっと通読したときは「フンフン確かにこうなるね」と流していたところを、
一個一個の数式を手計算で追いかけ、実例を構築してプログラムを書いて検証して、
とやっていく中で、思ったほどうまくいかないことにも結構遭遇しつつ理解を深められて良かったです。
単位根過程やベクトル自己回帰の話題はまだこのブログに出せていないので追い追い書いていきます。

今後の計画

手元のメモから転記するだけで記事になるネタは減ってきたのですが、
それでもコードや数式が入ったブログを書くこと自体には慣れてきて負荷が下がってきたので、
もうしばらく今のペースで記事を書きたいと思ってます。

また、データサイエンティストとして働きはじめて2年経ちましたので、
仕事やキャリアの話なども、これから目指す人の参考になるように書いていきたいですね。

2019-04-142021-01-05

pythonでMD5

業務でMD5を使って文字列をハッシュ化する機会があったのでそのメモです。
(最近ではMD5自体があまり安全では無く、暗号学的な強度が必要なときは使うべきでは無いのでご注意ください。)

python3では、標準ライブラリにhashlibというものがあり、これにMD5が実装されています。
ドキュメント: hashlib — セキュアハッシュおよびメッセージダイジェスト

使い方は簡単で、インポートして hashlib.md5(“text”) するだけです、
となれば楽なのですが、これでは動きません。こんなふうにエラーになります。


import hashlib
hashlib.md5("text")

# 以下出力
TypeError                                 Traceback (most recent call last)
 in ()
----> 1 hashlib.md5("text")

TypeError: Unicode-objects must be encoded before hashing

これは、bytes型でデータを渡す必要があるからです。
そのため、b"text"か"text".encode("utf-8")のようにする必要があります。

これで一応エラーは無くなりますが、まだ少し不十分です。
というのも、戻り値がHASH object という型で返ってきます。


import hashlib
text = "Hello, World!"
bytes_text = text.encode("utf-8")
print(hashlib.md5(bytes_text))

# 出力
<md5 HASH object @ 0x10e3f7620>

欲しいのは、 16進法で表記された文字列だと思うので、もう一手間かけて、hexdigestという関数を実行します。


import hashlib
text = "Hello, World!"
bytes_text = text.encode("utf-8")
hash_obj = hashlib.md5(bytes_text)
print(hash_obj.hexdigest())

# 出力
65a8e27d8879283831b664bd8b7f0ad4

# 1行でやる方法
print(hashlib.md5(b"Hello, World!").hexdigest())

# 出力は同じ
65a8e27d8879283831b664bd8b7f0ad4

これでできました。

sha1やsha256もhashlibに実装されているので同じように使えます。

2019-04-132019-04-11

pythonでARMAモデルの推定

以前、ARモデルの推定をstatsmodelsで行う方法を書きました。
pythonでARモデルの推定

今回行うのははARMAモデルの推定です。
結果に表示されるconstが、予想してた値にならず結構苦戦しました。
これ、ARモデルと仕様が違って、定数項ではなく過程の期待値が返されるんですね。

密かにライブラリのバグじゃ無いのかと思ってます。ちなみにバージョンはこれ。


$ pip freeze | grep statsmodels
statsmodels==0.9.0

ドキュメントはこちら
statsmodels.tsa.arima_model.ARMA

ARモデルとは使い方が色々と異なります。
例えば、ARモデルでは、次数を推定する関数(select_order)をモデルが持っていましたが、
ARMAモデルにはなく、
arma_order_select_ic
という別のところ(stattools)に準備された関数を使います。

今回もARMA過程を一つ固定し、データを準備して行いますが、過程の次数が大きくなるとなかなか上手くいかないので、
ARMA(1, 1)過程のサンプルを作りました。

$$
y_t = 50 + 0.8y_{t-1} + \varepsilon_{t} + 0.4\varepsilon_{t-1} \\
\varepsilon_{t} \sim W.N.(2^2)
$$
この過程の期待値は $50/(1-0.8)=250$です。

それではデータの生成から行います。


T = 200
phi_1 = 0.8
theta_1 = 0.4
c = 50
sigma = 2
# 過程の期待値
mu = c/(1-phi_1)

# 撹乱項生成
epsilons = np.random.normal(0, sigma, size=T)

arma_data = np.zeros(T)
arma_data[0] = mu + epsilons[0]

for t in range(1, T):
    arma_data[t] = c+phi_1*arma_data[t-1]+epsilons[t]+theta_1*epsilons[t-1]

# データを可視化し、定常過程である様子を見る (記事中では出力は略します)
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 1, 1)
ax.plot(y)
plt.show()

これで検証用データが作成できましたので、
早速ARMAモデルを構築して係数と撹乱項の分散がもとまるのかやってみましょう。
今回も次数の決定にはAIC基準を使いました。


# 次数の推定
print(sm.tsa.arma_order_select_ic(arma_data, max_ar=2, max_ma=2, ic='aic'))

# 結果
'''
{'aic':              0           1            2
0  1107.937296  920.311164  1341.024062
1   844.276284  818.040579   820.039660
2   821.624647  820.039106   821.993934, 'aic_min_order': (1, 1)}
'''

# ARMAモデルの作成と推定
arma_model = sm.tsa.ARMA(y, order=(1, 1))
result = arma_model.fit()

最後に結果を表示します。


# 結果の表示
print(result.summary())

# 出力
'''
                              ARMA Model Results                              
==============================================================================
Dep. Variable:                      y   No. Observations:                  200
Model:                     ARMA(1, 1)   Log Likelihood                -416.830
Method:                       css-mle   S.D. of innovations              1.937
Date:                Thu, 11 Apr 2019   AIC                            841.660
Time:                        00:18:59   BIC                            854.853
Sample:                             0   HQIC                           846.999
                                                                              
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const        251.5200      0.866    290.374      0.000     249.822     253.218
ar.L1.y        0.7763      0.049     15.961      0.000       0.681       0.872
ma.L1.y        0.4409      0.067      6.622      0.000       0.310       0.571
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            1.2882           +0.0000j            1.2882            0.0000
MA.1           -2.2680           +0.0000j            2.2680            0.5000
-----------------------------------------------------------------------------
'''

自己回帰モデル部分のラグ1の係数が0.7763 (正解は0.8)、
移動平均モデル部分のラグ1の係数は0.4409 （正解は0.4)なので、そこそこ正しいですね。
撹乱項の標準偏差(S.D. of innovations) も 1.937となり、正解の2に近い値が得られています。

constには50に近い値になるはずだと思っていたのですが、
ここは過程の式の定数項ではなく期待値が出力されるようです。

公式ドキュメント内で該当の説明を見つけることができていませんが、
色々なモデルで検証したのでおそらく間違い無いでしょう。

また、次の関数を使って、元のデータとモデルの予測値をまとめて可視化できます。
result.plot_predict()

今回のモデルは非常に予測しやすいものを例として選んだので、大体近しい値になっていますね。

2019-04-122019-04-10

pythonで累積和

数列に対して、最初の項からn番目の項まで足したものの数列が必要になることがあります。
日々の売上データからその日までの累計の売上データを出すような計算ですね。

イメージとしては、
1,2,3,4,5,…に対して、1,3,6,10,15,… を得るような操作です。

スライスと内包表記を組み合わせてもいいし、足し算くらいならfor文を回しても問題ないのですが、
numpyやpandas に専用の関数が用意されているのでその紹介です。

ドキュメントはこの辺
numpy.cumsum
pandas.Series.cumsum

早速実行してみました。


import numpy as np
import pandas as pd

ary = np.arange(1, 11)
print(ary)
# [ 1  2  3  4  5  6  7  8  9 10]

print(ary.cumsum())
# [ 1  3  6 10 15 21 28 36 45 55]

series = pd.Series(np.arange(1, 11))
print(series)
'''
0     1
1     2
2     3
3     4
4     5
5     6
6     7
7     8
8     9
9    10
dtype: int64
'''

print(series.cumsum())
'''
0     1
1     3
2     6
3    10
4    15
5    21
6    28
7    36
8    45
9    55
dtype: int64
'''

よく似た関数として、
numpyには累積の積を返すcumprodが用意されており、
さらにpandasの方にはcumprodに加えて、そこまでの最大値と最小値を得られるcummax, cumminが用意されています。

2019-04-112019-04-08

辞書型のオブジェクトをJSON文字列に変換する

仕事で必要になったのでメモです。

pythonには jsonという標準ライブラリがあり、
それを使うことで、配列や辞書型のオブジェクトをJSON文字列に簡単に変換することができます。

ドキュメントはこちら。
json — JSON エンコーダおよびデコーダ

インポートして dumps関数を使うだけなので早速やってみましょう。


import json
data = {
        "key1": "value1",
        "key2": "value2",
        "ary1": ["element1", "element2"],
    }
json_data = json.dumps(data)
print(json_data)

# 出力
# {"key1": "value1", "key2": "value2", "ary1": ["element1", "element2"]}

ちなみに、json.dumpsを使わず、str(data)で文字列に変換すると、結果は
"{'key1': 'value1', 'key2': 'value2', 'ary1': ['element1', 'element2']}"
になります。
JSONのルールでは文字列はダブルクオーテーションで囲まないといけないので、
これはJSONではありません。

厳密にJSON型を要求する関数やAPIに渡すときはこれだと受け付けられないので、
json.dumpsを使いましょう。

2019-04-092019-04-07

numpyの乱数生成関数の設計について

このブログに登場するコードでも頻繁にnumpyで乱数を生成していますが、
そのドキュメントは一回読んでおいたほうがいいよという話です。

Random sampling (numpy.random)

pythonの学習を始めた頃は、本に載っているサンプルコードや検索したらでてくる各サイトを参考に、
個々の関数の使い方を覚えて使ってました。

例えばrandであれば次のような感じ。


import numpy as np
# スカラーで一つ値が欲しい場合
np.random.rand()
# 配列で5個値が欲しい場合
np.random.rand(5)
# 引数を複数指定すると多次元配列で乱数が生成される
np.random.rand(2, 2, 2)

標準正規分布に従う乱数が欲しい時はrandnがあり、randと同じように使えます。

その一方で、normalという関数も準備されていて、
これは使い方が違います。

引数が
normal([loc, scale, size])となっていて、
最初に平均と標準偏差を与える必要があり、
欲しい乱数の個数は3個目の引数に与えます。

一方、randやrandnは複数個の乱数を得るには便利ですが、パラメータを渡すことができません。
(なので、得た乱数に対して、定数を足したり掛けたりして調整する)

特に難しくも複雑でもないので、本に出てきた通りに暗記して使っていましたが、
関数の設計が統一されてないので不便だとずっと思ってました。

それでこの度ドキュメントをみたのですが、各関数が、
Simple random data と、Distributionsというカテゴリに分けて定義されていることを知りました。
(あと二つ、PermutationsとRandom generatorがありますが割愛)
そして、(完全に揃ってる訳ではないのですが、)それぞれのカテゴリ内である程度使い方を揃えて作られていることがわかりました。

要は、Simple random dataのほうは欲しいデータの数だけ渡せばよく、
Distributionsの方は最初に分布のパラメーターを指定して、データの個数はsize引数で渡す。

今更なんですが、最初にpython勉強し始めた頃に、
公式ドキュメントもきちんと読んでおけば、もう少し楽に覚えられたかなと思いますね。

2019-04-082019-04-09

pythonでARモデルの推定

時系列解析で理論の話題が続いていましたが、用語の定義がだいぶ済んだので、
そろそろ実際にコードを書いてモデルの推定等をやってみます。
(とはいえAIC,BICなどのかなり重要なキーワードや、推定の理論面の話をまだ紹介してないのですが。)

今回はAR(2)過程のデータを生成し、statusmodelsを使って、そのパラメーターを推定します。
真のモデルが不明なサンプルデータではなくてもともと正解がわかってるデータで雰囲気をつかもうというのが主目的です。

今回使うのは、次の定常AR(2)過程です。
$$
y_t = 5 + 1.4y_{t-1} – 0.48y_{t-2} + \varepsilon_{t}, \ \ \varepsilon_{t} \sim W.N.(0.5^2)
$$

早速データを生成して可視化しておきましょう。(過去の記事の例と大して変わらないので、可視化の結果は略します。)


import statsmodels.api as sm
import matplotlib.pyplot as plt
import numpy as np
phi_1 = 1.4
phi_2 = -0.48
c = 5
sigma = 0.5
T = 500
# 過程の期待値
mu = c/(1-phi_1-phi_2)


# データの生成
ar_data = np.zeros(T)
ar_data[0] = mu + np.random.normal(0, sigma)
ar_data[1] = mu + np.random.normal(0, sigma)
for t in range(2, T):
    ar_data[t] = c + phi_1*ar_data[t-1] + phi_2*ar_data[t-2] \
                                + np.random.normal(0, sigma)
# データの可視化 (記事中では出力は略します)
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 1, 1)
ax.plot(ar_data)
plt.show()

あとはこのデータから、元のパラメータである、定数項の5や、自己回帰係数の1.4と-0.48、撹乱項の分散0.25(=0.5^2)などを
推定できれば成功です。

使うのはこれ。
statsmodels.tsa.ar_model.AR
fitの戻り値であるARResultsのドキュメントも読んでおきましょう。
statsmodels.tsa.ar_model.ARResults

早速やってみます。
なお、手順は以下の通りです。

ARモデルの生成
次数の決定(select_order, 今回はAICを使用)
決定した次数で推定
推定結果の確認


# モデルの生成
model = sm.tsa.AR(ar_data)
# AICでモデルの次数を選択
print(model.select_order(maxlag=6, ic='aic'))  # 出力:2

# 推定
result = model.fit(maxlag=2)

# モデルが推定したパラメーター
print(result.params)
# 出力
# [ 5.68403723  1.38825575 -0.47929503]
print(result.sigma2)
# 出力
# 0.23797695866882768

今回は正解にそこそこ近い結果を得ることができました。
乱数を使うので当然ですが実行するたびに結果は変わります。
また、T=500くらいのデータ量だと、モデルの次数が異なる値になることも多く、
結構上手くいかないことがあります。

実行するたびに全然結果が変わるのと、
実際の運用では同一条件のデータを500件も集められることは滅多にないので、
正直、ARモデルでの推定結果をあまり過信しすぎないように気をつけようと思いました。

2019-04-072019-04-06

LightsailのMySQLに接続する

このブログはAWSのLightsailを使って構築しています。
ちょっと訳あって、そのサーバーのDB(MySQL)に接続したくなったのですが、
少しつまずいたのでメモです。

mysqlは勝手に構築してくれているので初期パスワードに何が設定されているのか、わかりませんでした。

結論から言うと
user : root
pass : bitnami_application_password ファイルの中身
で接続できます。
このwordpressの管理画面と同じ初期パスワードだったんですね。

わかるまで、サーバーにログインし、

\$ mysql
\$ mysql -u user
\$ mysql -u root
など試しましたが、ログインできず。

パスワードとして心当たりがあるのが管理画面のパスワードだけだったので、

\$ mysql -u root -p
とコマンドを打った後、bitnami_application_passwordの中身を入力したら接続できました。

念のため設定変えておこう。

2019-04-062019-04-04

定常AR過程をMA過程で表現する

前回の記事がMA過程の反転可能性について述べたものだったので、今回はAR過程をMA過程で表現できることについてです。
参考:MA過程の反転可能性

AR過程全てではなく、定常AR過程に限った性質ですが、
定常AR過程はMA($\infty$)過程に書き直すことができます。

たとえば、 $|\phi|<1$の時、次のAR(1)過程は定常です。 $$ y_t = \phi y_{t-1} + \varepsilon_t , \varepsilon_t\sim W.N.(\sigma^2). $$ これは、 $$ y_t = \sum_{k=0}^{\infty}\phi^k\varepsilon_{t-k} $$ と書き直すことができます。また、MA過程は常に定常なので、MA過程に書きなおせるAR過程は定常です。ということで、AR過程が定常になる条件と、AR過程をMA過程に書きなおせる条件は同一になり、 AR特性方程式の解を調べることで判断ができます。

投稿