高斯相交

Question

我只是想画出两个高斯人的相交点。我有下面的代码。但它并没有画出确切的交叉口，我真的不知道为什么。这就像刚刚稍微偏离，但我工作了导出的解决方案，如果我们采取的对数减去高斯，是的，它似乎应该是正确的。有人能帮忙吗？非常感谢!

import numpy as np 
import matplotlib.pyplot as plt 

def plot_normal(x, mean = 0, sigma = 1):
    return 1.0/(2*np.pi*sigma**2) * np.exp(-((x-mean)**2)/(2*sigma**2))

# found online
def solve_gasussians(m1, s1, m2, s2):
  a = 1.0/(2.0*s1**2) - 1.0/(2.0*s2**2)
  b = m2/(s2**2) - m1/(s1**2)
  c = m1**2 /(2*s1**2) - m2**2 / (2.0*s2**2) - np.log(s2/s1)
  return np.roots([a,b,c])

s1 = np.linspace(0, 10,300)
s2 = np.linspace(0, 14, 300)

solved_val = solve_gasussians(5.0, 0.5, 7.0, 1.0)
print solved_val
solved_val = solved_val[0]
plt.figure('Baseline Distributions')
plt.title('Baseline Distributions')
plt.xlabel('Response Rate')
plt.ylabel('Probability')
plt.plot(s1, plot_normal(s1, 5.0, 0.5),'r', label='s1')
plt.plot(s2, plot_normal(s2, 7.0, 1.0),'b', label='s2')
plt.plot(solved_val, plot_normal(solved_val, 7.0, 1.0), 'mo')
plt.legend()
plt.show()

Answer 1

plot_normal函数中有一个小错误--分母中缺少平方根。适当版本：

def plot_normal(x, mean = 0, sigma = 1):
    return 1.0/np.sqrt(2*np.pi*sigma**2) * np.exp(-((x-mean)**2)/(2*sigma**2))

给出了预期结果：

​

还有两句话。

1. 请记住，一般情况下，您可以有两个方程的根(两个交点)，这就是您提供的参数的情况。
2. 据我所知，np.roots给出了近似结果，但是cat很容易得到准确的结果，将solve_gasussians函数重写为： solve_gasussians(m1，s1，m2，( s2)：二次方程的#系数ax^2 + bx +c=0 a= ( s1**2.0 ) - ( s2**2.0 ) b=2* (m1 * s2**2.0 - m2 * s1**2.0) c= m2**2.0 *S1*2.0-M1*2.0*S2*2.0-2*S1*2.0*S2*2.0* np.log(s1 )/s2) x1 = (-b +np.sqrt(b*2.0- 4.0 *a*c)/ (2.0 * a) x2 = (-b -np.sqrt(b**-2.0- 4.0 *a* c)) / (2.0 * a)返回x1，x2

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 = 5
std1 = 0.5
m2 = 7
std2 = 1

result = solve(m1,m2,std1,std2)

x = np.linspace(-5,9,10000)
plot1=plt.plot(x,[norm.pdf(_,m1,std1) for _ in x])
plot2=plt.plot(x,[norm.pdf(_,m2,std2) for _ in x])
plot3=plt.plot(result[0],norm.pdf(result[0],m1,std1) ,'o')

plt.show()

我会提供两条不请自来的建议，让你过得更轻松(就像他们对我做的那样)：

• 当您调整代码时，请尝试进行小的增量更改，并检查代码在每个步骤中是否仍然工作。
• 寻找现有的免费库。在本例中，来自规范的是原始代码中使用的很好的替代品。

Answer 3

错误就在这里。这一行：

def plot_normal(x, mean = 0, sigma = 1):
  return 1.0/(2*np.pi*sigma**2) * np.exp(-((x-mean)**2)/(2*sigma**2))

应该是这样：

def plot_normal(x, mean = 0, sigma = 1):
  return 1.0/np.sqrt(2*np.pi*sigma**2) * np.exp(-((x-mean)**2)/(2*sigma**2))

你忘了sqrt。

如果可用的话，使用预先存在的普通pdf会更明智，例如：

import scipy.stats
def plot_normal(x, mean = 0, sigma = 1):
  return scipy.stats.norm.pdf(x,loc=mean,scale=sigma)

也有可能精确地求解这些交叉口。This answer给出了高斯交叉口根的二次方程。使用maxima求解x给出了以下表达式。它虽然复杂，但不依赖于迭代方法，可以从更简单的表达式自动生成。

def solve_gaussians(m1,s1,m2,s2):
  x1 = (s1*s2*np.sqrt((-2*np.log(s1/s2)*s2**2)+2*s1**2*np.log(s1/s2)+m2**2-2*m1*m2+m1**2)+m1*s2**2-m2*s1**2)/(s2**2-s1**2)
  x2 = -(s1*s2*np.sqrt((-2*np.log(s1/s2)*s2**2)+2*s1**2*np.log(s1/s2)+m2**2-2*m1*m2+m1**2)-m1*s2**2+m2*s1**2)/(s2**2-s1**2)
  return x1,x2

总之，它提供了：

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

def plot_normal(x, mean = 0, sigma = 1):
  return scipy.stats.norm.pdf(x,loc=mean,scale=sigma)

#Use the equation from [this answer](https://stats.stackexchange.com/a/12213/12116) solved for x
def solve_gaussians(m1,s1,m2,s2):
  x1 = (s1*s2*np.sqrt((-2*np.log(s1/s2)*s2**2)+2*s1**2*np.log(s1/s2)+m2**2-2*m1*m2+m1**2)+m1*s2**2-m2*s1**2)/(s2**2-s1**2)
  x2 = -(s1*s2*np.sqrt((-2*np.log(s1/s2)*s2**2)+2*s1**2*np.log(s1/s2)+m2**2-2*m1*m2+m1**2)-m1*s2**2+m2*s1**2)/(s2**2-s1**2)
  return x1,x2

s = np.linspace(0, 14,300)
x = solve_gaussians(5.0,0.5,7.0,1.0)

plt.figure('Baseline Distributions')
plt.title('Baseline Distributions')
plt.xlabel('Response Rate')
plt.ylabel('Probability')
plt.plot(s, plot_normal(s, 5.0, 0.5),'r', label='s1')
plt.plot(s, plot_normal(s, 7.0, 1.0),'b', label='s2')
plt.plot(x[0],plot_normal(x[0],5.,0.5),'mo')
plt.plot(x[1],plot_normal(x[1],5.,0.5),'mo')
plt.legend()
plt.show()

给予：

​

​