<p>你可以用@duhalme建议的答案来求交,然后用这个点来定义积分极限的范围</p>
<p><a href="https://i.stack.imgur.com/7xtWy.png" rel="noreferrer"><img src="https://i.stack.imgur.com/7xtWy.png" alt="enter image description here"/></a></p>
<p>这里的代码看起来像</p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
norm.cdf(1.96)
def solve(m1,m2,std1,std2):
a = 1/(2*std1**2) - 1/(2*std2**2)
b = m2/(std2**2) - m1/(std1**2)
c = m1**2 /(2*std1**2) - m2**2 / (2*std2**2) - np.log(std2/std1)
return np.roots([a,b,c])
m1 = 2.5
std1 = 1.0
m2 = 5.0
std2 = 1.0
#Get point of intersect
result = solve(m1,m2,std1,std2)
#Get point on surface
x = np.linspace(-5,9,10000)
plot1=plt.plot(x,norm.pdf(x,m1,std1))
plot2=plt.plot(x,norm.pdf(x,m2,std2))
plot3=plt.plot(result,norm.pdf(result,m1,std1),'o')
#Plots integrated area
r = result[0]
olap = plt.fill_between(x[x>r], 0, norm.pdf(x[x>r],m1,std1),alpha=0.3)
olap = plt.fill_between(x[x<r], 0, norm.pdf(x[x<r],m2,std2),alpha=0.3)
# integrate
area = norm.cdf(r,m2,std2) + (1.-norm.cdf(r,m1,std1))
print("Area under curves ", area)
plt.show()
</code></pre>
<p>cdf用于获得高斯函数的积分,尽管可以定义高斯函数的符号形式并使用<code>scipy.quad</code>(或其他)。或者,您可以使用类似于<a href="http://rpsychologist.com/calculating-the-overlap-of-two-normal-distributions-using-monte-carlo-integration" rel="noreferrer">link</a>的蒙特卡罗方法(即生成随机数并拒绝任何超出您所需范围的随机数)。</p>