投稿

投稿日:

機械学習のプロダクト適用におけるリアルを語る会に参加しました

先日、ABEJAさんのオフィスで開催された、機械学習のプロダクト適用におけるリアルを語る会に参加しました。

AI/MLシステム開発の難しさ

発表者はディー・エヌ・エーの鈴木翔太さん。

企画からリリースまでの各段階において大変だったことを話されていました。
自分も経験があることが多くて全面的に納得です。

企画段階では、チームで何を作るのかのゴールの認識合わせが重要。
精度何%くらいを目指すのかと言った話もしておく必要があります。
また、全てを機械学習でやろうとしないことも大切ですね。

組織については、ML/AI開発チームにもソフトウェアエンジニアを入れて連携しているという話以前に、
MLチーム自体を作れてることを羨ましいと思いました。
(自分の職場には機械学習専門チームはなく、アナリストタスクと掛け持ちなので。)
懇親会で伺った話によるとやはりkaggler枠などの影響もあり、リファラルで採用が進んでいるそうです。

データ収集やアノテーション、データ管理が大変なのはどこでも語られている通り。
インフラをコスパ良く作るのも大変です。

また、誰が本番用コードを書くのかの話題も取り上げられました。
データサイエンティストが自分で書くのか、エンジニアに依頼するのか。
どちらも一長一短あります。
データサイエンティストのコードが汚いのは世界共通なんだろうか。

品質担保の話も取り上げられていたのは珍しいと感じました。
ただ、自分としてもここは悩んだ覚えがあります。
推論を含むシステムがエラーなく動くことをテストするのはいいとして、
予測である以上は結果が外れのこともあるわけで、
ここのデータに対して予測を外すのが正解みたいなテストを書くのも変で、
何がベストなのかよくわからなくなります。
そして最後はリリース後の評価やコストの話。

Deep Learningでリアルタイムに マーケット予測をしてみた

発表者はAlpaca Japanの梅澤慶介さん。

自分自身が個人トレーダーでもあるので密かに楽しみにしていたマーケット予測の事例。
ただ、やはり苦労されていたのはチームで開発するにあたって発生するあれこれのようでした。
メンバーの開発環境を再現できないとか、渡されたnotebookが動かないとか。
うちもここのメンバーのローカル環境揃ってないので良くわかります。
(自分一人に限ってもローカルPCで動いてたコードが開発用に建ててもらったインスタンスで動かないことがある。)
それをコンテナなどの技術で解決していったそうです。
Notebook禁止は衝撃。

機械学習のプラットフォームを自前で構築したり、
コードを共有してライブラリ化したりと、チームでの開発が順調に進んでいるようです。
ModelPackageみたいなの、僕も作りたいですね。

懇親会でも機械学習以外の部分含めて色々と質問ができ、非常に有意義でした。

智を集約しツラみを乗り越えたリピート推定の開発現場

発表者はABEJAの中川裕太さん。

題材はいくつかのミートアップですでに伺っているABEJA Insight for Retailのリピーター推定の話でした。
ただ、これまでと少し違った視点でのつらみの話があり興味深かったです。

「無限のビジネス可能性の前に立ち尽くす」ってのは新しい視点ですね。
たしかに、これができればいろんなことができそうっていうワクワクはありますが、
なんでもできるはなんにもできないと同じ、っていうのはあると思います。
(自分もかなり身に覚えがある)
顧客が片付けたいジョブに集中することが大事です。

研究と開発の時間軸が違って大変という話は、以前も伺った通りで
自分も直面するところなので良くわかります。
また、作ってる側から新しい手法が発見されるというも辛いところ。

コストと精度のトレードオフもありますね。
ここは、小売業向けのビジネスとしてのシビアな側面もありそうです。
(金融や医療、公的機関や大企業などお金があるところ向けでも無縁ではないと思いますが。)

今回はそれぞれの解決策もあり参考になりました。

