SciPyで重積分

もう結構古い記事なのですが、以前SciPyで定積分をやる方法を記事にしたことがあります。
参考: scipyで定積分

最近、2変数関数の積分をやる機会があったのでこの機会に重積分をSciPyで行う方法を紹介します。SciPyのintegrateモジュールには、重積分用の関数が複数あります。
dblquad (2変数関数の定積分)
tplquad (3変数関数の定積分)
nquad (一般のn変数の定積分)

dblquadの使い方

順番に説明していきます。まずは2重積分のdblquadです。関数の定義は次のようになっています。
scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

必須なのは、積分対象のfunc, 外側の積分区間のa, b、そして内側の積分区間を示す、gfun, hfunです。

gfunとhfunは名前からわかる通り、定数ではなく関数です。これにより内側の積分の積分区間を変数にすることができます。つまり以下のような積分区間の積分ができます。
$$\int_0^1\int_0^y xy \,dxdy$$

例に挙げたのでこれを実装してみましょう。ちなみに解は$1/8=0.125$です。funcの定義は、内側の関数を第1引数にする必要があるので注意してください。

from scipy.integrate import dblquad


def f1(x, y):
    return x*y


def x0(y):
    return 0


def x1(y):
    return y


print(dblquad(f1, 0, 1, x0, x1))
# (0.125, 5.515032205777789e-15)

想定通りですね。積分結果と推定誤差が返ってくるのは1変数の積分と同様です。

内側の積分区間も定数から定数までだよ、要するに長方形領域で積分したいよ、って場合はgfun, hfun に定数を返す関数を返してください。

tplquadの使い方

続いて、3変数向けのtplquadです。これもdblquadとかなり似てる感じで使えます。積分変数が一個増えているので上限加減の指定がもう一個ある感じです。
scipy.integrate.tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

たとえば次の積分をやってみましょう。

$$\int_0^1\int_0^z\int_0^{y+z} xyz \,dxdydz.$$

ちなみに答えは$17/144=0.11805555…$となるはずです。

引数の順番に注意が必要なので慎重にコーディングしてください。

from scipy.integrate import tplquad


def f2(x, y, z):
    return x*y*z


def x0(y, z):
    return 0


def x1(y, z):
    return y + z


def y0(z):
    return 0


def y1(z):
    return z


print(tplquad(f2, 0, 1, y0, y1,  x0, x1))
# (0.11805555555555557, 2.1916761217856673e-14)

バッチリですね。

nquadの使い方

最後に一般のn変数を積分できるnquadの使い方を紹介します。

引数の形式が先ほどの二つと少し違います。
scipy.integrate.nquad(func, ranges, args=None, opts=None, full_output=False)

funcにn変数関数を渡して、rangesに積分区間を渡すことになります。rangesは配列で、1変数目から順番に区間の下限上限の2値の配列を格納しておけば良いです。また、ここにも一応関数を使うことはできます。

これはシンプルな例で、定数関数1を超立方体区間で積分してみました。

from scipy.integrate import nquad


def f3(w, x, y, z):
    return 1

print(nquad(f3, [[-2, 2], [-2, 2], [-2, 2], [-2, 2]]))
# (256.0, 2.8421709430404007e-12)

$4^4=256$になりましたね。

ここで急にシンプルな例を出したのには事情がありまして、変数の数が多くなるとやはり積分は困難なようで、ちょっと複雑な例になると規定の反復回数をこなしても必要な精度に届かずWarningが出たりするケースが多々あります。

どうしても計算したい場合は limit パラメーター等をいじっての対応になりますのでドキュメントを参照しながら調整してみてください。
(僕も実運用で必要になったら改めて調査して紹介しようと思います。)

Pythonの関数から一部の引数を固定して新しい関数を作る

Pythonの多くのライブラリの様々な関数が非常に汎用的に使えるように作られているので多くの引数を受け取れるようになっています。しかし、そのほとんどの引数を固定して1変数関数として使いたいなぁと思うようなことがあります。PandasのDataFrameのapplyなど関数を引数として受け取る関数に渡す場合等ですね。
また、大量にある引数のほとんどを固定して一部だけ変えながら何度も実行する、といった場面も考えられます。

