ハミング距離(Hamming distance)

二つの文字列がどのくらい異なるかを表す距離関数として、以前レーベンシュタイン距離と言うのを紹介しました。
参考:pythonで編集距離(レーベンシュタイン距離)を求める

これよりももっとシンプルな関数として、ハミング距離(Hamming distance)と言うのがあるのでその紹介です。
これは文字数が同じ二つの文字列に対して、単純に異なる文字を数えたものです。
ハミング距離 – Wikipedia

自分の場合、文字数が同じ時しか定義できないなどの理由により、レーベンシュタイン距離に比べて使う頻度は低いです。
ただ、文字数が同じでさえあれば高速に計算できます。

pythonでの実装ですが、レーベンシュタイン距離の時に使った、
python-Levenshtein に含まれています。

試しにやってみましょう。


import Levenshtein
print(Levenshtein.hamming("toned", "roses"))
# 3

この例の toned と roses では、ハミング距離もレーベンシュタイン距離もどちらも3ですが、
文字数が同じであってもこの2種類の距離の値が異なる例はあります。

例えば次のようなケースです。


text1 = 'abcdefg'
text2 = 'bcdefga'
print(Levenshtein.distance(text1, text2))
# 2
print(Levenshtein.hamming(text1, text2))
# 7

Macでlsコマンドの出力結果の文字色を変更する

今回もMacの小ネタ。
Macのターミナルでlsコマンドを使うと、ディレクトリに色を塗ってくれますが、背景色との組み合わせによっては見辛いことがあります。
(具体的に言えば、黒背景で使っている時の通常ディレクトリの青は見辛いです)

これらの色は、 LSCOLORS という環境変数に適切に値を設定することで変更可能です。
設定内容は、 ls コマンドのマニュアル内で検索すると出てきます。

$ man ls
で、マニュアルを開いて、
/LSCOLORS
で検索しましょう。

自分のMacで実行した時の該当部分を載せておきます。

LSCOLORS The value of this variable describes what color to use for which attribute when col-
ors are enabled with CLICOLOR. This string is a concatenation of pairs of the format
fb, where f is the foreground color and b is the background color.

The color designators are as follows:

a black
b red
c green
d brown
e blue
f magenta
g cyan
h light grey
A bold black, usually shows up as dark grey
B bold red
C bold green
D bold brown, usually shows up as yellow
E bold blue
F bold magenta
G bold cyan
H bold light grey; looks like bright white
x default foreground or background

Note that the above are standard ANSI colors. The actual display may differ depend-
ing on the color capabilities of the terminal in use.

The order of the attributes are as follows:

1. directory
2. symbolic link
3. socket
4. pipe
5. executable
6. block special
7. character special
8. executable with setuid bit set
9. executable with setgid bit set
10. directory writable to others, with sticky bit
11. directory writable to others, without sticky bit

The default is “exfxcxdxbxegedabagacad”, i.e. blue foreground and default background
for regular directories, black foreground and red background for setuid executables,
etc.

初期設定はこれです。
exfxcxdxbxegedabagacad

この22文字が、2文字ずつ”文字色””背景色”のペアになっていて、マニュアル中の11種類に対応しています。
最初の ex が directory の配色で、 blue の文字と、デフォルトの背景色です。

これが読みにくいので、 g : cyan あたりに変更するには、
gxfxcxdxbxegedabagacad と設定すれば大丈夫です。

.bash_profile に以下のように設定しておきましょう。


export LSCOLORS=gxfxcxdxbxegedabagacad

濃い青が明るい水色になって(僕の環境では)とても読みやすくなります。

Note that the above are standard ANSI colors. The actual display may differ depending on the color capabilities of the terminal in use.
と注意されている通り、環境によって色は違うので設定を変更するときは試しながら行うことをお勧めします。
例えば今の僕の環境では d : brown は黄色いです。

Mac の起動時刻を調べる

