標準化レーベンシュタイン距離は距離関数なのか

以前の記事で、標準化レーベンシュタイン距離(標準化編集距離)というのを紹介し、
自分も使っていたのですが挙動に少し違和感があったので確認しました。

参考:標準化レーベンシュタイン距離

レーベンシュタイン距離はその名の通り、距離関数なのですが、
これを標準化してしまうとどうも距離関数っぽくない動きをしてるように思えたのです。

念の為、距離関数というもの自体の定義をおさらいしておきましょう。

集合$X$に対して、$d:X\times X \rightarrow \mathbb{R}$ が距離関数であるとは、
$x,y,z \in X$に対して次の条件が成り立つ時に言います。
1. $d(x,y) \geq 0$ (非負性)
2. $d(x, y) = 0 \Leftrightarrow x = y$ (同一律)
3. $d(x, y) = d(y, x)$ (対象律)
4. $d(x, z) \leq d(x, y) + d(y, z)$ (三角不等式)

この条件のうち、 1. 2. 3. は特に問題ないのですが、
標準化レーベンシュタイン距離 については、4. の三角不等式がちょっと怪しかったです。
で、反例を探してみたところ簡単に見つかりました。
$ld(*,*)$をレーベンシュタイン距離、$nld(*,*)$を標準化レーベンシュタイン距離とし、
x = ‘ab’, y= ‘aba’, z = ‘ba’ とおきます、
すると、
$ld(x, z) = 2$ なので、$nld(x, z) = 1$ ですが、
$ld(x, y) = ld(y, z) = 1$ なので、$nld(x, y) = nld(y, z) = \frac13$ です。

そのため、 $nld(x, z) > nld(x, y) + nld(y, z) = \frac23$ となり、
三角不等式を満たしません。

標準化レーベンシュタイン距離 は 標準化レーベンシュタイン という名前の距離関数と考えるのは誤りで、
レーベンシュタイン距離 という距離関数を標準化したもの(その結果距離関数ではなくなってしまったもの)と、
考える必要があります。

現状これで激しく困ったということはないのですが、
一部のライブラリにある、自分で作った距離関数を引数に渡せるようなものには、
標準化レーベンシュタイン距離は突っ込まない方が安全そうです。

pandasで指数平滑移動平均

昨日の記事が移動平均だったので、今日は指数平滑移動平均を扱います。
初めて知った日は衝撃だったのですが、pandasには指数平滑移動平均を計算する専用の関数が用意されています。
(pythonを使い始める前はExcel VBAでいちいち実装していたので非常にありがたいです。)

馴染みがない人もいると思いますので軽く紹介しておきます。
元のデータを${x_t}$とし、期間$n$に対して指数平滑移動平均${EWMA_t}$は次のように算出されます。
$$
\begin{align}\alpha &= \frac{2}{1+n}\\
EWMA_0 &= x_0\\
EWMA_t &= (1-\alpha)*EWMA_{t-1} + \alpha * x_t
\end{align}
$$

3番目の式を自分自身に逐次的に代入するとわかるのですが、
$EWMA_t$は、$x_t$から次のように算出されます。
$$
EWMA_t = \alpha\sum_{k=0}^{\infty}(1-\alpha)^k x_{t-k}
$$
$(1-\alpha)$の絶対値は1より小さいので、この無限級数の後ろの方の項は無視できるほど小さくなります。
結果的に、過程${x_t}$の最近の値に重みを置いた加重平均と見做せます。

さて、早速ですが計算してみましょう。

pandasのDataFrameおよび、Seriesに定義されているewm関数を使います。
pandas.DataFrame.ewm


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# データ作成
data = pd.Series(np.random.normal(0, 100, 200).cumsum() + 20000)
# 指数平滑移動平均の計算
data_ewm = data.ewm(span=10).mean()
# 可視化
plt.rcParams["font.size"] = 14
fig = plt.figure(figsize=(12, 7))
ax = fig.add_subplot(1, 1, 1)
ax.plot(data, label="元データ")
ax.plot(data_ewm, label="指数平滑移動平均")
plt.legend()
plt.show()

出力がこちら。

