Find parameters with maximum likelihood method from histograms

Question

For a problem I am asked to generate two random exponential distributions with different lambda's, mix them up, and then extract two lambdas with the method of moments. In the second part of the problem I am asked to find the two lambda's from the mixed dataset for each histogram using the maximum likelihood method and finding probabilities by integrating over the probability density function for each bin, then plot the maximum likelihood function near the obtained maxima. I am kinda stuck on how to get the two lambdas from each histogram using the maximum likelihood function. Below is my code:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from random import choice
from scipy.optimize import fsolve
#Defining initial conditions for the random exponential distribution
n=100
lambda1=10
lambda2=50
#Creating two datasets with different lambdas
dataset1=pd.DataFrame(np.random.default_rng().exponential(scale=(1/lambda1),size=n))
dataset2=pd.DataFrame(np.random.default_rng().exponential(scale=(1/lambda2), size=n))
#Mixing all the data
mixeddata=pd.concat([dataset1,dataset2])
#Creating histograms for the different amount of bins
bins=[5,10,25,50]
fig, (ax1,ax2,ax3,ax4)= plt.subplots(4)
ax1.hist(mixeddata, bins=5,edgecolor='black', density=True)
ax2.hist(mixeddata, bins=10,edgecolor='black', density=True)
ax3.hist(mixeddata, bins=25,edgecolor='black', density=True)
ax4.hist(mixeddata, bins=50,edgecolor='black', density=True)
ax1.set_ylim(0,1)
ax2.set_ylim(0,1)
ax3.set_ylim(0,1)
ax4.set_ylim(0,1)

plt.show()
#Using method of moments to find lambda1 and lambda2 from the mixed dataset
sample_mean = mixeddata.mean().values[0]
sample_variance = mixeddata.var().values[0]

def equations(params):
    lambda1, lambda2, p = params
    eq1 = p * (1/lambda1) + (1-p) * (1/lambda2) - sample_mean
    eq2 = p * (1/(lambda1**2)) + (1-p) * (1/(lambda2**2)) + p*(1-p) * ((1/lambda1) - (1/lambda2))**2 - sample_variance
    # Assuming equal mixing probabilities for simplicity
    eq3 = p - 0.5
    return [eq1, eq2, eq3]

# Initial guesses for lambda1, lambda2, and p
initial_guesses = [10, 50, 0.5]
solutions = fsolve(equations, initial_guesses)

lambda1_est, lambda2_est, p_est = solutions
print(f'Estimated Parameters: lambda1 = {lambda1_est}, lambda2 = {lambda2_est}, p = {p_est}')

I tried setting up a function but I get no results.

Answer 1

Don't use Pandas for this application and don't vertically truncate your graphs. Otherwise, you already have code that generates lambdas:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

# Defining initial conditions for the random exponential distribution
n = 1000
lambda1 = 10
lambda2 = 50

# Creating two datasets with different lambdas
dataset1 = np.random.default_rng().exponential(scale=1/lambda1,size=n)
dataset2 = np.random.default_rng().exponential(scale=1/lambda2, size=n)

# Mixing all the data
mixeddata = np.concatenate((dataset1, dataset2))

# Creating histograms for the different amount of bins
fig, axes = plt.subplots(4)
for (ax, bins) in zip(axes, (5, 10, 25, 50)):
    ax.hist(mixeddata, bins=bins, edgecolor='black', density=True)
plt.show()

# Using method of moments to find lambda1 and lambda2 from the mixed dataset
sample_mean = mixeddata.mean()
sample_variance = mixeddata.var()


def equations(params: np.ndarray) -> list[float]:
    lambda1, lambda2, p = params
    eq1 = p * (1/lambda1) + (1-p) * (1/lambda2) - sample_mean
    eq2 = p * (1/(lambda1**2)) + (1-p) * (1/(lambda2**2)) + p*(1-p) * ((1/lambda1) - (1/lambda2))**2 - sample_variance
    # Assuming equal mixing probabilities for simplicity
    eq3 = p - 0.5
    return [eq1, eq2, eq3]


# Initial guesses for lambda1, lambda2, and p
initial_guesses = (10, 50, 0.5)
solutions = fsolve(func=equations, x0=initial_guesses)

lambda1_est, lambda2_est, p_est = solutions
print(f'Estimated Parameters: lambda1={lambda1_est:.2f}, lambda2={lambda2_est:.2f}, p={p_est:.2f}')

Estimated Parameters: lambda1=10.27, lambda2=51.32, p=0.50