先日分け合って、使っているMacを何時に起動したのか調べる必要が発生しました。

結論から言うと、次のコマンドで調べることができます。
rebootなのに再起動だけでなく通常の起動もわかります。


$ last reboot
# --- 略 ---
reboot    ~                         Tue Jan  1 17:13
reboot    ~                         Thu Nov 29 19:28
reboot    ~                         Sun Nov 25 12:47
reboot    ~                         Thu Nov 22 14:51
reboot    ~                         Wed Nov 21 20:18

また、シャットダウン時刻を調べたい時は以下のコマンドです。 (出力略)


$ last shutdown

last とだけ打つと両方同時に表示されます。

もっと細かい使い方はそのうち調べようと想うのですが、
取り急ぎ今回の要件は満たしたのでメモしておきました。

pythonのfilter関数の使い方

前回の記事がmap関数の話だったので、今回は使い方のよく似たfilter関数です。

ドキュメントはこちら。
組み込み関数

基本的な構文は以下の形で、iterable の各要素に functionを適用して、
結果が新なものだけを取り出せます。
filter(function, iterable)

map関数と同様に、戻り値はリストではなくイテレーターになるので最初は少し戸惑いました。
試しに、整数のうち偶数だけ抽出するフィルターを書いてみます。


def f(n):
    return n % 2 == 0


m = filter(f, range(10))
print(m)
# <filter object at 0x114d178d0>
print(list(m))
# [0, 2, 4, 6, 8]
print(list(m))
# []

細かい説明は mapの記事に書いた通りなのですが、
イテレーターを使うために、リストに変換するなり、nextで取り出すなりする必要があり、
一度取り出すともう一回listに変換しようとしても空のリストしか返ってこなくなります。

なお、内包表記でもほぼ同じ処理が実装でき、自分はこちらを使うことが多いです。


[x for x in range(10) if f(x)]
[0, 2, 4, 6, 8]

pythonのmap関数の使い方

前の記事でさらっと使っていたmap関数の使い方の紹介です。

ドキュメントはこちら。
組み込み関数

MathematicaのMapとよく似た関数(と聞いて、「ああ、あれね」と通じる人がどのくらいいらっしゃるのかわかりませんが)であり、
配列などの各要素に関数を順番に適用することができます。

使い方は map(適用したい関数, 配列1, ) です。
僕はpython初心者の頃、関数を適用した結果のリストがすぐ戻ってくると期待していたのに、イテレーターが戻ってきたので結構戸惑いました。

例えば、引数を2乗する関数で試してみましょう。


def f(x):
    return x**2


m = map(f, range(10))
print(m)
# <map object at 0x1156f19b0>

リストにしたければ、型変換してあげる必要があります。
ただし、通常のリストと違い、イテレーターなので、一度最後のデータかまで取り出すと、次の値が取れなくなります。


# 1回目は与えられたリストに関数を適用した結果が戻る
print(list(m))
# [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

# 続けて全く同じように呼び出すと空のリストがかえる。
print(list(m))
# []

内包表記でほぼ同じ処理を実装できますが、違いは関数が実行されるタイミングです。
内包表記は、それが定義されたタイミングで計算され、
mapの場合は、値が必要になったタイミングで実行されます。

適用する関数にprint文を仕込んで、様子を見てみましょう。
最初に内包表記の場合、


def g(x):
    print(x)
    return x**2


l = [g(x) for x in range(10)]
# この時点で g が実行されているので、以下が出力される。
'''
0
1
2
3
4
5
6
7
8
9
'''
print(l)
# [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

次にmapの場合。


m1 =  map(g, range(10))
# まだgが実行されてないので、何も出力されない。

# nextを使って、最初の3個の値を取り出すと、その3この値に対してだけ関数gが実行される。
print(next(m1))
'''
0
0
'''
print(next(m1))
'''
1
1
'''
print(next(m1))
'''
2
4
'''