ここで一つ注意する点があります。
data_ewm = data.ewm(span=10).mean()
という風に、spanという変数名で期間$10$を渡しています。
ドキュメントを読んでいただくとわかるのですが、span=をつけないと、
comという別の変数に値が渡され、$\alpha$の計算が、
$\alpha=1/(1+com)$となり、結果が変わります。

また、spanやcomを使う以外にも、alpha=で$\alpha$のあたいを直接指定することも可能です。

pandasで移動平均や高値線、安値線を計算する

前回がローソク足だったので今回も市場データでよく使われるテクニックから移動平均を取り上げてみたいと思います。
ついでにHLバンド(ドンチャンチャンネル/高値線,安値線)も同様にもとまるので紹介します。

技術としては、window関数と呼ばれる種類の関数を使って算出します。

ドキュメントはこの辺り。
Window
pandas.DataFrame.rolling
pandas.Series.rolling

DataFrameとSeries両方に実装されていて、同じように使うことができます。
rolling() で 指定期間ごとに区切ったデータを作り、それに対して、 meanやmax,minなどの関数を適用して
平均や最大値、最小値を算出して配列として返します。

実際に見た方が早いと思うので、乱数でランダムウォークデータを生成し、
計算して可視化してみましょう。


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# データ作成
data = pd.DataFrame(np.random.normal(0, 100, 200).cumsum() + 20000)
# 移動平均、期間の高値/安値線の計算
ma = data.rolling(10).mean()
high_band = data.rolling(20).max()
low_band = data.rolling(20).min()
# 可視化
fig = plt.figure(figsize=(12, 7))
ax = fig.add_subplot(1, 1, 1)
ax.plot(data, label="元データ")
ax.plot(ma, label="移動平均")
ax.plot(high_band, label="高値線")
ax.plot(low_band, label="安値線")
plt.legend()
plt.show()

結果がこちら。

期間の最初の方はデータ不足により線が途切れています。
この辺りの制御は min_periods などの引数で細かく調整できるので、
データ量やその時の目的に応じて調整して使いましょう。

pythonでローソク足を描く

以前(pythonを勉強し始めた頃)は、matplotlibでローソク足をかけたはずなのですが、最近は方法が変わってしまったようなのでそのメモです。
なお、ここでサンプルに使うデータはすでにcsvフィアルか何かで保存されているものとします。

以前は matplotlib.finance というのをimport でき、これを使ってかけたのですが、version 2.0 からなくなってしまったようです。
matplotlib.finance
This module is deprecated in 2.0 and has been moved to a module called mpl_finance.

そしてさらに良くないことに、移動先の mpl_finance ですが、あまりしっかり保守されてない様子。

githubのリポジトリに下記の文言があります。
The code is provided as is and is basically un-maintained.

ただ、一応動くようなので動かしてみましょう。
anacondaには含まれていないようなので、インストールから必要です。
pip install mpl_finance
これで、
mpl-finance==0.10.0
が入りました。

さて、使い方ですがun-maintainedの宣言通り、 mpl-finance の公式ドキュメントらしきものは見当たらず、
上の、matplotlib.finance時代のドキュメントを読んで使わないといけないようです。

ローソク足を書く関数は次の4つあり、それぞれデータの渡し方が違います。
.candlestick2_ochl(ax, opens, closes, highs, lows, width=4, colorup=’k’, colordown=’r’, alpha=0.75)
.candlestick2_ohlc(ax, opens, highs, lows, closes, width=4, colorup=’k’, colordown=’r’, alpha=0.75)
.candlestick_ochl(ax, quotes, width=0.2, colorup=’k’, colordown=’r’, alpha=1.0)
.candlestick_ohlc(ax, quotes, width=0.2, colorup=’k’, colordown=’r’, alpha=1.0)

今回は手元のデータと相性が良いので .candlestick_ohlc を使います。
quotes に 日付、始値、高値、安値、終値、の5列のデータがデータ件数行だけ並んだ配列を渡してあげる必要があります。
ここで面倒なのは日付の渡し方で、float型で渡す必要があります。
ドキュメントに time must be in float days format – see date2numとある通り、専用の関数があるのでそれを使います。

