らくがき入門

環境と研究テーマが大幅に変わりました。だいたい何かに入門しています。

pandasを用いてフラグがついている列が先頭になるように行ごとにシフトする

やりたい処理は、すべての行でフラグ1が先頭にくるようにシフトしたい。

うまいやり方かどうかわからないけど、一応うまくいっているような気がする。

サンプルコード

import numpy as np
from numpy import nan
import pandas as pd
import warnings 
warnings.filterwarnings("ignore") # 実行に関係のない警告を非表示
# 左からゼロパディング
columns_list = ["a_{0:02d}".format(i) for i in range(0, 12)]
value_arr = np.array([[nan, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
                      [nan, nan, 1, 0, 0, nan, nan, nan, nan, nan, nan, nan], 
                      [nan, nan, nan, 1, 0, 0, 0, nan, nan, nan, nan, nan] 
                     ])

時系列のデータを想定したDataFrameを作成。

フラグが1になるカラムが行ごとにバラバラなので、すべての行で先頭のカラムが1になるようにシフトしたい。

df = pd.DataFrame(value_arr, columns=columns_list)
df
a_00 a_01 a_02 a_03 a_04 a_05 a_06 a_07 a_08 a_09 a_10 a_11
0 NaN 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1 NaN NaN 1.0 0.0 0.0 NaN NaN NaN NaN NaN NaN NaN
2 NaN NaN NaN 1.0 0.0 0.0 0.0 NaN NaN NaN NaN NaN

行ごとに処理を行う。

  1. 各行ごとにフラグ1がくるインデックスを取得して、その分だけ左方向にシフトする(作成されるのはSeries)

  2. 作成されたseriesを行方向に結合する

  3. 結合してできたDataFrameをtransposeして元のDataframeと形を揃える。

se_concat = pd.concat([df.ix[i].shift(val) for i, val 
                       in enumerate([-np.where(value_arr == 1)[1][j] 
                                     for j in range(len(df))])], axis=1)
se_concat.T
a_00 a_01 a_02 a_03 a_04 a_05 a_06 a_07 a_08 a_09 a_10 a_11
0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 NaN
1 1.0 0.0 0.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 1.0 0.0 0.0 0.0 NaN NaN NaN NaN NaN NaN NaN NaN

参考にした書籍・サイト

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

<Python, pandas> 縦にずらす。

numpyのarrayで複数の要素が配列内に存在するか判定する

numpy.arrayで複数の要素をリストで渡して、真偽値行列を作成してみる。

使うのは、numpy.in1dで、 指定したarray-likeな要素がある配列内に存在するかどうかを判定して、1次元の真偽値行列を返してくれる。

配列の形を合わせたいなら、reshape(arr.shape)で形を合わせることができる。

In [1]: import numpy as np

In [2]: arr = np.array([[1, 3, 5], [2, 3, 4], [2, 3, 3]])

In [3]: np.in1d(arr, [1, 3, 5])
Out[3]: array([ True,  True,  True, False,  True, False, False,  True,  True])

In [4]: np.in1d(arr, [1, 3, 5]).reshape(arr.shape)
Out[4]: 
array([[ True,  True,  True],
       [False,  True, False],
       [False,  True,  True]])

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

LightGBMをインストールする

LightGBMをubuntu18.04にインストールします。

LightGBMとは

Microsoftが開発した勾配ブースティング(Gradient Boosting)のフレームワーク。勾配ブースティングのフレームワークには、他にXGBoostとかも割と有名なのは知っていましたが、LightGBMは知りませんでした。知ったきっかけは、Kaggleで適当に興味あるコンペのKernelを読み漁っているとよく登場してたため調べてみる気になりました。

勾配ブースティング

そもそもブースティングは、複数の弱学習器を用意して、それぞれの学習器を直列接続するイメージの手法。前の弱学習器で学習した内容を現在の弱学習器に継承しながら学習を進めていきます。

勾配ブースティングは、各ステップごとの弱学習器のでの損失関数の最小化問題に対して勾配降下法を用いるのが由来です。

勾配ブースティングの弱学習器では、決定木が採用されることが多いです。

LightGBMのインストール

LightGBMの公式の通りインストールしていくだけですが、一応手順を示します。 下準備としてsetuptoolsは準備しておきましょう。 今回は、Github経由でのインストールを採用しました。

Cmakeのインストール

Cmakeはソフトウェアのビルド自動化ツールで、Windows,mac, linuxとクロスプラットホームで使えます。 LightGBMのインストール時に必要なのでインストールします。

$ cd ~/Downloads
$  wget https://cmake.org/files/v3.12/cmake-3.12.2.tar.gz
$ tar xvf cmake-3.12.2.tar.gz
$ cd cmake-3.12.2
$ ./configure
$ make
$ sudo make install 
# インストールできているか確認
$ cmake
Usage

  cmake [options] <path-to-source>
  cmake [options] <path-to-existing-build>

Specify a source directory to (re-)generate a build system for it in the
current working directory.  Specify an existing build directory to
re-generate its build system.

Run 'cmake --help' for more information.

cmakeして上の表示がされていれば、インストールできています。 うまく行かない場合は、C++コンパイラがない可能性があるのでgccとかを入れる必要があるかも?

LightGBMのインストール

$ git clone --recursive https://github.com/Microsoft/LightGBM.git
$ LightGBM/python-package
$ python setup.py install
# grepしてインストールできているか確認
$ pip list --format=columns | grep -i lightgbm
lightgbm                           2.2.1 

上のように、lightgbm 2.2.1 のようにversion表示されていればインストール成功です。

参考にしたサイト

qiita.com

blog.amedama.jp

はじめてのパターン認識

はじめてのパターン認識

2つのnumpy arrayからDataFrameを作る

2つの同じ長さのnumpy arrayを用いてそれらを列に持つpandasのDataFrameを作成します。

In [1]: import numpy as np

In [2]: import pandas as pd

In [3]: arr1 = np.array(["a", "b", "c", "d"])

In [4]: arr2 = np.array([1, 2, 3, 4])

In [5]: dict_ = dict(zip(arr1, arr2))

In [6]: df = pd.DataFrame(list(dict_.items()), columns=["col1", "col2"])

In [7]: df
Out[7]: 
  col1  col2
0    a     1
1    b     2
2    c     3
3    d     4

手順は、

  1. 1つのarrayをkey、もう一つのarrayをvalueとする辞書を作成
  2. 作成した辞書をlist.items()でkeyとvalueのarrayのタプルを1つの値ごとに格納したリストに変換
  3. 2.で作成したリストをDataFrameに変換

という流れです。

参考にしたサイトは以下のサイトです。

stackoverflow.com

PandasのDataFrameから特定の値を持つ行を削除する

言われてみたら簡単なんだけど、意外に思いつかなかった。

DataFrameの特定の行に含まれている値を指定して、それ以外を抽出するイメージ。

In [1]: import pandas as pd

In [2]: df = pd.DataFrame([[1, 2], [2, 3], [3, 4]], columns=["a", "b"])

In [3]: df
Out[3]: 
   a  b
0  1  2
1  2  3
2  3  4

In [4]: df[df.a != 2]
Out[4]: 
   a  b
0  1  2
2  3  4

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 第2版 ―NumPy、pandasを使ったデータ処理

pythonで地図上の2時点間の位置関係を求める;

2時点間の緯度、経度、高度が与えられたときに簡易的に2時点間のローカルな位置関係を求めます。 緯度・経度・高度といった情報ではローカルな2時点間の関係性が分かりづらく扱いづらいので、変換します。 具体的には、2時点間の距離とある時点からもう一時点を見たときの方位角、仰角を求めることを目標とします。

以下、国土地理院のサイトと独立行政法人電子航法研究所の資料を参考に作成しています。

https://www.enri.go.jp/~fks442/K_MUSEN/1st/1st060428rev2.pdf

地球上の2時点間の距離

高校などで扱った平面世界での2時点間の距離は、簡単に求めることができたと思いますが、地球上の2時点間ではそれほど簡単には求めることができません。 なぜなら、地球が球体であり、局面上にある2時点間の距離を求めることになるからです。

ジオイドとは

地球は、時点による遠心力の影響で、赤道方向が少し膨らんだ楕円体のような形をしています。 測地学では、世界の界面の平均位置に最も近い「重力の等ポテンシャル面」を「ジオイド」と定めて、これを地球の形状ということにしています。標高はジオイドを基準として定められており、標高はジオイドから測った高さになります。

GPSを用いて標高を求めるには

現在広く活用されているGPSでは、緯度・経度・高度といった幾何学的な位置を求めることができますが、重力を考慮した設計になっていないため標高を直接求めることができません。GPSを用いて標高を求めるためには、ジオイド高がj必要になります。

標高 = 楕円体高 ー ジオイド

ジオイド高は国土地理院が提供している計算サイトで計算することができます。

ジオイド計算

WGS84

WGS84とは、GPSの基準座標系で、直交座標系です。これはECEFと呼ばれ、地球の重心を原点として地球の自転軸の北極方向をZ軸、Z軸に垂直にグリニッジ子午線の方向をX軸として、これらの軸と直行するように右手系でY軸とした直交座標系です。

ECEFから東をX, 北をY、上空をZとした地平直交座標(ENU)に変換します。 ECEFからENUへの変換は、行列による回転と原点移動で実現できます。詳しくは上記で示した電子航法研究所の資料をご確認ください。

githubにサンプルコードを提示してあります。coordinate_test.pyの9, 10行目で取り扱いたい2時点の緯度・経度・ジオイド高のサンプルデータを与えています。結果表示では、2時点間の距離およびある時点(コードではRx)からみたもう1時点(コードではSat)の方位角、仰角を求めることができます。

github.com

時系列データにおける定常過程と単位根過程の違いとADF検定(Pythonによる実装例あり)

前回の記事では、「時系列データの定常性を確認する」というタイトルで記事を書きました。与えられた時系列データが定常過程であれば、ARMAモデルなどの基礎的なモデルに落とし込むことができます。

今回はよく用いられる単位根検定である拡張 Dickey-Fullar(augmented Dickey-Fullar test;ADF test)を取り上げます。

使用言語はPythonです。 今回も沖本本を参照しています。

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

単位根過程とは

前回の記事で、定常過程に関しては取り上げているのでそちらを参照してください。

leisurelab.hatenablog.com

単位根過程とは、原系列 y_tが非定常過程であり、差分系列 \Delta y_t = y_t - y_{t-1}が定常過程であるとき、過程は単位根過程であるといわれる。

単位根過程は、別名がいくつかあり、差分系列が定常となるため差分定常過程と呼ばれることがあります。

さて今回Pythonのコード中の検証データとして使用しているランダムウォークに関しても記述しておきます。 なぜ今回ランダムウォークを検証データとして使用しているかと言うと、ランダムウォークが単位根過程の代表選手だからです。 ランダムウォークの定義は、

過程 y_t
 \displaystyle
y_t = \delta + y_{t-1} + \varepsilon_t
,  \varepsilon_t \sim \text{iid}(0, \sigma ^2)
と表現されるとき、 y_tランダムウォークと呼ばれる。定数項 \deltaはドリフト率と呼ばれる。

上記の定義より、繰り返し代入すると、
 \displaystyle
y_t = \delta t + \nu_t
と表現できます。ただし、 \nu_t = \varepsilon_1 + \varepsilon_2 + \cdots + \varepsilon_tという撹乱項の線形和で表現され、これは確率的トレンドと呼ばれます。 ドリフト率 \deltaは線形トレンドの傾きを表しているので、単位根過程は線形トレンドと確率的トレンドを併せ持つ過程と言えます。

定常過程と単位根過程との性質の違い

定常過程と単位根過程の違いを以下の表にまとめました。

性質 定常過程 単位根過程
トレンド もたない 線形トレンドを持つ
長期予測 過去の観測値が減衰していき過程の期待値に近づく 過去の観測値の影響が消えない
予測誤差(MSE) 予測期間が長くなるにつれて過程の分散(有限値)に近づく 長期的に線形的に増大していく
インパルス応答関数 ショックは一時的な影響力を持つ ショックは恒久的な影響力を持つ

表のトレンドに関して詳しく説明します。
定常過程(弱定常性に従う過程)は前回の記事で説明したように、期待値と自己共分散が時間によらず一定であるという仮定を置いています。 この仮定より、明らかに定常過程はトレンドを持ちません(時間によらず一定であるから)。 定常過程は平均回帰性を持ち、過程が長期的には必ず平均の方向に戻っていくことを意味します。
ただし、定常過程であっても線形トレンドを記述できるモデルとして、トレンド定常過程と言われるものが存在します。
トレンド定常過程は、文字通り定常過程にトレンド項を加えることで線形トレンドを記述できるようにしたモデルを指します。
それでは、トレンドを持つトレンド定常過程と単位根過程は線形トレンドを記述できるという点で同じなわけですが、前述したように単位根過程は確率的トレンドも持つ点で異なります。 単位根過程は不確実性を線形的に増大させていると言えます。

インパルス応答関数に関してですが、VARモデルに関して取り上げる際に書こうと思っているので一旦保留とします。 簡単に概要だけ説明すると、ある時点の影響が将来どの程度影響をもつかを定量的に表現できるツールです。

単位根検定

単位根検定として有名なADF検定について説明します。 ADF検定は真の過程をAR(p)モデルと仮定し、仮定が単位根AR(p)過程であるという帰無仮説を、過程が定常AR(p)過程であるという対立仮説に対して検定するものです。 これが厳密な定義ですが、帰無仮説「単位根がある」、対立仮説「単位根がない(=定常である)」と言ってしまってもいいと思います。

この検定を用いる場合、ある有意水準に対して帰無仮説を棄却できれば、対立仮説が支持されて与えられた時系列データが定常過程であると言って良さそうです。 ただし、帰無仮説を棄却できなかった場合は、帰無仮説を積極的に支持したことにはならないことに注意してください。これは統計学的非対称性によるのですが長くなりそうなのでまた別の機会に取り上げます。

したがって、与えられた時系列データが単位根過程であるかは、冒頭で取り上げた単位根過程の定義である過程の差分系列が定常課程であることを用います。 したがって、手順としては以下のようになります。

  1. 与えられたデータに対してADF検定→
    帰無仮説が棄却されれば、与えられたデータは定常過程
  2. 帰無仮説を棄却できなければ差分系列を求め、その系列に対してADF検定→
    帰無仮説が棄却できれば、与えられたデータは単位根過程
  3. 帰無仮説が棄却できなければ、与えられたデータは単位根過程ではない非定常過程

このようにしてデータを分類することになります。

この手順をPythonで実装した例を示します。

gist7fb5ceb03fae95dc3f0774debf880340