sympy 最大似然：从书中澄清脚本

如何解决sympy 最大似然：从书中澄清脚本

下面的代码来自“Python for Probability,Statistics,and Machine Learning”一书。需要澄清的是绘图部分。问题是脚本中的“logJ”没有定义。但是，这本书提供了这是绘制图形的代码。您如何更正（编码）脚本的绘图部分，以便绘制显示的输出？

    # coin-flipping experiment using a Bernoulli distribution
    from scipy.stats import bernoulli
    import matplotlib.pyplot as plt
    import sympy
    import numpy as np

    #Simulating
    p_true = 1/2.0
    fp = bernoulli(p_true)
    xs = fp.rvs(100)
    print(xs[:30])


    x,p,z = sympy.symbols('x,z',positive = True)
    phi = P**x*(1-p)**(1-x)
    L = np.prod([phi.subs(x,i) for i in xs])
    print(L)

    logL = sympy.expand_log(sympy.log(L))
    sol,= sympy.solve(sympy.diff(logL,p),p)
    print(sol)

    #This is the plotting part,which does not work because LogJ is not defined
    fig,ax = plt.subplots()
    x = np.linspace(0,1,100)
    ax.plot(x,map(sympy.lambdify(p,logJ,'numpy'),x),'k-',lw=3)
    ax.plot(sol,logJ.subs(p,sol),'o',color = 'gray',ms=15,label = 'Estimated')
    ax.plot(p_true,p_true),'s',color = 'k',label = 'Actual')
    ax.set_xlabel('$p$',fontsize = 18)
    ax.set_ylabel('Likelihood',fontsize = 18)
    ax.set_title('Estimate not equal to true value',fontsize = 18)
    ax.legend(loc=0)

预期输出：

解决方法

通过一些更改（logL 到 logJ，并将地图制作成列表）显示图表：

# coin-flipping experiment using a Bernoulli distribution
from scipy.stats import bernoulli
import matplotlib.pyplot as plt
import sympy
import numpy as np
%matplotlib
#Simulating
p_true = 1/2.0
fp = bernoulli(p_true)
xs = fp.rvs(100)
print(xs[:30])


x,p,z = sympy.symbols('x,z',positive = True)
phi = p**x*(1-p)**(1-x)
L = np.prod([phi.subs(x,i) for i in xs])
print(L)

logJ = sympy.expand_log(sympy.log(L))
sol,= sympy.solve(sympy.diff(logJ,p),p)
print(sol)

#This is the plotting part,which does not work because LogJ is not defined
fig,ax = plt.subplots()
x = np.linspace(0,1,100)
ax.plot(x,list(map(sympy.lambdify(p,logJ,'numpy'),x)),'k-',lw=3)
ax.plot(sol,logJ.subs(p,sol),'o',color = 'gray',ms=15,label = 'Estimated')
ax.plot(p_true,p_true),'s',color = 'k',label = 'Actual')
ax.set_xlabel('$p$',fontsize = 18)
ax.set_ylabel('Likelihood',fontsize = 18)
ax.set_title('Estimate not equal to true value',fontsize = 18)
ax.legend(loc=0)

此外，我在 Jupyter 笔记本中测试时包含了 %matplotlib