matplotlib.dates.date2num(d)

また、この関数は引数がdatetime型なので、元々が2019-05-07 のような文字列になっているならば、
datetime型に変換しておく必要があります。
それにはpandasの to_datetimeを使います。
pandas.to_datetime
(いつもならそれぞれ1記事使ってるようなテクニックですね。to_datetimeの方は便利なのでそのうち専用記事書くかも。)

前置きが長くなりましたが、ここまでの情報でできるので日経平均のcsvデータからローソク足を書いてみましょう。


import pandas as pd
import mpl_finance
import matplotlib.pyplot as plt
from matplotlib.dates import date2num

# データの読み込み
df = pd.read_csv("./日経平均データ.csv")
print(df.head())

'''
        date      open      high       low     close
0  2019-3-11  21062.75  21145.94  20938.00  21125.09
1  2019-3-12  21361.61  21568.48  21348.81  21503.69
2  2019-3-13  21425.77  21474.17  21198.99  21290.24
3  2019-3-14  21474.58  21522.75  21287.02  21287.02
4  2019-3-15  21376.73  21521.68  21374.85  21450.85
'''

# dateの型変換
# まずdatetime型にする
print(df["date"].dtypes)  # object
df["date"] = pd.to_datetime(df["date"])
print(df["date"].dtypes)  # datetime64[ns]
# 続いて float型へ
df["date"] = matplotlib.dates.date2num(df["date"])
print(df["date"].dtypes)  # float64

print(df.head())
'''
       date      open      high       low     close
0  737129.0  21062.75  21145.94  20938.00  21125.09
1  737130.0  21361.61  21568.48  21348.81  21503.69
2  737131.0  21425.77  21474.17  21198.99  21290.24
3  737132.0  21474.58  21522.75  21287.02  21287.02
4  737133.0  21376.73  21521.68  21374.85  21450.85
'''

# 可視化
fig = plt.figure(figsize=(13, 7))
ax = fig.add_subplot(1, 1, 1)
mpl_finance.candlestick_ohlc(ax, df.values)
plt.show()

こうして出来上がるチャートが次です。

正直これ単体では手間の割に可視化するメリットがないなーというのが正直なところです。
ただ、matplotlibの仕組みに乗っかっているので、
自分のオリジナルの指標などを追加していくことができます。

matplotlibのpcolorとpcolormesh

先日のトピックモデルの記事中で、試しにヒートマップでの可視化を試みた時、
matplotlibのpcolorって関数を使用しました。

参考:pythonでトピックモデル(LDA)

matplotlibでヒートマップを描こうと思うと少々無理やりな実装になるものも含めて、
imshowや、contourf、pcolorなど複数の方法が考えられ、結構迷いますが一番自然に書ける気がして採用しました。
しかしどうやらこのpcolor、あまり評判がよろしくないようです。

公式ドキュメントを見ても次ように記載があります。

matplotlib.axes.Axes.pcolor

Hint

pcolor() can be very slow for large arrays. In most cases you should use the similar but much faster pcolormesh instead. See there for a discussion of the differences.

要するに pcolormesh を使う方が良いようです。

ドキュメントはこちら。
matplotlib.axes.Axes.pcolormesh

Differences between pcolor() and pcolormesh()
の部分を読んでも、あまり pcolorにメリットを感じないので、
いっそのこと pcolor 自体を pcolormesh のエイリアスか何かに変えてしまっても良さそうなのですが、
戻り値の型が違うこともありますし、何か同じように使えない事情もあるようですね。

とりあえず今後はこれまでpcolorを使っていた場面ではpcolormeshを使うようにしようと思います。
ちなみに、x軸y軸が等間隔の場合はimshowの方がさらに速いそうです。

imshow
If X and Y are each equidistant, imshow can be a faster alternative.

標準化レーベンシュタイン距離

以前の記事で、レーベンシュタイン距離を計算できるライブラリを紹介しました。
参考:pythonで編集距離(レーベンシュタイン距離)を求める

