交感极大似然:从书中澄清剧本

Question

下面的代码来自于“概率、统计和机器学习的Pythonfor概率论、统计和机器学习”(Pythonfor概率论、统计和机器学习)一书。需要澄清的是绘图部分。问题是没有定义脚本中的"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, p, 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, logJ.subs(p, p_true), 's', color = 'k', ms=15, 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)

预期产出：

​

​

Answer 1

通过一些更改(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, p, 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, logJ.subs(p, p_true), 's', color = 'k', ms=15, 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)

​

​

此外，我还在木星笔记本上测试了%matplotlib。