月: 2023年10月

いつも使っているJupyterに関する話です。Jupyterの中では、InteractiveShellというオブジェクトが使われていて、これを利用することでさまざまな高度な操作を行うことができます。

多くのことができるのですが、最近これを使って変数を操作する機会があったのでその辺を紹介していきます。

そもそも、InteractiveShell オブジェクトとは？という話ですが、これはIPython（Jupyter Labのバックエンドとして動作する対話的なPythonシェル）の中心的なクラスです。このオブジェクトを通じて、現在のセッションの状態を取得したり、マジックコマンドを実行したりすることができます。

このInteractiveShellはget_ipython() という関数を使うことでそのnotebook自身で取得し、アクセスすることができます。

マジックコマンドの動的な実行

このInteractiveShellオブジェクトを使うと例えば、マジックコマンドをプログラムで実行したりできます。

例えば、%pwd っていう現在のでディレクトリを取得するマジックコマンドがありますが、これが次のようにして実行できます。

ip = get_ipython()
ip.run_line_magic('pwd', '')

(pwdをメソッドとして実行するメリットはなかなか薄いので例が適切でなくてすみません。)

マジックコマンドをそのままマジックコマンドととして実行するのと比べると、結果を文字列として取得して以降の処理で使えるとか、実行するコマンドをプログラムで動的に変化させて実行させることができるといった利点があります。

もちろん、 run_cell_magic() って関数もあるのでセルマジックも実行できますよ。

変数の取得や設定

今回紹介するもう一つの使い方は変数の取得や追加です。名前空間の操作といった方がわかりやすいかもしれません。

InteractiveShellオブジェクトを使って、upyter Labのセッション内で定義されている変数を取得したり、新しい変数を設定したりすることができます。

変数の取得は user_ns 属性を使います。現在定義されている変数が全部辞書として入っています。デフォルトで入ってるJupyter/IPython固有の変数も多いので、自分が宣言した覚えのないものもあると思います。

例えば次のようにして使います。

x = 10
y = "Hello"

ip = get_ipython()
print(ip.user_ns['x'])
# 10
print(ip.user_ns['y'])
# "Hello"

変数の取得としてはあまり用途ないかもしれないですね。どちらかというと既に宣言されているかどうかのチェックとかで使えるかもしれません。

逆に、次のようにして新しい変数を定義することもできます。

ip = get_ipython()
ip.push({'z': "新しい変数"})

print(z)
# 新しい変数

この例だけ見てると、 z=”新しい変数” って直接書くのと同じやん、と思えますが、この方法の利点は変数名の方もプログラムで動的に指定できることです。