これらの性質により、上手く使えば実行時間やメモリの節約が可能になるそうです。
(それを実感するほど上手く使えたことはほとんどないのですが)
ただ、pythonを使っていく上で、イテレーターの扱いに慣れておくのは有益なので、学んでおいて損はなかったと思ってます。

pythonのthisモジュールに定義されている変数について

以前紹介した The Zen of Python の記事の中で、
thisというモジュールに仕掛けられたEaster Eggを紹介しました。

その記事中ではインポートした瞬間に表示される文字列にしか焦点を当ててませんでしたが、
このモジュールに含まれている関数や定数についてもちょっと面白いので紹介します。

このthisの中で、どんな変数やメソッドが定義されているのか、dir関数で見てみましょう。


import this
# 略

dir(this)
'''
['__builtins__',
 '__cached__',
 '__doc__',
 '__file__',
 '__loader__',
 '__name__',
 '__package__',
 '__spec__',
 'c',
 'd',
 'i',
 's']
'''

“__”で始まる特別な値以外に、c, d, i, s の4つの変数が含まれています。
これらのうち、iとcはそれぞれ整数ですがあまり意味はありません。
ちょっと面白いのは、dとsです。

まず、 dの方は次の辞書です。


print(this.d)
# 以下出力
{'A': 'N', 'B': 'O', 'C': 'P', 'D': 'Q', 'E': 'R', 'F': 'S', 'G': 'T', 'H': 'U', 'I': 'V', 'J': 'W', 'K': 'X', 'L': 'Y', 'M': 'Z', 'N': 'A', 'O': 'B', 'P': 'C', 'Q': 'D', 'R': 'E', 'S': 'F', 'T': 'G', 'U': 'H', 'V': 'I', 'W': 'J', 'X': 'K', 'Y': 'L', 'Z': 'M', 'a': 'n', 'b': 'o', 'c': 'p', 'd': 'q', 'e': 'r', 'f': 's', 'g': 't', 'h': 'u', 'i': 'v', 'j': 'w', 'k': 'x', 'l': 'y', 'm': 'z', 'n': 'a', 'o': 'b', 'p': 'c', 'q': 'd', 'r': 'e', 's': 'f', 't': 'g', 'u': 'h', 'v': 'i', 'w': 'j', 'x': 'k', 'y': 'l', 'z': 'm'}

そして、 sは次の意味不明な文字列が入ってます。


print(this.s)
Gur Mra bs Clguba, ol Gvz Crgref

Ornhgvshy vf orggre guna htyl.
Rkcyvpvg vf orggre guna vzcyvpvg.
Fvzcyr vf orggre guna pbzcyrk.
Pbzcyrk vf orggre guna pbzcyvpngrq.
Syng vf orggre guna arfgrq.
Fcnefr vf orggre guna qrafr.
Ernqnovyvgl pbhagf.
Fcrpvny pnfrf nera'g fcrpvny rabhtu gb oernx gur ehyrf.
Nygubhtu cenpgvpnyvgl orngf chevgl.
Reebef fubhyq arire cnff fvyragyl.
Hayrff rkcyvpvgyl fvyraprq.
Va gur snpr bs nzovthvgl, ershfr gur grzcgngvba gb thrff.
Gurer fubhyq or bar-- naq cersrenoyl bayl bar --boivbhf jnl gb qb vg.
Nygubhtu gung jnl znl abg or boivbhf ng svefg hayrff lbh'er Qhgpu.
Abj vf orggre guna arire.
Nygubhtu arire vf bsgra orggre guna *evtug* abj.
Vs gur vzcyrzragngvba vf uneq gb rkcynva, vg'f n onq vqrn.
Vs gur vzcyrzragngvba vf rnfl gb rkcynva, vg znl or n tbbq vqrn.
Anzrfcnprf ner bar ubaxvat terng vqrn -- yrg'f qb zber bs gubfr!

ただこの文字列、薄目にみてみると、各単語の文字数が The Zen of Python が似ています。

