CODEWITHSUNDEEP
numpymatplotlibscipy
user3398600
user3398600

Reputation: 805

sympy Maximum-Likelihood: clarification of script from book

the code below is from the book "Python for Probability, Statistics, and Machine Learning. The clarification needed is for the plotting section. The problem is that "logJ" in the script is not-defined. However, the book provides this as the code to plot the graph. How do you correct (code) the plotting part of the script so it plots the output shown?

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

Expected output:

enter image description here

Upvotes: 0

Views: 104

Answers (1)

Jayvee
Jayvee

Reputation: 10875

With a couple of changes (logL to logJ, and map made into list) displays the graph:

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

enter image description here

Also, I included %matplotlib as I tested in a Jupyter notebook

Upvotes: 1

Related Questions