まぁ、 exec(“z=’新しい変数'”) ってすれば変数名を動的に変えられるので固有のメリットというわけではないのですが、execは本当にリスクの大きな手段なのでそれよりは安全に使えるという利点があると思います。

この変数の名の指定は、関数やマジックコマンドの中で使うと、例えばユーザーが指定した変数名に結果を格納するといった応用ができます。

最近、別の場所で書いているnote記事の中で自分が自作したマジックコマンドを紹介したのですが、そこで使っているクエリの結果を引数で指定した変数に格納するのにこのip.pushを使っています。

参考: Snowflakeに手軽にSQLを打てるJupyterマジックコマンドを作ってみた

ip.push を使わずに、
ip.user_ns[‘z’] = “新しい変数”
みたいに、user_nsオブジェクトに直接突っ込んでも同じことができますので合わせて覚えておいてください。

InteractiveShellオブジェクトでできることのごく一部ではありますが、今回の記事としては以上となります。

はじめに

ちょっと重い記事が続いてたので今回は小ネタです。

二つの区間データ、要するに数値の区間であれば最小値と最大値のペア、時刻の区間であれば開始時間と終了時間のデータが二つあったときに、それらの区間に重複があるかどうかをサクッと判定するアルゴリズムを紹介します。昔どこかでみた覚えがあるのですが、久々に実装で必要になったときにちょっとだけ悩んだのでそのメモです。

結論

結果だけ知りたい人向けに先に結論だけ書いときます。

区間[s1, e1]と[s2, e2]に重複があるかどうかは次の2つの方法で判定できます。

方法1: $(s1 \leq e2) \land (s2 \leq e1)$ ならばその2区間には重複があります。

方法2: $\max(s1, s2) \leq \max(e1, e2)$ ならばその2区間には重複があります。

両方の方法のメリットデメリットですが、方法1の方が計算速度が速いです。大小比較2個と論理演算ですからmax関数呼び出すより有利ですね。ただ、方法2の方は、応用として、二つの区間に重複があった場合、$\max(e1, e2) – \max(s1, s2)$で重複区間の長さを算出できます。(値が負になったら重複無し。)

(max関数がそんな極端に遅いってこともなく、わざわざ時間測定しない限り体感することもない速度差ではあるのであくまでも参考にどうぞ。)

方法1の導出

これだけで終えてしまうと記事量としてあまりに少ないのと、結論だけ書いてるとすぐ忘れそうなので、真面目に導出を説明しておきます。まずは方法1の方からです。

まず、二つの区間の相対的な位置関係は次の6パターンが存在します。(説明の単純化のため。s1,s2,e1,e2は全て異なる値とします。)

1. 区間1全体が区間2より小さい。つまり、 $s1 < e1 < s2 < e2$.
2. 区間1全体が区間2より大きい。つまり、 $s2 < e2 < s1 < e1$.
3. 区間1が区間2に内包される。つまり、$s2 < s1 < e1 < e2$.
4. 区間1が区間2を内包する。つまり、$s1 < s2 < e2 < e1$.
5. 区間1の後半と区間2の前半が重なる。つまり、$s1 < s2 < e1 < e2$.
6. 区間1の前半と区間2の後半が重なる。つまり、$s2 < s1 < e2 < e1$.

上記の6パターンのうち、1,2 が区間に重複がなく、3,4,5,6は区間に重複があります。

ここで、「そうか！3,4,5,6のどれかを満たすことを確認すればいいんだ！」と判断して4パターンの論理式を書いてorで繋ぐ、としないことがコツです。

区間に重複がないパターンの方が2種類しかなくて判定が単純なんですよね。そして、区間の定義から$s1<e1$と$s2<e2$はわかってるので、あと確認するのはe1とs2, e2とs1の代償関係だけなんですよ。

ということで、 $(e1 < s2) \lor (e2 < s1)$ であれば二つの区間に重複はないと言えます。

これの否定を考えると、ドモルガンの法則から方法1としてあげた式が導かれます。

方法2の導出

続いて方法2の導出です。

これは、先ほど挙げた3,4,5,6の4パターンの不等式を眺めると見つかる法則性なのですが、左の2辺はs1,かs2で、右の2辺はe1, e2なんですよね。

これを素直に数式に落とすと、$\max(s1, s2) \leq \max(e1, e2)$ となります。等号が成立するのはただ1点だけが重複する場合。

そして、区間が重複する場合は、$\max(s1, s2)$は重複区間の開始点であり、$\max(e1, e2)$は重複区間の終了点なので、この2項の差を取ると重複部分の長さも得られます。

まとめ

そんなに難しい計算ではなく、覚えてさえればサクッと実装できる式ではありますが、重複のパターンって4種類あるよなぁと考え始めてしまうと意外に手間取ります。

結果を丸暗記しなくても、区間が重複しないパターンは2個しかなくてそれを否定したら簡単だってことを頭の片隅にでも置いといてもらえると幸いです。

エントロピーや相互情報量の記事が続いていますが、今回の記事で計算の実装方法を紹介して一旦区切りとします。

エントロピーも相互情報量も数式はそこまで難しくないので、numpyで定義通り計算しても良いのですが、エントロピー関係はSciPyに、相互情報量はscikit-learnに用意されているので今回の記事ではそれを使っていきます。

計算対象のデータは、[“a1”, “a2”, “a1”, “a1”, “a2”] みたいにローデータの一覧で保有している場合もあれば、”a1″が3個で”a2″が2個のようにカウント済みのものがある場合もあると思うのでそれぞれ説明していきます。

エントロピーの計算

まず一番基本的なエントロピーの計算からです。これは、scipy.stats.entropy メソッドを使います。
参考: scipy.stats.entropy — SciPy v1.11.3 Manual

基本的な引数はpkなので、確率の一覧を渡すのが想定されていますが、和が1でないなら1になるように正規化してくれるのでサンプルがある場合は個数を渡しても大丈夫です。また、base引数で対数関数の底を指定でき、デフォルトが$e$なので、情報理論で使う場合は$2$を指定しましょう。

やってみます。

import numpy as np  # データ作りに利用
import pandas as pd  # データ作りに利用
from scipy.stats import entropy


pk = np.array([1/2, 1/3, 1/6])  # 確率の一覧が得られた場合。
print(entropy(pk, base=2))
# 1.459147917027245

count_list = np.array([3, 2, 1])  # データの個数の場合
print(entropy(count_list, base=2))
# 1.4591479170272446

# カウント前のデータの一覧がある場合
data_sr = pd.Series(["a1", "a1", "a1", "a2", "a2", "a3"])
# value_counts()で数えあげたものをentropyに渡す
print(entropy(data_sr.value_counts(), base=2))
# 1.4591479170272446

結合エントロピーの計算

次は結合エントロピーです。エントロピーを単純に2次元に拡張したやつですね。(条件付きエントロピーではないので注意してください、

例えば次のような例を考えましょうか。

結合エントロピーの場合はですね、元のカウントデータを2次元から1次元に並び替えて渡します。

matrix_data = np.array([[4, 1], [2, 3]])
print(matrix_data)
"""
[[4 1]
 [2 3]]
"""

# ravel か flattenで1次元化して計算する
print(entropy(matrix_data.ravel(), base=2))
# 1.8464393446710157

# 標本データがある場合
df = pd.DataFrame({
    "A": ["a1", "a1", "a1", "a1", "a1", "a2", "a2", "a2", "a2", "a2"],
    "B": ["b1", "b1", "b1", "b1", "b2", "b1", "b1", "b2", "b2", "b2"],
})

# カウントしたデータを使う
print(df.groupby(["A", "B"]).size())
"""
A   B
a1  b1    4
    b2    1
a2  b1    2
    b2    3
dtype: int64
"""

print(entropy(df.groupby(['A', 'B']).size(), base=2))
# 1.8464393446710157

条件付きエントロピー

次は条件付きエントロピーです。残念なことなのですが、メジャーなライブラリでは条件付きエントロピー専用の関数は提供されていません。

そこで、$H(A|B) = H(A, B) – H(B)$などのエントロピー間の関係式を使って計算することになります。相互情報量も含めて、$H(A|B) = H(A) – I(A; B)$などで計算してもいいのですが、SciPyで完結できるので最初の式のほうが良いでしょう。

先ほどの表データをサンプルとします。$H(B)$については、表データを縦に足し合わせてBだけのカウントデータを作って計算します。

data_B = matrix_data.sum(axis=0)
print(data_B)
# [6 4]

# H(B)の計算
entropy_B = entropy(data_B, base=2)
print(entropy_B)
# 0.9709505944546688

# H(A, B)の計算
joint_entropy = entropy(matrix_data.ravel(), base=2)
print(joint_entropy)
# 1.8464393446710157

# H(A|B) = H(A, B) - H(B)
conditional_entropy_A_given_B = joint_entropy - entropy_B
print(conditional_entropy_A_given_B)
# 0.8754887502163469

# 標本データがある場合
entropy_B = entropy(df["B"].value_counts(), base=2)
joint_entropy = entropy(df.groupby(["A", "B"]).size(), base=2)

conditional_entropy_A_given_B = joint_entropy - entropy_B
print(conditional_entropy_A_given_B)
# 0.8754887502163469

以上で、3種類のエントロピーが計算できました。

相互情報量

最後に相互情報量の計算方法です。

$I(A; B) =H(A)-H(A|B)$など複数の表現方法があるので、ここまでに計算してきた値から算出することもできます。

entropy_A = entropy(df["A"].value_counts(), base=2)
print(entropy(df["A"].value_counts(), base=2) - conditional_entropy_A_given_B)
# 0.12451124978365313

ただ、scikit-learnに専用のメソッドがあるのでこちらの使い方も見ておきましょう。
参考: sklearn.metrics.mutual_info_score — scikit-learn 0.18.2 documentation

引数は、mutual_info_score(labels_true, labels_pred, contingency=None)
となっており、標本データを受け取るのが標準的な使い方で、その第一,第二引数をNoneにしてcontingency引数にカウントデータを渡すこともできます。(contingencyがNoneでない場合はこれが優先されて、先の二つの引数が無視されます。)

1点注意しないといけないのは、entropyと違って対数の底が指定できず、自然対数に固定されてしまうことです。底を$2$で考えたい場合は、$\ln{x}/\ln{2} = \log_2{x}$を使って変換が必要です。

from sklearn.metrics import mutual_info_score


# np.log(2)で割ることを忘れない
# カウントした表データがある場合
print(mutual_info_score(None, None, contingency=matrix_data)/np.log(2))
# 0.12451124978365345

# 標本データがある場合
print(mutual_info_score(df["A"], df["B"])/np.log(2))
# 0.12451124978365345

計算の都合上超軽微な誤差がありますが、それ以外は想定通りの値が得られていますね。

以上で、相互情報量も計算できるようになりました。