どちらかというと、ライブラリよりもアルゴリズム側の説明の続きなのですが、
標準化されたレーベンシュタイン距離(normalized Levenshtein distance)というものも提案されています。
これは、二つの文字列のレーベンシュタイン距離を、文字数が多い方の文字数で割った値として定義されます。
固有名詞の名寄せなどでレーベンシュタイン距離を使う場合、
こちらを使った方がうまく行くことが多いようです(個人的な経験から。)

以前紹介した、 python-Levenshtein
実装されているんじゃないかと期待してしばらく調べていたのですが、どうやらこれは実装されてないようです。
特にLevenshtein.ratio という関数に期待したのですがこれは全然違いました。

ということで自分で実装しましょう。


import Levenshtein


def normalized_distance(text0, text1):
    dist = Levenshtein.distance(text0, text1)
    max_len = max(len(text0), len(text1))
    return dist / max_len

関数名は normalized_levenshtein_distance にしたかったが流石に長すぎるので少し短縮。
ただ、これでも長いですね。(pep8的につらい)

これによって、例えば通常のレーベンシュタイン距離では、
“バニラ” と “アイス” の組み合わせも “チョコレート” と “チョコレートアイス” もどちらも距離は3でしたが、
標準化した距離を使うことで、
前者は距離1、後者は 1/3(=0.333…)と算出されるようになり、前者の方が離れてると見なせます。

使ってみた結果がこちら。


print(Levenshtein.distance("アイス", "バニラ"))  # 3
print(Levenshtein.distance("チョコレートアイス", "チョコレート"))  # 3
print(normalized_distance("アイス", "バニラ"))  # 1.0
print(normalized_distance("チョコレートアイス", "チョコレート"))  # 0.3333333333333333

scipyで距離行列を計算する

前の記事でちらっと pdist関数が登場したので、scipyで距離行列を求める方法を紹介しておこうと思います。

距離行列の説明はwikipediaにあります。
距離行列 – Wikipedia

要するに、N個のデータに対して、(i, j)成分がi番目の要素とj番目の要素の距離になっているN*N正方行列のことです。

これは次の二つの関数を使って計算できます。
scipy.spatial.distance.pdist
scipy.spatial.distance.squareform

numpyで適当に5点取ってやってみましょう。


import numpy as np
from scipy.spatial.distance import pdist
from scipy.spatial.distance import squareform

# 出力する桁数を抑える
np.set_printoptions(precision=3)
# 乱数生成
X = np.random.randint(-5, 6, size=(5, 2))
print(X)
'''
[[-2 -4]
 [-3 -4]
 [ 2 -1]
 [ 4 -2]
 [-1 -2]]
'''

y = pdist(X)
print(y)
'''
[1.    5.    6.325 2.236 5.831 7.28  2.828 2.236 3.162 5.   ]
'''

M = squareform(y)
print(M)
'''
[[0.    1.    5.    6.325 2.236]
 [1.    0.    5.831 7.28  2.828]
 [5.    5.831 0.    2.236 3.162]
 [6.325 7.28  2.236 0.    5.   ]
 [2.236 2.828 3.162 5.    0.   ]]
'''

変数Mに距離行列が入っているのが確認できますね。
(1,2)成分が [-2 -4] と [-3 -4] の距離の1,
(1,3)成分が [-2 -4] と [ 2 -1] の距離の5など、きちんと距離が入っています。

なお、 pdistの戻り値自体は、正方行列の形をしておらず、squareformで正方行列に変形する必要があるので注意です。
pdist の戻り値(コード中では変数y)には、上三角行列部分の値が、1行めから順番に入っています。
配列の長さは$(N*(N-1)/2) = 5*4/2 = 10$ です。

pdist 関数は metric という引数に様々な距離関数の値をとることができ、いろんな種類の距離行列を作ることができます。
(デフォルトはユークリッド距離)

ここまで紹介した、pdistは1つのデータセットに対して距離行列を生成しましたが、
cdistという関数を使うと、2個のデータセットからそれぞれ1個ずつ要素を抽出して距離行列を作ることもできます。

こちらはsquareformを使わなくても初めから行列の形で結果が得られます。