実はこれはちょっとした暗号になっていて、
this.s の各文字を、this.dの辞書で変換すると、 The Zen of Python が現れます。

書き方はいろいろあると思いますが次のような形でやってみましょう。


print("".join(map(lambda x: this.d.get(x, x), this.s)))
# 以下出力
The Zen of Python, by Tim Peters

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!

モジュール this の実装自体が
The Zen of Python に反してややこしい値になっているというちょっとした遊び心のようです。

ちなみに、 this モジュールのソースコードも読んでみたのですが、 mapではなく、内包表記でやっているようですね。
sやdを定義したあとに、次のように書いてありました。
よく考えなくてもこれで十分ですね。


print("".join([d.get(c, c) for c in s]))

DataFrameのsampleメソッドのドキュメントを読む

超高頻度で使っているメソッドなのに、公式ドキュメントを読んだことがなかった、
pandas.DataFrame.sample
についてドキュメンを読んでみました。

元々の目的はデータ件数が不明なデータフレームからn個のサンプルが欲しい時に、
df.sample(n)とすると、データフレームの件数が少ないとエラーになるのが面倒だし、
事前にlen(df)して、条件分岐するのが面倒なので都合のいいオプションを探していました。


df.sample(200)
# ValueError: Cannot take a larger sample than population when 'replace=False'

これを回避するために、いつも何かしら工夫をしていますが、正直無駄に行数が増えてる気がしています。
(コードのイメージ)


n = 200
if len(df) >= n:
    df.sample(n)
else:
    df

dfの行数がnより大きいならn件返して欲しくて、行数が少ないなら全部のデータをそのまま渡すというのを、
if文を使わずにsample()の引数で実現したかったわけです。
(結論からいうと、そのような引数は用意されていませんでした。)

エラーメッセージの指示に従ってreplace=True を設定すると、
同じ行が複数回サンプリングされるのを許すのでnが行数より大きくても大丈夫になります。
(ただ、これは自分が元々やりたかったのとは違う。)

dfを150行のデータフレームとすると次のような感じ。


print(len(df))  # 150
print(len(df.sample(200, replace=True)))  # 200

目的のオプションは見つかりませんでしたが、その代わり、今まで知らなかった引数が使えることを知りました。
正確なステートメントは公式ドキュメントの方に任せて、ざっくりと書くと以下のようなものが使えます。

frac : サンプリングするデータの数を、個数ではなく割合で指定する。
replace : Trueにすると、同じデータを重複してサンプリングできる。
weights : サンプリングされる確率を重み付けできます。
random_state : 乱数固定。
axis : 1を指定すると、行ではなく列をサンプリングできる。

決定木系の機械学習モデルを自分で実装するときなどに便利そうですね。

matplotlibのグラフ間の間隔の調整

普段の分析でもこのブログでも、matplotlibの一つのfigureに複数のグラフを描写することがよくあります。

その時、稀に困るのがタイトルやラベルが他のグラフに重なってしまう時です。
それ以外にも微妙に感覚が詰まりすぎてたり広すぎたりで見た目が悪いなぁと思うことがあります。

その場合、subplots_adjustを使って、間隔を調整できます。
ドキュメントはこちら。
matplotlib.pyplot.subplots_adjust
引数の意味はこちらにもあります。
matplotlib.figure.SubplotParams

次の6種類の引数を渡せますが、間隔の調整に使うのは、hspaceとwspaceです。

The parameter meanings (and suggested defaults) are:

left = 0.125 # the left side of the subplots of the figure
right = 0.9 # the right side of the subplots of the figure
bottom = 0.1 # the bottom of the subplots of the figure
top = 0.9 # the top of the subplots of the figure
wspace = 0.2 # the amount of width reserved for space between subplots,
# expressed as a fraction of the average axis width
hspace = 0.2 # the amount of height reserved for space between subplots,
# expressed as a fraction of the average axis height
Copy to clipboard
The actual defaults are controlled by the rc file