マイクロサービスのDAGにして共通部分を使い回すなど、
実際、やってみるとかなり大変なのですがよく実現されたと思います。
スケーリングなどによるコスト最適化も、
アイデアは単純ですがおそらく実現は結構難しかったのではないかと。
精度の最適化の方にまだまだ課題も残っているようなので、今後に期待です。

まとめ

つらみの共有的な話を聞くことが増えていますが、
他社の人たちもそれぞれ苦戦しながら課題に立ち向かわれているのを知ると励みになりますね。
(普段は下手すると自分だけが効率悪く苦労しているような気がしてくるので。)

若干人数少なめな回だったのもあって、
懇親会でも発表者の方含めてゆっくり話ができて非常に有意義な時間を過ごせました。

機械学習の活用における課題は、
データサイエンティストだけで解決できるものだけではないと思うので、
この辺の知見はきちんと職場で共有して他チームの協力を得て取り組んでいきたいです。

投稿日:

Mac book に R環境構築

普段のデータ分析に使う環境はSQL, python, Tableau, Excelなどで完結させていたのですが、
業務上の都合でRも使うことになったので改めて学習を始めました。

Rを使うこと自他は初めてではなくpython使い始める前に何度か勉強したのですが、
どうにも好きになれなかったんで最近遠ざけていたのですがいよいよそうは言ってられなくなりましたね。

とりあえず、私物のMacにはインストールしていなかったので環境構築から始めます。
まずはRのインストールから。
配布サイトからダウンロードする方法もありますが、
homebrewを使う方法を採用しました。
いくつかのページを見ると、何かリポジトリを tapしないといけないように書いてあるところもあるのですが、
何もせずに brew install コマンドだけで入りました。


~$ brew install r
==> Installing dependencies for r: mpfr, gcc, libpng, openblas, pcre and readline
#(略) 依存パッケージとr本体がインストールされる

# 実行と終了
~$ r

R version 3.5.3 (2019-03-11) -- "Great Truth"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin18.2.0 (64-bit)

R は、自由なソフトウェアであり、「完全に無保証」です。
一定の条件に従えば、自由にこれを再配布することができます。
配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。

R は多くの貢献者による共同プロジェクトです。
詳しくは 'contributors()' と入力してください。
また、R や R のパッケージを出版物で引用する際の形式については
'citation()' と入力してください。

'demo()' と入力すればデモをみることができます。
'help()' とすればオンラインヘルプが出ます。
'help.start()' で HTML ブラウザによるヘルプがみられます。
'q()' と入力すれば R を終了します。

> q()

続いて、開発環境であるRStudioもインストールします。
これも配布元からダウロードする手もあるのですが、
homebrewで入れました。

これはアプリケーションフォルダに入れてGUIで使いたいので、brew cask を使います。


~$ brew cask install rstudio
==> Tapping homebrew/cask
Cloning into '/usr/local/Homebrew/Library/Taps/homebrew/homebrew-cask'...
remote: Enumerating objects: 4110, done.
remote: Counting objects: 100% (4110/4110), done.
remote: Compressing objects: 100% (4103/4103), done.
remote: Total 4110 (delta 24), reused 406 (delta 5), pack-reused 0
Receiving objects: 100% (4110/4110), 1.32 MiB | 187.00 KiB/s, done.
Resolving deltas: 100% (24/24), done.
Tapped 1 command and 4002 casks (4,116 files, 4.2MB).
==> Caveats
rstudio depends on R. The R Project provides official binaries:

  brew cask install r

Alternatively, the Homebrew-compiled version of R omits the GUI app:

  brew install r

==> Satisfying dependencies
==> Downloading https://download1.rstudio.org/desktop/macos/RStudio-1.2.1335.dmg
######################################################################## 100.0%
==> Verifying SHA-256 checksum for Cask 'rstudio'.
==> Installing Cask rstudio
==> Moving App 'RStudio.app' to '/Applications/RStudio.app'.
🍺  rstudio was successfully installed!

