Maximum Likelihood function using Sympy returns empty list

Question

I wanted to create a function that would return estimator calculated by Maximum Likelihood Function. The function I made is below:

def Maximum_Likelihood(param, pmf):
    
    i = symbols('i', positive=True)
    n = symbols('n', positive=True)
    
    Likelihood_function = Product(pmf, (i, 1, n)) 
    # calculate partial derivative for parameter (p for Bernoulli)
    deriv = diff(Likelihood_function, param) 
    equation_to_solve = Eq(deriv,0) # equate with 0
    
    # solve above equation and return parameter (p for Bernoulli)
    return solve(equation_to_solve, param) 

Param means parameter for which I want to know estimator and pmf is a probability mass function.

And for example, I want to get an estimator for parameter p in Bernoulli distribution. How the Maximum Likelihood should look is:

My code. Imports:

import numpy as np
import sympy as sym

from sympy.solvers import solve

from sympy import Product, Function, oo, IndexedBase, diff, Eq, symbols

Now, using Sympy I defined it:

def Maximum_Likelihood(param, pmf):
    
    i = symbols('i', positive=True)
    n = symbols('n', positive=True)
    
    Likelihood_function = Product(pmf, (i, 1, n)) 
    deriv = diff(Likelihood_function, param) 
    equation_to_solve = Eq(deriv,0) 
    
    return solve(equation_to_solve, param) 

and Bernoulli example:

x = IndexedBase('x')
i = symbols('i', positive=True)
n = symbols('n', positive=True)
formula = (p**x[i])*((1-p)**(1-x[i]))

Likelihood_function = Product(formula, (i, 1, n))
Likelihood_function

When I want to get the outcome of Maximum_Likelihood(param, pmf):

param = p 
pmf = formula
print(Maximum_Likelihood(param, pmf))

I get "[]". I want to obtain estimator of p that should look like:

Could you please take a look at it and advice what I do wrong. Thank you!

Answer 1

For some reason the diff of a product doesn't actually evaluate the derivative but you can use doit to force evaluation:

In [14]: solve(Eq(diff(Product(p**x[i]*(1 - p)**(1 - x[i]), (i, 1, n)), p), 0).doit(), p)
Out[14]: 
⎡  n       ⎤
⎢ ___      ⎥
⎢ ╲        ⎥
⎢  ╲       ⎥
⎢  ╱   x[i]⎥
⎢ ╱        ⎥
⎢ ‾‾‾      ⎥
⎢i = 1     ⎥
⎢──────────⎥
⎣    n     ⎦

That says that the MLE estimate for p is just the sample mean of the data I guess.