分析ノート

個人的な好みの話ですが、僕は割合の可視化は円グラフより100%積み上げ棒グラフの方が好みです。前職時代はこの種の可視化はもっぱらTableauで描いていたのですが転職して使えなくなってしまったので、Pythonでサクッと描ける方法をメモしておきます。

最初、matplotlibのbarやbarhでbottomやleftを逐一指定してカテゴリごとに描いていく方法を紹介しようと思っていたのですが、pandas.DataFrameのメソッドを使った方が簡単だったのでそちらを先に紹介します。

データの準備

サンプルデータが必要なので適当にDataFrameを作ります。地域ごとに、どの製品が売れているか、みたいなイメージのデータです。

import pandas as pd


data = {
    '製品A': [20, 30, 50],
    '製品B': [60, 70, 40],
    '製品C': [20, 0, 10]
}
df = pd.DataFrame(data, index=["地域1", "地域2", "地域3"])

print(df)
"""
     製品A  製品B  製品C
地域1   20   60   20
地域2   30   70    0
地域3   50   40   10
"""

こちらを地域別に、売れている製品の内訳を可視化していきます。
これも個人的な好みの話ですが、棒グラフは縦向きではなく横向きにします。(ラベルが日本語だとその方が自然な結果になりやすいからです。)

シンプルにやるには、jupyter環境であれば以下のコードだけで目的のグラフが表示されます。

# データの正規化
df_normalized = df.div(df.sum(axis=1), axis=0)
# 横向きの積み上げ棒グラフの描画
df_normalized.plot(kind='barh', stacked=True)

少し脇道に逸れますが、みなさん、この行ごとに正規化するdf.div(df.sum(axis=1), axis=0) ってやり方ご存知でした？昔、「Pandasのデータを割合に変換する」って記事を書いたときは知らなかった方法でした。このときは転置して列ごとに正規化してもう一回転置するという手順を踏んでましたね。

記事の本題に戻ると、plotメソッドのstackedっていう引数をTrueにしておくと棒グラフを積み上げて表示してくれるため目的のグラフになります。
参考: pandas.DataFrame.plot.barh — pandas 2.1.3 documentation

ただ、これで出力すると、グラフの棒の並び順が上から順に地域3, 地域2, 地域1 となり、ちょっと嫌なのと、あとラベル等のカスタマイズも行えた方が良いと思うのでもう少し丁寧なコードも紹介しておきます。(こだわりなければさっきの2行で十分です。)

import matplotlib.pyplot as plt


# データの正規化
df_normalized = df.div(df.sum(axis=1), axis=0)
# 行の順序を逆にする
df_normalized = df_normalized.iloc[::-1]

# 横向きの積み上げ棒グラフの描画
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
df_normalized.plot(kind="barh", stacked=True, ax=ax)

# タイトルとラベルの設定(任意)
ax.set_title("地域別の販売製品割合")
ax.set_ylabel("地域")
ax.set_xlabel("割合")
plt.show()

出来上がりがこちら。

はい、サクッとできましたね。

ちなみに、pandasのplotメソッドを使わない場合は棒グラフの各色の部分の左端の位置を指定しながら順番に描くため、次のコードになります。

import numpy as np


# 積み上げ棒グラフを作成するための基準位置
left = np.zeros(len(df_normalized))

fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)


# 各カテゴリごとにバーを追加
for column in df_normalized.columns:
    ax.barh(
        df_normalized.index,
        df_normalized[column],
        left=left,
        label=column,
        height=0.5
    )
    left += df_normalized[column]

# タイトルとラベルの設定
ax.set_title("地域別の販売製品割合")
ax.set_ylabel('地域')
ax.set_xlabel('割合')

# 凡例の表示
ax.legend()

plt.show()

明らかに面倒なコードでしたね。pandasがいい感じに調整してくれていた、棒の太さ(横向きなので今回の例では縦幅)なども自分で調整する必要があります。
Pythonでコーディングしてるのであれば大抵の場合は元データはDataFrameになっているでしょうし、もしそうでなくても変換は容易なので、Pandasに頼った方を採用しましょう。

いつも使っている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

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

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