Atlas/Statistics Inference/Series/Series6.ipynb

Import Required Libraries

In [21]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb

Exercise 1: Binomial Distribution MLE

Problem Statement

Consider the binomial distribution of the number k of successes in n independent trials:

\sum_{i=1}^{n} X_i \sim \text{Binomial}(n, p)

The log-likelihood function is:

\log L(p) = \log \binom{n}{k} + k \log(p) + (n-k) \log(1-p)

Part 1: Plot the Log-Likelihood Function

Given: n = 100, k = 10

In [22]:

# Parameters
n = 100
k = 10

In [23]:

# Create a sequence of p values
p_values = np.linspace(0.01, 0.99, 100)

# Define the log-likelihood function
def log_likelihood(p, n, k):
    return np.log(comb(n, k)) + k * np.log(p) + (n - k) * np.log(1 - p)

# Calculate log-likelihood for each p
log_lik_values = [log_likelihood(p, n, k) for p in p_values]

# Update parameters for MLE calculation
k = 15
n = 100

# Calculate MLE
p_MLE = k / n
print(f"MLE of p: {p_MLE}")

MLE of p: 0.15

In [24]:

# Plot the log-likelihood function
plt.figure(figsize=(10, 6))
plt.plot(p_values, log_lik_values, 'b-', linewidth=2)
plt.xlabel('p', fontsize=12)
plt.ylabel('Log-Likelihood', fontsize=12)
plt.title('Log-Likelihood Function for Binomial(100, p) with k=10', fontsize=14)
plt.grid(True, alpha=0.3)
plt.show()

In [25]:

# Summary results
print("=" * 50)
print("Summary of Results")
print("=" * 50)
print(f"Sample size (n):              {n}")
print(f"Observed successes (k):       {k}")
print(f"MLE of p:                     {p_MLE}")
print(f"Comparison p:                 0.10")
print(f"\nLog-likelihood at p = {p_MLE}:   {log_likelihood(p_MLE, n, k):.4f}")
print(f"Log-likelihood at p = 0.10:   {log_likelihood(0.10, n, k):.4f}")
print("=" * 50)
print(f"\nThe MLE p̂ = {p_MLE} maximizes the log-likelihood function")
print("This is simply the sample proportion of successes")

==================================================
Summary of Results
==================================================
Sample size (n):              100
Observed successes (k):       15
MLE of p:                     0.15
Comparison p:                 0.10

Log-likelihood at p = 0.15:   -2.1974
Log-likelihood at p = 0.10:   -3.4209
==================================================

The MLE p̂ = 0.15 maximizes the log-likelihood function
This is simply the sample proportion of successes

Summary of Results

Let's create a summary table of our findings:

In [26]:

# Recalculate log-likelihood values with updated k
p_values = np.linspace(0.01, 0.99, 100)
log_lik_values = [log_likelihood(p, n, k) for p in p_values]

# Create the plot
plt.figure(figsize=(12, 7))
plt.plot(p_values, log_lik_values, 'b-', linewidth=2, label='Log-Likelihood')

# Add vertical lines for p = 0.10 and p = p_MLE
plt.axvline(x=0.10, color='red', linewidth=2, linestyle='--', label='p = 0.10')
plt.axvline(x=p_MLE, color='green', linewidth=2, linestyle='--', label=f'p = {p_MLE} (MLE)')

# Add points at specific p values
plt.plot(0.10, log_likelihood(0.10, n, k), 'ro', markersize=10)
plt.plot(p_MLE, log_likelihood(p_MLE, n, k), 'go', markersize=10)

# Labels and formatting
plt.xlabel('p', fontsize=12)
plt.ylabel('Log-Likelihood', fontsize=12)
plt.title('Log-Likelihood with Binomial Distributions', fontsize=14)
plt.legend(loc='upper right', fontsize=10)
plt.grid(True, alpha=0.3)
plt.show()

Part 3: Compare Log-Likelihood at Different p Values

Let's visualize the log-likelihood function with the MLE and compare it to p = 0.10.

In [27]:

# Update parameters for MLE calculation
k = 15
n = 100

# Calculate MLE
p_MLE = k / n
print(f"MLE of p: {p_MLE}")

MLE of p: 0.15

Part 2: Compute the MLE

Given: \sum_{i=1}^{n} x_i = 15 with n = 100 trials

Finding the MLE

To find the maximum likelihood estimator, we take the derivative of the log-likelihood with respect to p and set it equal to zero:

\frac{d}{dp} \log L(p) = \frac{k}{p} - \frac{n-k}{1-p} = 0

Solving for p:

\frac{k}{p} = \frac{n-k}{1-p} k(1-p) = p(n-k) k - kp = pn - pk k = pn \hat{p}_{MLE} = \frac{k}{n}

In [28]:

# Create a sequence of p values
p_values = np.linspace(0.01, 0.99, 100)

# Calculate log-likelihood for each p
log_lik_values = [log_likelihood(p, n, k) for p in p_values]

In [29]:

# Define the log-likelihood function
def log_likelihood(p, n, k):
    return np.log(comb(n, k)) + k * np.log(p) + (n - k) * np.log(1 - p)

In [30]:

# Define the log-likelihood function
def log_likelihood(p, n, k):
    return np.log(comb(n, k)) + k * np.log(p) + (n - k) * np.log(1 - p)

In [31]:

print(f"MLE of p: {p_MLE}")

print(f"Log-likelihood at p = {p_MLE}:   {log_likelihood(p_MLE, n, k):.4f}")

MLE of p: 0.15
Log-likelihood at p = 0.15:   -2.1974