これで完了です。
念の為ですが、homebrew何も問題が起きてないことも見ておきました。


~$ brew doctor
Your system is ready to brew.
投稿日:

和分過程とARIMA過程

前回に引き続き今回も言葉の定義です。出典はいつもの沖本本から。

和分過程
d-1階差分をとった系列は非定常過程であるが、d階差分をとった系列が定常過程に従う過程は、
d次和分過程、もしくはI(d)過程と呼ばれます。
また、I(0)過程は定常過程で定義されます。

ARIMA過程
d階差分をとった系列が定常かつ反転可能なARMA(p,q)過程に従う過程は、
次数(p,d,q)の自己回帰和分移動平均過程、もしくはARIMA(p,d,q)過程と呼ばれます。

和分過程やARIMA過程の和分次数dは、AR特性方程式におけるz=1という解の個数に等しいことが知られているそうです。

前回の記事で定義した単位根過程は I(1)過程になり、
単位根過程の差分系列が定常で反転可能なARMA(p,q)過程の時、
それはARIMA(p,1,q)過程になります。

単位根過程の例として代表的なものにランダムウォークがあります。
言葉だけなら市場系のデータについて学んでいるとしょっちゅう出てきますね。

過程$y_t$が、
$$
y_t = \delta + y_{t-1} + \varepsilon_t , \ \ \ \varepsilon_t\sim iid(0, \sigma^2)
$$
と表現される時、$y_t$はランダムウォーク(random walk)と呼ばれます。
ただし、 $y_0=0$とします。また、 $\delta$はドリフト率と呼ばれるそうです。

定義式をよく見ると、AR(1)過程の形をしており、
AR特性方程式の解は$1$なので、これが定常でないこともわかります。
また、$\Delta y_t = y_t – y_{t-1} = \delta + \varepsilon_t$
であり、これは定常過程なので、ランダムウォークが単位根過程であることがわかります。

投稿日:

差分系列と単位根過程

久々に沖本先生の経済・ファイナンスデータの計量時系列分析に戻って、
時系列分析の用語を二つ紹介します。

まずは差分系列(5ページ)。
原系列$y_t$に対して、1点離れたデータとの差を取ることで生成できる系列、$y_t-y_{t-1}$を、
差分系列、または階差系列と呼び、$\Delta y_t$と書きます。

続いて単位根過程(105ページ)。
原系列$y_t$が非定常過程であり、差分系列$\Delta y_t=y_t-y_{t-1}$が定常過程である時、
過程は単位根過程である(unit root process)と言われます。

名前の由来は、単位根過程を誤差項が定常過程であるAR過程を用いて表現した時に、
AR特性方程式が$z=1$という解を１つ持つことらかきているそうです。

単位根過程の定義においては、差分系列が定常過程であることしか定義されておらず、
これがAR過程やARMA過程であることなどは要求されていないので注意です。

(お恥ずかしい話ですが、自分はデータサイエンティストに転職するずいぶん前に時系列分析を雑に勉強して、
差分とったらAR過程になるのが単位根過程だと誤って理解していたことがあります。)

このブログの通常の流れでは、ここから単位根過程のサンプルを一個構築して
プロットした図を出すのですが、実はすでに登場しているのでそちらにリンクしておきます。
参考:自己回帰過程のサンプル
この記事で構築したAR(3)過程の3個目の例が単位根過程です。

投稿日:

statsmodelsで重回帰分析

前回に引き続きpythonで重回帰分析を行い、回帰係数と切片を求める方法の紹介です。
今回はstatsmodelsを使います。

ドキュメントはこちら。
statsmodels.regression.linear_model.OLS

データは同じものを使い、結果が一致することを確認したいので
保存してたものを読み込みます。


import numpy as np
import statsmodels.api as sm

# データの読み込み
npzfile = np.load("sample_data.npz")
X = npzfile["X"]
y = npzfile["y"]