from scipy.spatial.distance import cdist
XA = X[:3]
XB = X[3:]
print(cdist(XA, XB))

# 出力
'''
[[6.325 2.236]
 [7.28  2.828]
 [2.236 3.162]]
'''

最初の例においても、
pdistとsquareform を使うよりも、 cdist(X, X) した方が便利かもしれません。
結果は同じになります。

pandasでカテゴリー変数を数値に変換する

pandasの関数を使ってカテゴリー変数を数値に変換する方法のメモです。
one-hot表現とは異なり、[“a”, “b”, “c”] みたいなデータを、 “a”は0,”b”は1,”c”は2みたいな感じで数字に変換します。

scikit-learnのデータセットにおける、正解ラベルのようなデータ型にする、と言ったほうがわかりやすいかもしれません。
Rにはカテゴリカル変数の型があるのでそちらの方がイメージしやすい人もいるかも。

使う関数は、pandas.factorizeです。

この関数に変換したいリストを渡すと、数値に変換したリストと、何番めの数がもともと何の値だったのかを示す配列を返してくれます。

実際にやってみましょう。
numpyでランダムに10個のデータを作って、それを数値化します。


import numpy as np
import pandas as pd

data = np.random.choice(["a", "b", "c"], 10)
print(data)
# ['b' 'c' 'a' 'c' 'c' 'c' 'c' 'a' 'b' 'a']
labels, uniques = pd.factorize(data)
print(labels)
# [0 1 2 1 1 1 1 2 0 2]
print(uniques)
# ['b' 'c' 'a']

うまくいきました。
uniques の値をアルファベット順にしたい時は、
factorize する時に sort=Trueを指定するとできます。


labels, uniques = pd.factorize(data, sort=True)
print(labels)
# [1 2 0 2 2 2 2 0 1 0]
print(uniques)
# ['a' 'b' 'c']

できました。

numpyで配列内の値を特定の値に制限する

前の記事の最後で、
異常時の前処理として、1〜99パーセンタイルでクリップするって話を少し書いたので、
それをnumpyで実現する関数の紹介です。
と言っても、わざわざ専用関数を使わなくても容易に実装できるのですが、せっかく用意されているのがあるので知っておくと便利です。

ドキュメントはこちら。
numpy.clip

np.clip(配列, 最小値, 最大値)と指定すると、
配列の値のうち、区間[最小値, 最大値]からはみ出た値を、その範囲に収まるように区切ってくれます。

ためしに、標準正規分布に従う値20個を生成して、[-1, 1]の範囲にクリッピングしてみましょう。


import numpy as np
data = np.random.randn(20)
print(data)
'''
[-1.71434343  0.33093523 -0.0648882   1.34432289 -0.15426638 -1.05988754
 -0.41423379 -0.8896041   0.12403786  1.40810052  0.61199047  1.60193951
 -0.72897283 -0.00861939 -0.38774556  0.40188148  0.08256356  1.61743754
 -0.12320721  1.45184382]
'''
data = np.clip(data, -1, 1)
print(data)
'''
[-1.          0.33093523 -0.0648882   1.         -0.15426638 -1.
 -0.41423379 -0.8896041   0.12403786  1.          0.61199047  1.
 -0.72897283 -0.00861939 -0.38774556  0.40188148  0.08256356  1.
 -0.12320721  1.        ]
'''

-1 より小さかった値は-1に、 1より大きかった値は1になりました。

ちなみに下記のコードでも同じことができます。


data[data < -1] = -1
data[data > 1] = 1

前回紹介した、percentileと組み合わせて使うことで、
nパーセンタイルからmパーセンタイルにクリップするということも簡単に実現できます。

試しに 5〜95パーセンタイルにクリップしたのが次のコードです。


data = np.random.randn(10)
print(data)
'''
[-0.41127091 -1.34043164  0.09598778 -1.19662011 -0.04607188 -0.02745831
  0.23184919  0.85601106  0.58430572  0.88205005]
'''
c_min, c_max = np.percentile(data, [5, 95])
print(c_min, c_max)
'''
-1.2757164503037743 0.8703325037861378
'''
data = np.clip(data, c_min, c_max)
'''
[-0.41127091 -1.27571645  0.09598778 -1.19662011 -0.04607188 -0.02745831
  0.23184919  0.85601106  0.58430572  0.8703325 ]
'''
print(data)

