複数の確率密度関数の重なった部分の面積を計算する

Question

集団Aと集団Bから平均と分散を計算し、ヒストグラムと近似する確率密度関数func Aとfunc Bがわかっています。参考までに両者のグラフイメージをコードで示します。

# locとscaleの値は仮想のもの
y1 = stats.norm.pdf(x,loc=2,scale=0.5)
y2 = stats.norm.pdf(x,loc=3,scale=0.8)
plt.plot(x,y1,label='func A')
plt.plot(x,y2,label='func B')
plt.legend()
plt.show()

やりたいことは上記のグラフの重なり合った部分の面積を計算することです。
まだ、コードにしていませんが、次のようなプロセスで計算します。
(1)fucn Aとfunc Bの交点を推定する。
func A - Func B　という関数を定義して答えとなり得るxを複数生成してこの関数に入力する、
この結果が限りなくゼロに近くなるxを交点のxとする。
(2)func Bの極小値から交点までのxを範囲にfunc Bの面積を計算する。
func Bについてxの極小値から交点までの面積（積分値）を計算する、計算はscipy.integrate.quadを用いて計算する。極小値は、期待値（func Bの結果）が0.0001となるxを答えとなりそうなxを複数生成して結果が最も近似するものを探索することで推測する。
(3)func Aの極大値と交点までのxの範囲を対象にfunv Aの面積を計算する。
計算方式は項番(2)と同じとする。

教えてほしいこと

項番(1)の交点を求める関数はあるのでしょうか？
今は試行錯誤で近似値を計算する方針ですが、できればpythonに実装されている関数でうまく処理したく思っています。

項番(2)および(3)における極小値・極大値を求める関数はあるのでしょうか？
今は上記と同様に試行錯誤で計算する方針ですが、pythonで実装済みの関数を使いたいと考えています。

科学計算系なので、scipyにあるのではないかと想定して公式サイトで探そうとしたのですが、キーワードがわからず検索できていません。よって、関数そのものでなくてもキーワードでも構いません。

Answer 1

どこまでの精度が必要か分かりませんが、提示のコードからだと、xの範囲でy1とy2のリストについて小さい方の値をとったy3というリストを作って面積を求めるのではだめですか？

x = np.linspace(-10,10,100)
y1 = stats.norm.pdf(x,loc=2,scale=0.5)
y2 = stats.norm.pdf(x,loc=3,scale=0.8)
y3 = np.minimum(y1,y2)
plt.plot(x,y1,label='func A')
plt.plot(x,y2,label='func B')
plt.plot(x,y3,label="A*B")
plt.legend()
plt.show()

==追記==
積分範囲はy1,y2の関数で取りうる値域をそのまま使うのではだめですか？
その場合、Xのmax/minを計算しなくても済みますが。

# locとscaleの値は仮想のもの
x = np.linspace(0,10,5000)  #5000点に増やした
y1 = stats.norm.pdf(x,loc=2,scale=0.5)
y2 = stats.norm.pdf(x,loc=3,scale=0.8)
y3 = np.minimum(y1,y2)
abInteg=integrate.simps(y3, x) #AxBの面積を求める

plt.plot(x,y1,label='func A')
plt.plot(x,y2,label='func B')
plt.fill(x,y3,label="A*B")
AxB='AxB area={}'.format(abInteg)
plt.text(5,0.6,AxB)
plt.legend()
plt.show()

Comments

[R.Shigemori]

回答ありがとうございます。
ご提案の方法を使うと交点を求める必要がなさそうですね。（y3を求める関数を定義して、その関数に対して積分計算を行うとうまくいきそうな気がします。）残るは、極小値と極大値の求め方ですね。実は、同じ作業を30ほどの地域で行う必要があるので、試行錯誤は避けたいです。

[退会済みユーザー]

極小値、極大値は微分の符号反転ですかね。

[R.Shigemori]

dkato0077さん
微分というよりも逆関数にニュアンスとしては近いかと思っています。
確率密度関数は、Xを与えると期待値を計算してくれるものなので、期待値を与えるとXを求めることもできるはずです。いろいろと式変換すれば計算式が得られると思うのですが、できればpythonに実装済みの関数を使うことが希望です。

[TaroToyotomi]

ここでいう確率密度関数の極小値と極大値の定義が分かりません。極大値は確率密度関数のピーク値の部分でいいんですよね？ただ、そうなると極小値は正規分布の場合はどこになるんでしょう？

[R.Shigemori]

すみません。言葉の使い方が悪かったようです。積分するにはxの範囲の定義が必要です。x軸についてこの積分するスタート地点を極小値、エンド地点を極大値と深く考えずに呼んでしまいました。どちらとも無限を取り得るので、質問に記載した通り期待値が0.0001になるxが真の定義です

[R.Shigemori]

simps関数は知りませんでした。
確かにこれを使えば難しいことを考えずに答えが得られそうです

Answer 2

(1)fucn Aとfunc Bの交点を推定する。

例えばy = func A - Func Bとしてyの符号に注目する：y_sign = np.sign(y)。符号が反転している辺りを交点とする。np.nonzero(y_sign[1:] - y_sign[:-1])。これを交点とする。