つぎに、回帰分析に入りますが、statsmodelsで回帰分析する場合には、一点注意があります。
それは、上で読み込んだ Xとyをそのままモデルに食わせると、切片(定数項)を0として結果を返してくることです。
それを回避するために、 X に対して、全ての値が1の列を追加する必要があります。
それは、次の関数を使って実行しています。
statsmodels.tools.tools.add_constant


# 配列Xに1だけからなる列を追加する。
X0 = sm.add_constant(X)
# 回帰分析実行
st_model = sm.OLS(y, X0)
results = st_model.fit()

# 切片と回帰係数
print(results.params)
# [ 4.94346432  2.00044153 -2.99255801  3.98231315  0.01708309]

無事にscikit-learnで行った時と全く同じ結果を得ることができましたね。
これだけだと、add_constantをやらなくてよかった分scikit-learnの圧勝に見えますが、
statsmodels を使うと一つメリットがあります。
それは、各係数が0でないことの検定ができることです。
resultsの sumary()って関数を実行してみましょう。


print(results.summary())

# 以下出力です。
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.999
Model:                            OLS   Adj. R-squared:                  0.999
Method:                 Least Squares   F-statistic:                 2.136e+04
Date:                Mon, 15 Apr 2019   Prob (F-statistic):          2.24e-139
Time:                        00:59:22   Log-Likelihood:                -143.94
No. Observations:                 100   AIC:                             297.9
Df Residuals:                      95   BIC:                             310.9
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9435      0.108     45.726      0.000       4.729       5.158
x1             2.0004      0.020    101.534      0.000       1.961       2.040
x2            -2.9926      0.018   -166.158      0.000      -3.028      -2.957
x3             3.9823      0.017    231.575      0.000       3.948       4.016
x4             0.0171      0.019      0.914      0.363      -0.020       0.054
==============================================================================
Omnibus:                        2.170   Durbin-Watson:                   1.994
Prob(Omnibus):                  0.338   Jarque-Bera (JB):                1.862
Skew:                          -0.208   Prob(JB):                        0.394
Kurtosis:                       2.477   Cond. No.                         7.12
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

出力の真ん中にある通り、定数項および回帰係数のt値と、p値が出力されています。
そして、x4 は P値が 0.05より大きいので、x4の係数が0であるという帰無仮説を棄却できていないことがわかります。
もともと、x4の回帰係数は0としてサンプルデータを作っていたので、これは良い結果です。
そのほか、AICなどの値も出してくれています。

scikit-learnでこれらの情報を収集するのは手間なので、
目的によってはstatsmodelsが便利な場面もあります。

投稿日:

scikit-learnで重回帰分析

今回と次回でpythonで重回帰分析を実行する方法を二つ紹介します。
今回はscikit-learnのLinearRegressionを使う方法です。

ドキュメントはこちら。
sklearn.linear_model.LinearRegression

最初に検証用のダミーデータを作ります。
$x_{i,j}$を -10 ~ 10の一様分布からサンプリングし、次の式で$y_i$を作ります。
$x_{i,3}$の係数が0になっていることから分かる通り、$x_{i,3}$は$y_i$には何の関係もない値です。
また、ノイズとして正規分布に従う乱数加えておきます。
$$
y_i = 5 + 2x_{i,0} -3x_{i,1} + 4x_{i,2} + 0x_{i,3} + \varepsilon_i,\\
\varepsilon_i \sim N(0,1) , \ \ \ i = 0,1,\cdots, 99
$$

サンプルデータを作って保存するコードがこちら。


import numpy as np
X = np.random.uniform(-10, 10, size=(100, 4))
y = 5 + X@[2, -3, 4, 0] + np.random.normal(0, 1, 100)
np.savez("sample_data.npz", X=X, y=y)

早速、回帰分析して回帰係数と定数項ついでに決定係数を求めてみましょう。
(回帰分析の目的が予測モデルを作ることであれば、データを訓練用と評価用に分けるべきなのですが、
今回は回帰係数を求める方法の紹介が主目的なので分けていません。)