無理やり作った例で恐縮ですが、調整せずにグラフ間で重なってしまった例と、
subplots_adjustで調整した例を作りました。


import matplotlib.pyplot as plt
x = [1, 2, 3, 4, 5]
y = [1, 10, 100, 1000, 10000]
# 悪い例
fig = plt.figure(figsize=(9, 6),  facecolor="white")
for i in range(3):
    for j in range(3):
        ax = fig.add_subplot(3, 3, 1 + j + 3*i)
        ax.set_title("複数行の\nタイトル")
        ax.plot(x, y)
plt.show()
# 間隔を調整した例
fig = plt.figure(figsize=(9, 6),  facecolor="white")
fig.subplots_adjust(hspace=0.6, wspace=0.4)
for i in range(3):
    for j in range(3):
        ax = fig.add_subplot(3, 3, 1 + j + 3*i)
        ax.set_title("複数行の\nタイトル")
        ax.plot(x, y)
plt.show()

出力はそれぞれこちらです。

hspace=0.6 と wspace=0.4 の値ですが、これはドキュメント読んで何かしら計算して決めるより、
適当な値を入れて何度か試すのがおすすめです。
初期値はそれぞれ0.2なので、それより大きい値を入れると広くなります。

単位根過程の見せかけの回帰

久しぶりに時系列解析の話題です。

以前の記事で、単位根過程という用語を紹介しました。
参考:差分系列と単位根過程
それに関して、見せかけの回帰(spurious regression)という現象があるのでそれを紹介します。

いつもの沖本先生の本から引用します。(127ページ)

定義 (見せかけの回帰)
単位根過程$y_t$を定数と$y_t$と関係のない単位根過程$x_t$に回帰すると、$x_t$と$y_t$の間に有意な関係があり、
回帰の説明力が高いように見える現象は見せかけの回帰と言われる。

この現象を確認するために、実際に単位根過程を2個生成して回帰をやってみましょう。
単位根過程は単純に次の二つのランダムウォークでやります。
それぞれ独立した系列$\varepsilon_{1,t}$と$\varepsilon_{2,t}$から生成しているので、これらは無関係な過程です。
$$
x_t = x_{t-1} + \varepsilon_{1,t} \ \ \ \varepsilon_{1,t}\sim N(0,1)\\
y_t = y_{t-1} + \varepsilon_{2,t} \ \ \ \varepsilon_{2,t}\sim N(0,1)
$$

回帰係数の有意性もみたいので、statsmodelsで実装しました。


import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
# データ生成
T = 200
X = np.cumsum(np.random.randn(T)).reshape(-1, 1)
y = np.cumsum(np.random.randn(T))