この例だけ見てもありがたみを感じないのですが、実際のデータを決定木などにかける時、
ほんの数件のデータだけ極端な外れ値になっていたりすると、
いい感じの範囲にデータを収めることができるので便利です。

また、scikit-learnなどのライブラリのコードを見てみると、
値を 0より大きく1より小さい範囲に収める目的などでも使われています。
ここなど
n以上m以下、ではなくnより大きいmより小さい、で区切る時は便宜上、eps=1e-15のような非常に小さい値を用意して、
[n+eps, m-eps]で代用するようですね。
こういう書き方も参考になります。

numpyのpercentile関数の仕様を確認する

中央値や四分位数を一般化した概念に分位数ってのがあります。
その中でも特にq/100分位数をqパーセンタイルといい、numpyに専用の関数が用意されています。
numpy.percentile

データの可視化や外れ値の除外で使うためにこれの仕様を確認したのでそのメモです。

そもそも僕が何を疑問に思ったのかを説明したほうがいいと思うので、いくつか例を紹介します。

まずわかりやすい例で50パーセンタイル。
これは、奇数個の値があればその中央の値、偶数個の値に対しては、真ん中の二つの値の中点を返します。


import numpy

# 5個の値の3番目の数を返す
data_1 = np.array([3, 12, 3, 7, 10])
print(np.percentile(data_1, 50))  # 7.0

# 6個の値の3番目の数と4番目の数の平均を返す
data_2 = np.array([3, 12, 3, 7, 10, 20])
print(np.percentile(data_2, 50))  # 8.5

同様にして、区切りのいい値がある時のパーセンタイルは非常にわかりやすい。
11個の値があれば、それぞれ順番に 0パーセンタイル, 10パーセンタイル, … 90パーセンタイル, 100パーセンタイルです。


data_3 = np.random.randint(0, 2000, 11)
print(data_3)
# 出力
# [1306  183 1323  266  998 1263 1503 1986  250  305 1397]
for p in range(0, 101, 10):
    print(p, "パーセンタイル・・・", np.percentile(data_3, p))
# 出力
'''
0 パーセンタイル・・・ 183.0
10 パーセンタイル・・・ 250.0
20 パーセンタイル・・・ 266.0
30 パーセンタイル・・・ 305.0
40 パーセンタイル・・・ 998.0
50 パーセンタイル・・・ 1263.0
60 パーセンタイル・・・ 1306.0
70 パーセンタイル・・・ 1323.0
80 パーセンタイル・・・ 1397.0
90 パーセンタイル・・・ 1503.0
100 パーセンタイル・・・ 1986.0
'''

ここまではわかりやすいのですが、自分が疑問に思ったのは、
もっと中途半端なパーセンタイルです。

(例)この出力の40.16ってどうやって算出された?


data_4 = np.array([15, 52, 100, 73, 102])
print(np.percentile(data_4, 17))
#  出力
# 40.16

この疑問放置したままなのが気持ち悪かったので、
これまでパーセンタイルや四分位数、そしてこれらを使う箱ひげ図などを使わなかったのですが、
とあるタスクの中で箱ひげ図を使いたくなったのでこの機会に仕様を確認しました。

といっても、numpyの該当ページにもNote.として記されていますし、
wikipediaにも普通に載ってます。
分位数
あと、pを1刻みで動かして適当なデータに対してパーセンタイル算出してプロットしたら明快にわかりました。

要は、中途半端な値に対しては、隣接の2つの値を線形補完して求めるそうです。
上の例で言えば、
15が0パーセンタイル、52が25パーセンタイルなので、17パーセンタイルは
$(52-15)*17/25+15=40.16$ と計算されています。
仕様がわかったのでこれからはバシバシ使おう。

機械学習を行う時、異常時の前処理として、1〜99パーセンタイルでクリップすると有効なことがあるという話を最近聞いたので、
それも試してみたいです。