Bernoulli

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
from scipy.optimize import minimize

The log-likelihood for the Bernoulli distributon is

\[\ell(p | \mathbf{X}) = \sum_{n=1}^N X_n \log{p} + (1 - X_n) \log{1-p}\]

To find the maximum likelihood estimator for \(p\), we first find \(\frac{d\ell}{dp}\):

\[\frac{d\ell}{dp} = \sum_{n=1}^N X_n / p - (1 - X_n) / (1 - p)\]

Then set \(\frac{d\ell}{dp}\) equal to \(0\)

\[(1 - p) \sum_{n=1}^N X_n = p \sum_{n=1}^N (1 - X_n)\]

Last, solve for \(p\)

\[\hat{p} = N^{-1}\sum_{n=1}^N X_n\]

def ll_bernoulli(p, x):
    N = np.size(x)
    sx = np.sum(x)
    return -sx * np.log(p) + (sx - N) * np.log1p(-p)

rng = np.random.default_rng()
data = rng.binomial(1, 0.3, size = 1001)
minimize(ll_bernoulli, rng.uniform(), args = (data,), 
         method = "L-BFGS-B", bounds = [(1e-5, 1 - 1e-5)])

  message: CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH
  success: True
   status: 0
      fun: 624.1970903914475
        x: [ 3.157e-01]
      nit: 4
      jac: [ 6.821e-05]
     nfev: 12
     njev: 6
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>

np.mean(data)

0.3156843156843157

Exponential

The log-likelihood for the Exponential distributon is

\[\ell( \lambda | \mathbf{X}) = N \log{\lambda} - \lambda \sum_{n=1}^N X_n\]

To find the maximum likelihood estimator for \(\lambda\), we first find \(\frac{d\ell}{d\lambda}\):

\[\frac{d\ell}{d\lambda} = N / \lambda - \sum_{n=1}^N X_n\]

Then set \(\frac{d\ell}{d\lambda}\) equal to \(0\)

\[N / \lambda = \sum_{n=1}^N X_n\]

Last, solve for \(\lambda\)

\[\hat{\lambda} = N / \sum_{n=1}^N X_n\]

def ll_exponential(l, x):
    N = np.size(x)
    return -N * np.log(l) + l * np.sum(x)

data = rng.exponential(0.2, size = 1001) # answer should be close to 5
minimize(ll_exponential, rng.uniform(), args = (data,), 
         method = "L-BFGS-B", bounds = [(1e-5, np.inf)])

  message: CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH
  success: True
   status: 0
      fun: -601.80496199843
        x: [ 4.959e+00]
      nit: 9
      jac: [-1.137e-05]
     nfev: 20
     njev: 10
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>

1 / np.mean(data)

4.958998267986209

Other Formats

Exponential