# 定数項を作成
X0 = sm.add_constant(X)
# 回帰分析実行
st_model = sm.OLS(y, X0)
results = st_model.fit()
# 結果表示
print(results.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.171
Model:                            OLS   Adj. R-squared:                  0.167
Method:                 Least Squares   F-statistic:                     40.85
Date:                Sat, 25 May 2019   Prob (F-statistic):           1.15e-09
Time:                        23:16:39   Log-Likelihood:                -423.18
No. Observations:                 200   AIC:                             850.4
Df Residuals:                     198   BIC:                             856.9
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -5.2125      0.327    -15.960      0.000      -5.856      -4.568
x1            -0.1826      0.029     -6.391      0.000      -0.239      -0.126
==============================================================================
Omnibus:                        5.521   Durbin-Watson:                   0.215
Prob(Omnibus):                  0.063   Jarque-Bera (JB):                3.122
Skew:                          -0.033   Prob(JB):                        0.210
Kurtosis:                       2.391   Cond. No.                         26.3
==============================================================================
'''

結果を見ると、係数のp値が非常に小さくなっていて、$x_t$と$y_t$の間に関係があるように見えます。

定常過程同士であればこのような現象は発生しないことを見るために、
$x_t$と$y_t$の差分系列をとって、同じように回帰してみましょう。
(差分系列はどちらも定常過程です。)


X_diff = X[1:] - X[:-1]
y_diff = y[1:] - y[:-1]
# 定数項を作成
X0_diff = sm.add_constant(X_diff)
# 回帰分析実行
st_model = sm.OLS(y_diff, X0_diff)
results = st_model.fit()
# 結果表示
print(results.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                 -0.005
Method:                 Least Squares   F-statistic:                   0.03113
Date:                Sat, 25 May 2019   Prob (F-statistic):              0.860
Time:                        23:48:11   Log-Likelihood:                -265.29
No. Observations:                 199   AIC:                             534.6
Df Residuals:                     197   BIC:                             541.2
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0195      0.065     -0.298      0.766      -0.148       0.109
x1            -0.0117      0.066     -0.176      0.860      -0.142       0.119
==============================================================================
Omnibus:                        2.378   Durbin-Watson:                   2.107
Prob(Omnibus):                  0.305   Jarque-Bera (JB):                2.338
Skew:                           0.263   Prob(JB):                        0.311
Kurtosis:                       2.923   Cond. No.                         1.01
==============================================================================
'''

p値が明らかにでかい値になりましたね。

今回の記事ではグラフは省略しましたが、
単純にプロットしてみてみると確かに単位根過程同士では何らかの関連を持っているように見え、
定常過程同士では無関係に見えます。

スタージェスの公式によるヒストグラムのビンの数の決定

データの傾向を見るときにヒストグラムを描くことは頻繁にありますが、
そのとき課題になるのが、ビンの数を何本にするかです。

いつもmatplotlibのデフォルトである10本でとりあえず試したり、
適当に変えながら何パターンか試したりしています。
Tableauであれば、本数ではなく幅での指定ですね。

ただ、何かしら参考指標が欲しいとは思っていたので調べたところ、複数の方法が提案されていました、
その中でスタージェスの公式(Sturges’ formula)が良さそうだったので試してみました。

定義はwikipediaのものを採用しましょう。
(他のサイトを見ると、小数点以下の扱いで微妙に異なるパターンがあります。)
ヒストグラム – Wikipedia

スタージェスの公式によると、n個のデータがあるとき、ビンの数kの目安は次の式で得られます。

$$k = \lceil \log_{2}{n} + 1\rceil$$

注意として、スタージェスの公式はその導出の背景に、二項分布が正規分布で近似できるという性質を使っています。
そのため、nが小さすぎる場合にはあまり参考になりません。
(そのような時はそもそも、ヒストグラム自体があまり有効ではないです。)
また、 データの分布が二項分布/正規分布と大きく異なる時もうまくいきません。

それではいくつかのデータで試してみましょう。
次のコードは、7種類の件数に対して、ランダムにデータを取得し、
スタージェスの公式で得られたビンの数ののヒストグラム(中央列)と、
それよりビンが2本少ないヒストグラム(左列)、2本多いヒストグラム(右列)を描写します。


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt


def sturges_formula(data_size):
    return int(np.floor(1+np.log2(data_size)))


fig = plt.figure(figsize=(15, 28), facecolor="white")
for i in range(7):
    data_size = int(1.5 * 2 ** (i+4))
    data = np.random.randn(data_size)
    bin_count = sturges_formula(data_size)

    for j in range(3):
        ax = fig.add_subplot(7, 3, 3*i+j+1)
        ax.set_title("データ件数:{d}件, ビン:{b}本".format(
            d=data_size, b=bin_count-2+2*j
        ))
        ax.hist(data, bins=bin_count-2+2*j, rwidth=0.8)

plt.show()

乱数を使っているので毎回結果は変わりますが、出力の一例がこちら。

左列は確かに若干少ないかなという気がします。
ただ、データ件数が多くなると、もう少し多い数(右列)でもいいので、やはり目安として使うのが良さそうですね。