lambda式などを作ってラップした新しい関数を実装してもいいのですが、 functoolsという標準ライブラリにその専用のpartial というメソッドが用意されています。
参考: functools.partial(func/*args**keywords)

このpartialを使うと、引数の一部を固定した引数の少ない新しい関数を作ってくれます。

一個目の引数に元になる関数を渡し、2個目以降の引数に渡したものが、元の関数の固定引数として使われます。keyword引数で渡せばそのkeyword引数が固定されます。

一引数の固定の方は先頭から順番に固定されるので注意してください。つまり2番目以降の引数を固定したい場合はそれらはキーワード引数として指定する必要があります。

サンプル

引数を順番に表示するだけの単純な関数を作ってやってみましょう。

from functools import partial


# 3つの引数を表示するだけの関数
def sample_func(a, b, c):
    print("a=", a)
    print("b=", b)
    print("c=", c)


# テスト実行
sample_func(1, 2, 3)
"""
a= 1
b= 2
c= 3
"""

# a = 10, b = 20 を固定した新しい関数が作られる。
partial_f = partial(sample_func, 10, 20)


# 3個目の引数 c = 50だけ渡して実行できる。
partial_f(50)
"""
a= 10
b= 20
c= 50
"""

# キーワード引数で固定することもできる。
partial_f2 = partial(sample_func, a=100, c=200)

# b の値だけ渡して実行できる
partial_f2(b=-5)
"""
a= 100
b= -5
c= 200
"""

キーワード引数を固定した関数を、位置引数で使う場合は注意が必要です。
たとえば、次のようにaを固定して生成した関数に、残り2個の引数を位置引数で渡すと、aを2回渡した扱いになってエラーが起きます。

# aを固定
partial_f3 = partial(sample_func, a=1)

# bとcのつもりで残り2個の引数を渡すとエラー
try:
    partial_f3(2, 3)
except Exception as e:
    print(e)
# sample_func() got multiple values for argument 'a'

# bとcもキーワード引数で渡す。
partial_f3(b=2, c=3)
"""
a= 1
b= 2
c= 3
"""

まとめ

ほぼ小ネタのような内容でしたが、自作関数をベースに一部の振る舞いを固定した簡易的な関数を作るとか、apply等の1変数関数を受け取るメソッドに渡したいとかそういう場面で役に立つことがあるテクニックとしてpartialを紹介しました。

scipyのstats配下の各種メソッドであれば、それぞれがパラメーターを固定するfrozenメソッドを持ってるとか、引数が多いなら引数を辞書にまとめて**(アスタリスク2個)で展開すればいいとか、ラップした関数を自分で実装したらいいとか、代用手段も多いのですが、partialを使うとその辺の記述がシンプルになるので機会があれば使ってみてください。

SciPyでニュートン法を利用する

前回の記事の二分法に続いて、もう一つ求根アルゴリズムを紹介します。
参考: 二分法を用いて関数の根を求める

今回紹介するのはニュートン法です。これは微分可能な関数$f(x)$の根を求めることができるアルゴリズムです。
参考: ニュートン法 – Wikipedia

詳しい説明は上記のWikipediaにあるので、ざっくりと概要を説明します。

この方法の背景にあるのは、滑らかな関数をある点の近くだけ着目してみるとほぼ直線になり、接線で近似できるということをベースのアイデアにしています。

つまり、微分可能な関数$f$があって$f(x)=0$だとします。その根$x$の近くに点$x_0$を取ると、$x_0$の近くでは、$f$と$f$の接線ってかなり近いよね、それなら$f$の根と$f$の$x_0$における接線の根って近いよね、っていうのが基本的なアイデアです。

関数$f$の$f(x_0)$における接戦は次の式で書けます。

$$y=f'(x_0)(x-x_0)+f(x_0).$$

$f(x)=0$は解けない場合でも、この接線の根は容易に算出することができ、

$$x_0-\frac{f(x_0)}{f'(x_0)}$$

と求まります。

この値は元の$x_0$よりも真の根に近いことが期待され、これをもう一回$x_0$とおいて同じ操作を繰り返せば真の根にたどり着く、というのがニュートン法です。

Wikipediaから画像拝借しますが、図で見るとイメージしやすいですね。

注意しないといけないのは、初期値$x_0$は真の根$x$の十分近くに取らないといけない点です。十分近くを見れば関数をその接線で近侍できるよね、というのがアイデアの前提なので、根が近くになかったらその前提が崩れてしまいこのアルゴリズムは真の根に収束しなくなってしまいます。

ニュートン法のメリットとデメリット

先に紹介した二分法と比べて、ニュートン法のメリットデメリットを説明していきます。

1番のメリットは収束の速さです。二分法に比べてより少ない計算回数で効率的に会を探索することができます。

また、初期値として与える点が1点だけで良いというのもメリットです。二分法の場合は初期値は区間で設定する必要がありましたからね。

その一方で複数のデメリットもあります。実装していて一番不便に感じるのはその関数だけでなく微分も必要ということでしょうか。もちろん微分不可能な関数ではニュートン法は使えません。

また、初期値が真の解に十分近くない場合や、微分した値が$0$に近い場合、うまく収束せずにアルゴリズムが失敗してしまう、という点も大きなデメリットです。

SciPyによる実装

SciPyではscipy.optimizeというモジュールで実装されています。newtonという専用メソッドを使うか、root_scalarという汎用的なメソッドで(method=’newton’)を指定して使うことになります。二分法と同じですね。

参考:
scipy.optimize.newton — SciPy v1.12.0 Manual
root_scalar(method=’newton’) — SciPy v1.12.0 Manual

二分法の時と同じように、$\sin$関数の根$\pi$を探索させてみましょう。微分は$\cos$なのでこれを使います。

from scipy import optimize
import numpy as np


root1 = optimize.newton(np.sin, x0=3, fprime=np.cos)
print(root1)
# 3.141592653589793

root_result = optimize.root_scalar(np.sin, method="newton", x0=3, fprime=np.cos)
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 6
     iterations: 3
           root: 3.141592653589793
"""

print(root_result.root)
# 3.141592653589793

簡単ですね。

注目するのは、iterationsの部分です。たった3回のイテレーションで収束していて、関数が実行されたのは、fとfの微分合わせて6回だけです。
二分法の時は39回もイテレーションが必要だったのと大違いです。そして実はこの例では解の精度もニュートン法の方が高くなっています。

ニュートン法が失敗する例

初期値が真の解の近くにないと失敗するという話がありましたのでそちらも見ておきます。

例えば、タンジェントの逆関数、$\arctan$で試してみましょう。(sin, cosは根が無限にあって、根から遠い実数を用意できないので関数を変えます。)

$\arctan(x)$の微分は$\frac{1}{1+x^2}$です。

やってみました。

def f(x):
    return np.arctan(x)


def fprime(x):
    return 1/(1+x**2)


# 初期値が1なら収束する
root_result_1 = optimize.root_scalar(np.sin, method="newton", x0=1, fprime=fprime)
print(root_result_1)
"""
      converged: True
           flag: 'converged'
 function_calls: 12
     iterations: 6
           root: 0.0
"""

# 初期値が2だと失敗し、結果のflagが'convergence error'になる。
root_result_2 = optimize.root_scalar(np.sin, method="newton", x0=2, fprime=fprime)
print(root_result_2)
"""
      converged: False
           flag: 'convergence error'
 function_calls: 100
     iterations: 50
           root: 1.854706857103781
"""


# optimize.newton の方だと例外が上がる。
try:
    optimize.newton(f, fprime=fprime, x0=2)
except Exception as e:
    print(e)
# Derivative was zero. Failed to converge after 10 iterations, value is -6.999943395317963e+168.

失敗した時の振る舞いがそれぞれ違うので、どちらのコードを使うかで注意深く扱う必要がありますね。optimize.root_scalarはコード自体は正常に終了しますがフラグが立ち、optimize.newtonの方は例外があがります。

ちなみに、エラーの中で出てくるvalue の値、 -6.999943395317963e+168 は次のように自分でニュートン法を実装しても同じ値が出て来ます。

x0 = 2  # 初期値
for i in range(12):
    x0 = x0 - f(x0)/fprime(x0)
    print(i+1, "回目: x0=", x0)

"""
1 回目: x0= -3.535743588970453
2 回目: x0= 13.950959086927496
3 回目: x0= -279.34406653361754
4 回目: x0= 122016.9989179547
5 回目: x0= -23386004197.933937
6 回目: x0= 8.590766671950415e+20
7 回目: x0= -1.1592676698907411e+42
8 回目: x0= 2.110995587611039e+84
9 回目: x0= -6.999943395317963e+168
10 回目: x0= inf
11 回目: x0= nan
12 回目: x0= nan
"""

絶対値が大きくなり続けていて全く収束に向かっていないのがわかりますね。

まとめ

元の関数だけではなく導関数も必要だったり、初期値の設定段階である程度解の目星をつけておかないといけないなどのデメリットはありますが、速度や精度の面で優秀でしかもロジックもわかりやすい手法なので、何か機会があればニュートン法の活用を検討してみてください。

二分法を用いて関数の根を求める

1変数連続関数の根(値が0になる点)を求める、二分法というアルゴリズムとそれをScipyで実装する方法を紹介します。

アルゴリズムの内容

二分法というのは中間値の定理をベースとした求根アルゴリズムです。アイデアは非常に単純で、連続関数$f$に対して、$f(x_1)$と$f(x_2)$の符号が異なるように、$x_1, x_2$を選びます。この時点で、中間値の定理より区間$(x_1, x_2)$に根があることがわかりますのでさらに細かくみていきます。次は、$x_1, x_2$の中点$x_M = \frac{x_1+x_2}{2}$を取り、$f(x_M)$の符号を調べます。$f(x_M)$と$f(x_1)$の符号が同じであれば、$x_1$を$x_M$で置き換え、逆に$f(x_M)$と$f(x_2)$の符号が同じであれば、$x_2$を$x_M$で置き換えると、区間の幅が半分になった$(x_1, x_2)$が得られますが根はこの中にあることがわかります。これを繰り返すと、根が存在する範囲を狭めていくことができ、$f(x_M)$の絶対値が0になるか、もしくは十分0に近づいたらその値を数値的に求めた根とします。

以上が、一般的な二分法のアルゴリズムの説明です。ただし、後に紹介するSciPyではどうやら区間が十分狭くなったかイレーション回数が上限に達したか等の基準でループを打ち切り、$f(x)$の値を確認していないようです。

二分法のメリットとデメリット

方法が単純でわかりやすい、というのが個人的に感じている1番のメリットです。

また、連続関数であれば使えるため、ニュートン法などのアルゴリズムのように元の関数の微分を必要とせず、微分が難しい関数や微分不可能な関数でも使えます。

また、根が存在しうる区間を狭めながら探索するため、最初の区間の幅と繰り返し回数により、結果の精度を保証できることも大きな利点です。

逆にデメリットとしては、ニュートン法等と比較して収束が遅いこととか、初期値として関数が異符号になる2点を探して与える必要があること、もしその区間に複数の根が存在した場合にどれに収束するか不確定なことなどが挙げられます。

ただし僕の経験上では、ある程度根の目処がついていたり、単調な関数に対して使う場面が多くこれらのデメリットを深刻に感じることは少ないです。

SciPyによる実装

SciPyではscipy.optimizeというモジュールで実装されています。bisectという専用メソッドを使うか、root_scalarという汎用的なメソッドで(method=’bisect’)を指定して使うことになります。

参考:
scipy.optimize.bisect — SciPy v1.11.4 Manual
root_scalar(method=’bisect’) — SciPy v1.11.4 Manual

試しに、$\sin$関数の 3と4の間にある根を探させて見せましょう。既知の通りそれは円周率$\pi$になるはずです。

from scipy import optimize
import numpy as np


root1 = optimize.bisect(np.sin, 3, 4)
print(root1)
# 3.1415926535901235

root_result = optimize.root_scalar(np.sin, bracket=[3, 4], method='bisect')
# 結果は複数の情報を含むRootResults形式で戻る。
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 3.1415926535901235
"""
print(type(root_result))
# class 'scipy.optimize._zeros_py.RootResults'

# 根の値へのアクセス方法
print(root_result.root)
# 3.1415926535901235

いかにも円周率ぽい結果が得られましたね。root_scalarの方では収束したことを示すフラグや、イテレーション回数なども得られています。

失敗事例1. 区間の両端での関数の値が同符号の場合

二分法は初期設定を誤ってると失敗するのでその場合のSciPyの挙動も見ておきましょう。失敗パターンの一つは、最初に指定した区間の両端で符号が一致していた場合です。もちろん関数の形によってはその区間内に根がある可能性もあるのですが、存在は保証されなくなります。

また$\sin$関数で、その間に根が存在しない区間$(1, 2)$と、実は両端で同符号だけど根が存在する区間$(1, 7)$でやってみましょう。bisectとroot_scalarで全く同じエラーメッセージ出るのでbisectの方だけ載せます。

try:
    optimize.bisect(np.sin, 1, 2)
except Exception as e:
    print(e)
# f(a) and f(b) must have different signs

try:
    optimize.bisect(np.sin, 1, 7)
except Exception as e:
    print(e)
# f(a) and f(b) must have different signs

はい、$f(a)$と$f(b)$は違う符号にせよとのことでした。根が存在しない1個目の例はさておき、2個目の例は根は区間内に2個存在するのですが探さずにエラーになりました。

失敗事例2. 連続関数ではなかった場合

もう一つ失敗するのは、関数が連続関数ではないケースです。

例えば$\tan$の、$\frac{\pi}{2}$近辺の挙動で見てみましょう。数学的に厳密な話をすると、$\frac{\pi}{2}$では$\tan$は定義されないので、$\tan$は数学的には連続関数(定義域内のすべての点で連続)なのですが、数値計算的には不連続と考えた方が都合が良いです。

話が脇に逸れたので実例の話に移ります。実は、bisectメソッドは結果を返して来ちゃうんですよね。そしてそれが全然根ではないということも見ておきましょう。

root = optimize.bisect(np.tan, 1, 2)
# pi/2に近い結果が得られている
print(root)
# 1.5707963267941523

# 元の関数に代入した結果は全く0に近くない
np.tan(root)
# 1343445450736.3804

root_scalarの方だったら、結果のフラグ等もあるのでアラート等あげてくれるのかと期待したのですがそういう機能はなさそうです。

root_result = optimize.root_scalar(np.tan, bracket=[1, 2], method='bisect')
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 1.5707963267941523
"""
np.tan(root_result.root)
# 1343445450736.3804

要するに、SciPyに渡す関数が本当に連続関数であるかどうかは利用者が責任を持たないといけないということです。また、結果が本当に根なのかどうかは代入して確認した方が良いでしょう。

まとめ

以上が二分法の説明とSciPyで利用する方法、その注意点でした。

argparseで引数を受け取る

はじめに

今週の記事もPythonスクリプトで引数を受け取って使う話です。前回はsys.argvつかって受け取る方法を紹介していましたが、今回は便利な専用モジュールのargparseを紹介します。

これを使うと、引数を変数に自動的に格納したり、オプション引数やフラグを作成したり、ヘルプ機能を自動的に作ってくれたりします。

よくUnix/Linux コマンドでは -o filename みたいな感じで出力先ファイルを指定できたりしますが、これをsys.argvで実装しようとすると、配列を全部見て-oがあるかどうか確にして、その次の値をfilenameとして取得して、みたいな結構面倒な処理を自分で作る必要があります。-oが複数出てきたらどうするかとか、-oの次にファイル名がなかった場合のハンドリングとか色々考えないといけないのでとても面倒です。こういう手間を削減してくれます。

順に使い方書いていきますが、ドキュメントはこちらです。
参考: argparse — コマンドラインオプション、引数、サブコマンドのパーサー — Python 3.12.1 ドキュメント

基本的な使い方

ざっくりいうと、argparseは次の3手順で使います。

  • ArgumentParserオブジェクトの作成
  • 必要な引数をパーサーオブジェクトに追加する
  • 引数を解析して結果を取得する

一回単純なサンプル作ってやってみましょう。上記ドキュメントの例をそのまま使います。
sample.py というファイル名で次のスクリプトを作成し、実行権限を `$ chmod u+x sample.py` でつけておきます。

#!/usr/bin/env python
import argparse


# パーサーオブジェクトの作成
parser = argparse.ArgumentParser(
    prog="ProgramName",
    description="What the program does",
    epilog="Text at the bottom of help"
)

# 必要な引数の追加
parser.add_argument("filename")  # 位置引数
parser.add_argument("-c", "--count")  # 値を取るオプション 
parser.add_argument("-v", "--verbose", action="store_true")  # on/off フラグ

# 引数の解析
args = parser.parse_args()
print(args.filename, args.count, args.verbose)

だいたいイメージできると思うのですが、./sample.py を実行する時、最初の引数が filename に格納されて、 -c か –count で指定した値が count変数に格納され、 -v を選択したかどうかがTrue/False で verbos に入ります。 色々やってみましょう。

$ ./sample.py test.txt
test.txt None False

$ ./sample.py test.txt -c 4 -v
test.txt 4 True

$ ./sample.py -v --count abc test.txt
test.txt abc True

$ ./sample.py test1.txt test2.txt 
usage: ProgramName [-h] [-c COUNT] [-v] filename
ProgramName: error: unrecognized arguments: test2.txt

$ ./sample.py                    
usage: ProgramName [-h] [-c COUNT] [-v] filename
ProgramName: error: the following arguments are required: filename

$ ./sample.py --help
usage: ProgramName [-h] [-c COUNT] [-v] filename

What the program does

positional arguments:
  filename

options:
  -h, --help            show this help message and exit
  -c COUNT, --count COUNT
  -v, --verbose

Text at the bottom of help

はい、最初の3例が正しくコマンドを打ったケースでしたが、だいたいイメージ通りに引数を受け取れていることが確認できると思います。
4つ目は位置引数を過剰に設定、5つ目は逆に指定しませんでしたが、それぞれちゃんとエラー文を出してくれていますね。
6個目の例は–helpをつけていますが、なんと自動的にヘルプメッセージを作成して表示してくれています。

コマンド名が ProgramName になっていますが、これはパーサーを作成したときのprog 引数をプログラム名として使っているからです。progを省略すると、ファイル名が使われます。

これは大事なことなのですが、プログラム名=ファイル名のことが多いと思うので、基本的に省略した方がいいと思います。(さっきの例は公式ドキュメントをただ真似しただけ。)

descriptionでプログラム中身の説明、epilogでヘルプの最後に表示するメッセージを指定できますが、これらもどちらも省略可能です。ただ、descriptionは何か書いていておいた方がいいと思います。

ここから細かく仕様を見ていきます。

引数の種類

引数の種類としては、コマンドの後に何番目に渡されたかどうかで扱いが決まる位置引数と、-(ハイフン)付きの名前で始まるオプション引数があります。

argparseは接頭辞の”-“を特別な文字として扱って、これによって挙動を変えています。

上の例でもわかりますが、次のように複数の名前を指定することもできますし、1種類だけの名前でも良いです。このとき注意しないといけないのは、参照するときの変数名です。

-c みたいな短い名前だけの時は – をとって c として参照しますが、-c, –count という2種類の名前を指定した場合は、一番最初に登場する長い名前、が採用されます。長い名前というのは文字列の長さの話ではなく、 – ではなく、 — で始まる引数ということです。
つまり、次のように3つの名前をつけたら、長い名前の中で最初に登場した countが採用されるということです。

#!/usr/bin/env python
import argparse


# パーサーオブジェクトの作成
parser = argparse.ArgumentParser()

# 必要な引数の追加
parser.add_argument("-b")
parser.add_argument("-c", "--count", "--cnt")

# 引数の解析
args = parser.parse_args()
# 短い名前 -b しかないのでbでアクセス
print(args.b)

# --count と --cnt が長い名前だが、先に登場したcountの方が優先
print(args.count)

add_argument の引数

add_argument には様々な引数を指定でき、各種の設定を行うことができます。
全部紹介するのも大変なので一部抜粋して紹介しますが、公式ドキュメントの該当欄の一読をお勧めします。
参考: add_argument() メソッド

default ・・・ コマンドラインに対応する引数が存在せず、さらに namespace オブジェクトにも存在しない場合に利用されるデフォルト値。
type ・・・ データ型。 int や float、ユーザー定義の型など色々指定できる。省略すると文字列(str)。
choices – 引数として許される値のシーケンス。
help ・・・ 引数の説明。-h や –help を使用した時に使われる。
nargs ・・・ 受け取れるコマンドライン引数の数。後で説明します。
action ・・・ コマンドラインにこの引数があったときの動作。後で説明します。

だいたいはイメージ通りの挙動をしてくれるのですが、nargsとactionについてはこの後説明します。

受け取れるコマンドライン引数の数について

nargs という値を使って、受け取れる引数の数を指定できます。

nargs の指定は正規表現風になっています。 整数Nを指定すればその個数、?なら1個か0個で、0だったらdefalut値が使われます。*とすると任意の数受け取れます。また、+だとこちらも任意の数受け取れますが、0個だった場合にエラーが起きます。

例えば、任意の数のファイルのデータを入力として、1個のファイルに結果を書き出すようなコマンドがあったとしましょう。(というより、tar コマンドでアーカイブ作る時はそういう指定しますよね。 )

次のような形です。

#!/usr/bin/env python
import argparse


# パーサーオブジェクトの作成
parser = argparse.ArgumentParser()

# 必要な引数の追加
parser.add_argument("-o", "--out_file")  # nargsを省略しているので1
parser.add_argument("-i", "--in_file", nargs="+")

# 引数の解析
args = parser.parse_args()
print(args.in_file)
print(args.out_file)

# 以下実行結果
$ ./sample.py -o out.txt  -i in1.txt in2.txt in3.txt
['in1.txt', 'in2.txt', 'in3.txt']
out.txt

— in_file の方は複数の結果を受け取れるようにしたので、Python上は配列で結果が来るようになりましたね。

actionによる動作の指定について

actionを使って、オプション引数が存在したときの挙動を指定できます。

デフォルトは store でこれは要するに変数を値に格納するという挙動です。さっきまで見てるのがこれですね。

ただし、Linux/Unixコマンドではこのような値を受け取る引数ばかりではありません。皆さんがよく使う ls コマンドの -l とか -a は別に何か引数を受け取ったりせず、その存在の有無だけが重要ですよね。

この記事の冒頭のコードの `parser.add_argument(“-v”, “–verbose”, action=”store_true”) # on/off フラグ`
もまさにそうで、 -v の有無だけが問題になります。これを実現しているのが、action=”store_true”の部分です。

要するに -v が見つかったら verboseにTrueを格納するよ、という挙動になります。
逆に見つからなかったらFalseが格納されます。

これと逆にオプションがあったらFalseでなかったらTrueになるのが、”store_false”です。

このほか、キーワード引数の登場回数を数えて格納する”count”とか、複数回登場したら結果を都度配列に追加していく”append”などもあります。

これらも一通り公式ドキュメントの一読をお勧めします。
参考: action

ヘルプの作成について

最後にヘルプ機能についてです。自動的に、-h と –helpがヘルプ機能として実装されます。

これはもう実際に試していただくのが一番早いのですが、description等で指定されたプログラムの説明や、受け取れるコマンドライン引数の情報などが表示でき大変便利です。

気をつけないといけないのは、 -h と –help を上書きしないようにすることですね。もちろんどうしてもこれらの引数名を別用途で使いたいとか、自作のヘルプメッセージを実装したいとか事情があれば話は別ですが、普通はデフォルトのヘルプを使った方が良いと思います。

argparseをコマンドライン引数以外の文字列のパースに使う

最後にちょっとマニアックな使い方を紹介します。
このargparseですが、何も指定せずに、 parser.parse_args() すると コマンドライン引数をパースしにいきますが、ここで配列を渡すとその配列をパースします。

sample_str = “-i filename -c 5” みたいな文字列があった時にsample_str.split()して配列に分解して、 parser.parse_args(sample_str.split())と渡すとそれをコマンドライン引数と見立ててパースしてくれるのです。

そんな技術いつ使うねん、と思われるかもしれませんが、僕はマジックコマンドを作る時などに使ってます。
参考: Snowflakeに手軽にSQLを打てるJupyterマジックコマンドを作ってみた|ホンディー | ライフイズテック  (このブログ書いてる人のnote記事です。)

これをやると、parse_args は渡された配列をパースしてるのでコマンドラインから渡した引数は全部無視する点には注意してください

まとめ

長くなりましたが、以上がargparseの説明になります。argparseはシンプルに利用することもできますし、多くの引数を活用して細かいカスタマイズもでき、大変柔軟にツールを作ることができます。 自前ツールを作成する際の大変有益な武器になりますので是非触ってみてください。

Pythonファイルをコマンドラインで実行したときの引数の受け取り方。(sys.argvを使う場合)

久々にコマンドラインツールを作った時に、引数の扱い方をド忘れしてしまっていたのでそのメモです。

ちなみに、argpaseっていう大変便利なライブラリもあるのでそれの使い方も紹介したいのですが、それは別記事に回して今回はもっとシンプルな方について書きます。

記事タイトルにも書いていますが、基本的には、sys.argvというのを使います。これは関数ではなく配列型のデータで、argv[0]にそのスクリプトのファイル名、argv[1]以降に、コマンドライン引数が入ります。

参考: sys — システムパラメータと関数 — Python 3.12.1 ドキュメント

また、Pythonのバージョン3.10 以降では、sys.orig_argv ってのも追加されています。

動かしてみるのが一番確認しやすいので、やってみましょう。
次のようなファイルを作ります。ファイル名は sample.py としました。

import sys


for i in range(len(sys.argv)):
    print(f"sys.argv[{i}]=", sys.argv[i])

for i in range(len(sys.orig_argv)):
    print(f"sys.orig_argv[{i}]=", sys.orig_argv[i])

ではこれを実行してみましょう。一旦実行権限とかつけてないので、明示的にpythonコマンドとして実行します。

引数はただ表示するだけで何も処理に影響はないので、デタラメにつけてみました。

% python sample.py file.txt -l aaa.txt "bbb" --name=xyz
sys.argv[0]= sample.py
sys.argv[1]= file.txt
sys.argv[2]= -l
sys.argv[3]= aaa.txt
sys.argv[4]= bbb
sys.argv[5]= --name=xyz
sys.orig_argv[0]= /Users/{macのユーザー名}/.pyenv/versions/3.11.1/bin/python
sys.orig_argv[1]= sample.py
sys.orig_argv[2]= file.txt
sys.orig_argv[3]= -l
sys.orig_argv[4]= aaa.txt
sys.orig_argv[5]= bbb
sys.orig_argv[6]= --name=xyz

はい、argvの方は、[0]番目にスクリプト名が表示され、[1]以降に引数が格納されていましたね。”bbb”のところに注目ですが、ダブルクオーテーションは外されています。

orig_argv は、ファイル名の前の、pythonパスまで入っています。the original command line arguments というから、コマンドに打った python がそのまま入ってるかと思ったらフルパスに展開されています。

ちなみに、sys.argv[0]= sample.py と sys.orig_argv[1]= sample.py のスクリプトパスですが、これはコマンドラインで打ったファイルパスがそのまま表示されいます。どういうことかというと、現在はカレントディレクトリで実行したからこのように表示されているだけで、違う場所から相対パスや絶対パスで指定して実行したらこの中身は変わります。フルパスで指定したらフルパスが入ります。

続いて、このファイルに実行権限をつけて頭のpythonを外してみましょう。
ファイルの先頭に、 #!/usr/bin/env python のシバンを挿入して、chmod で実行権限つけときます。

これで実行すると次のようになります。(引数少し減らしました)

% ./sample.py file.txt -l aaa.txt
sys.argv[0]= ./sample.py
sys.argv[1]= file.txt
sys.argv[2]= -l
sys.argv[3]= aaa.txt
sys.orig_argv[0]= /Users/{macのユーザー名}/.pyenv/versions/3.11.1/bin/python
sys.orig_argv[1]= ./sample.py
sys.orig_argv[2]= file.txt
sys.orig_argv[3]= -l
sys.orig_argv[4]= aaa.txt

sys.argv の方はさっきと大して変わらないですね。カレントディレクトリにパスが通ってないので、スクリプトの指定が./ファイル名 になったのが反映されています。

orig_argv の方も先の結果とほぼ変わりませんが、これはちょっと意外でした。the original command line arguments というから、コマンドラインで python 省略したらここからも省略されると思ってたら相変わらずpython本体のパスが登場しています。

orig_argv が何を意図して追加されたのかがよくわからないのですが、とりあえず argv の方を使っておけば良さそうです。

[1]以降の引数たちは冒頭でも書いたargpaseってライブラリで取ることが多く、sys.argvで取得するのはよほど単純なスクリプトの場合に限られるかなとも思います。となると、もしかしたらargv[0]のファイル名の方がよく使うかもしれないですね。

おまけ

記事の本題と逸れるのですが、argv[0]でスクリプト名が取れるとしたら、スクリプトのフルパスやカレントディレクトリの取得方法が気になるかもしれないのでその撮り方もメモしておきます。

結論を言うと、ファイル名は __file__ って特別な変数に格納されています。
また、カレントディレクトリは os.getcwd() で取れます。

Pythonで100%積み上げ棒グラフを描く

個人的な好みの話ですが、僕は割合の可視化は円グラフより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に頼った方を採用しましょう。

InteractiveShell オブジェクトを利用してマジックコマンドを実行したり変数を操作したりする

いつも使っている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個しかなくてそれを否定したら簡単だってことを頭の片隅にでも置いといてもらえると幸いです。

Pythonによる各種エントロピーや相互情報量の計算

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

エントロピーも相互情報量も数式はそこまで難しくないので、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次元に拡張したやつですね。(条件付きエントロピーではないので注意してください、

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

b1b2
a141
a223

結合エントロピーの場合はですね、元のカウントデータを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_truelabels_predcontingency=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

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

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