import numpy as np
from sklearn.linear_model import LinearRegression

# データの読み込み
npzfile = np.load("sample_data.npz")
X = npzfile["X"]
y = npzfile["y"]

# モデルのインスタンス生成
model = LinearRegression()
#学習
model.fit(X, y)
# LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

# 回帰係数
print(model.coef_)
# [ 2.00044153 -2.99255801  3.98231315  0.01708309]

# 切片
print(model.intercept_)
# 4.943464324307898

# 決定係数 R^2
model.score(X, y)
# 0.9988893054003364

各コードの後ろにつけているコメントが僕の環境で実行した時の結果です。
回帰係数も切片もそれぞれほぼ正解に近い値が算出されていますね。

注:サンプルデータを乱数で生成しているので、データ生成からやり直せば結果は変わります。

scikit-learnを使うと、非常に手軽に回帰分析ができることがわかりました。
次回はstatsmodelsで同じことをやってみます。

投稿日:

numpyの配列をファイルに保存する

日常の実際で必要になったことはないのですが、
numpyに配列(ndarray)をファイルに保存する機能があるのでその紹介です。

(実用上、ファイルに保存したくなるようなデータは pandasのデータフレームに変換したあと、
to_csvやto_excelを使ってcsvかエクセルにしてしまうことが多いです。)

ドキュメントはこの辺り
numpy.save
numpy.load
numpy.savez

簡単に言ってしまえば、
numpy.save(ファイル名, 保存したい配列) で保存して、
numpy.load(ファイル名) で読み込める。
numpy.savez(ファイル名, (名前=)保存したい配列, (名前=)保存したい配列, …) で、
1ファイルに複数配列保存できる、と言った使い方になります。(名前=)は省略可能。

実際にやってみましょう。
まず検証用にダミーデータを作ります。


import numpy as np

# 保存のテスト用データを作成
data0 = np.arange(12).reshape(3, 4)
data1 = data0 ** 2
print(data0)
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(data1)
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

まずは、save関数で配列を一つ保存し、それを別の変数に読み込んでみましょう。
特に難しいことはなく非常に直感的な使い方で上手く動きます。


# .npyフォーマットで保存
np.save("data0.npy", data0)

# 読み込み
data0_load = np.load("data0.npy")
print(data0_load)
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''

続いて、savez関数で複数配列を1ファイルに保存し、復元してみます。
こちらは、loadした時に返されるオブジェクトが保存した配列そのものではないので注意です。
まずは名前をつけずに保存するパターン。


# savez を使うと複数の配列を.npzフォーマットの1ファイルにまとめて保存できる
np.savez("data0_1.npz", ary0, ary1)

# .npzファイルを読み込むと、配列ではなく、NpzFile オブジェクトが返される。
npzfile = np.load("data0_1.npz")
# .filesプロパティを見ると、中に二つの配列を持っていることがわかる
print(npzfile.files)
# ['arr_0', 'arr_1']

# それぞれの配列の取り出し
print(npzfile["arr_0"])
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(npzfile["arr_1"])
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

また、名前をつけて保存すると次のようになります。


# 名前(例ではxとy)をつけて保存することも可能
np.savez("named_data0_1.npz", x=ary0, y=ary1)
# .npzファイルを読み込むと、配列ではなく、NpzFile オブジェクトが返される。
npzfile2 = np.load("named_data0_1.npz")
# .filesプロパティを見ると、中に二つの配列を持っていることがわかる
print(npzfile2.files)
# ['x', 'y']

# 保存時の名前で取り出すことが可能
print(npzfile2["x"])
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(npzfile2["y"])
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

勝手にarr_0とか命名されるよりも、自分で名前をつけておいた方が混乱がなさそうですね。

投稿日:

記事数が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年経ちましたので、
仕事やキャリアの話なども、これから目指す人の参考になるように書いていきたいですね。

投稿日:

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に実装されているので同じように使えます。

投稿日:

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()

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