Вероятность перекрытия двух нормальных распределений с scipy

Question

У меня есть две scipy.stats.norm (средняя, ​​std) .pdf (0) нормальная кривая распределения, и я пытаюсь выяснить дифференциальную (перекрывающуюся) две кривые.

Как рассчитать его с помощью scipy в Python? Благодаря

Answer 1

Вы можете использовать ответ, предложенный @duhalme, чтобы получить пересекаются, а затем использовать эту точку, чтобы определить диапазон интегральных пределов,

Если код это выглядит,

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() 

Здесь cdf используется для получения интеграла от гауссова, хотя символическая версия гаусса может быть определена и используется scipy.quad (или что-то еще). В качестве альтернативы вы можете использовать метод Монте-Карло, подобный этому: link (т. Е. Генерировать случайные числа и отклонять любые вне диапазона, который вы хотите).

Answer 2

Ответ Эд замечательный. Однако я заметил, что это не работает, когда есть две или бесконечные (полностью перекрывающиеся) точки соприкосновения. Вот версия кода, который обрабатывает эти два случая.

Если вы также хотите продолжить просмотр графиков распределений, вы можете использовать код Эд.

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.stats import norm 

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 

result = solve(m1,m2,std1,std2) 
# 'lower' and 'upper' represent the lower and upper bounds of the space within which we are computing the overlap 
if(len(result)==0): # Completely non-overlapping 
    overlap = 0.0 

elif(len(result)==1): # One point of contact 
    r = result[0] 
    if(m1>m2): 
     tm,ts=m2,std2 
     m2,std2=m1,std1 
     m1,std1=tm,ts 
    if(r<lower): # point of contact is less than the lower boundary. order: r-l-u 
     overlap = (norm.cdf(upper,m1,std1)-norm.cdf(lower,m1,std1)) 
    elif(r<upper): # point of contact is more than the upper boundary. order: l-u-r 
     overlap = (norm.cdf(r,m2,std2)-norm.cdf(lower,m2,std2))+(norm.cdf(upper,m1,std1)-norm.cdf(r,m1,std1)) 
    else: # point of contact is within the upper and lower boundaries. order: l-r-u 
     overlap = (norm.cdf(upper,m2,std2)-norm.cdf(lower,m2,std2)) 

elif(len(result)==2): # Two points of contact 
    r1 = result[0] 
    r2 = result[1] 
    if(r1>r2): 
     temp=r2 
     r2=r1 
     r1=temp 
    if(std1>std2): 
     tm,ts=m2,std2 
     m2,std2=m1,std1 
     m1,std1=tm,ts 
    if(r1<lower): 
     if(r2<lower):   # order: r1-r2-l-u 
      overlap = (norm.cdf(upper,m1,std1)-norm.cdf(lower,m1,std1)) 
     elif(r2<upper):   # order: r1-l-r2-u 
      overlap = (norm.cdf(r2,m2,std2)-norm.cdf(lower,m2,std2))+(norm.cdf(upper,m1,std1)-norm.cdf(r2,m1,std1)) 
     else:     # order: r1-l-u-r2 
      overlap = (norm.cdf(upper,m2,std2)-norm.cdf(lower,m2,std2)) 
    elif(r1<upper): 
     if(r2<upper):   # order: l-r1-r2-u 
      print norm.cdf(r1,m1,std1), "-", norm.cdf(lower,m1,std1), "+", norm.cdf(r2,m2,std2), "-", norm.cdf(r1,m2,std2), "+", norm.cdf(upper,m1,std1), "-", norm.cdf(r2,m1,std1) 
      overlap = (norm.cdf(r1,m1,std1)-norm.cdf(lower,m1,std1))+(norm.cdf(r2,m2,std2)-norm.cdf(r1,m2,std2))+(norm.cdf(upper,m1,std1)-norm.cdf(r2,m1,std1)) 
     else:     # order: l-r1-u-r2 
      overlap = (norm.cdf(r1,m1,std1)-norm.cdf(lower,m1,std1))+(norm.cdf(upper,m2,std2)-norm.cdf(r1,m2,std2)) 
    else:      # l-u-r1-r2 
     overlap = (norm.cdf(upper,m1,std1)-norm.cdf(lower,m1,std1))