marimo-team/marimo-learn
3
1# /// script2# requires-python = ">=3.10"3# dependencies = [4# "marimo",5# "matplotlib==3.10.8",6# "scipy==1.17.1",7# "numpy==2.4.3",8# "polars==1.24.0",9# "plotly==5.18.0",10# ]11# ///12 13import marimo14 15__generated_with = "0.18.4"16app = marimo.App(width="medium", app_title="Maximum Likelihood Estimation")17 18 19@app.cell(hide_code=True)20def _(mo):21 mo.md(r"""22 # Maximum Likelihood Estimation23 24 _This notebook is a computational companion to ["Probability for Computer Scientists"](https://chrispiech.github.io/probabilityForComputerScientists/en/part5/mle/), by Stanford professor Chris Piech._25 26 Maximum Likelihood Estimation (MLE) is a fundamental method in statistics for estimating parameters of a probability distribution. The central idea is elegantly simple: **choose the parameters that make the observed data most likely**.27 28 In this notebook, we'll try to understand MLE, starting with the core concept of likelihood and how it differs from probability. We'll explore how to formulate MLE problems mathematically and then solve them for various common distributions. Along the way, I've included some interactive visualizations to help build your intuition about these concepts. You'll see how MLE applies to real-world scenarios like linear regression, and hopefully gain a deeper appreciation for why this technique is so widely used in statistics and machine learning. Think of MLE as detective work - we have some evidence (our data) and we're trying to figure out the most plausible explanation (our parameters) for what we've observed.29 """)30 return31 32 33@app.cell(hide_code=True)34def _(mo):35 mo.md(r"""36 ## Likelihood: The Core Concept37 38 Before diving into MLE, we need to understand what "likelihood" means in a statistical context.39 40 ### Data and Parameters41 42 Suppose we have collected some data $X_1, X_2, \ldots, X_n$ that are independent and identically distributed (IID). We assume these data points come from a specific type of distribution (like Normal, Bernoulli, etc.) with unknown parameters $\theta$.43 44 ### What is Likelihood?45 46 Likelihood measures how probable our observed data is, given specific values of the parameters $\theta$.47 48 /// note49 **Probability vs. Likelihood**50 51 - **Probability**: Given parameters $\theta$, what's the chance of observing data $X$?52 - **Likelihood**: Given observed data $X$, how likely are different parameter values $\theta$?53 ///54 55 To simplify notation, we'll use $f(X=x|\Theta=\theta)$ to represent either the PMF or PDF of our data, conditioned on the parameters.56 """)57 return58 59 60@app.cell(hide_code=True)61def _(mo):62 mo.md(r"""63 ### The Likelihood Function64 65 Since we assume our data points are independent, the likelihood of all our data is the product of the likelihoods of each individual data point:66 67 $$L(\theta) = \prod_{i=1}^n f(X_i = x_i|\Theta = \theta)$$68 69 This function $L(\theta)$ gives us the likelihood of observing our entire dataset for different parameter values $\theta$.70 71 /// tip72 **Key Insight**: Different parameter values produce different likelihoods for the same data. Better parameter values will make the observed data more likely.73 ///74 """)75 return76 77 78@app.cell(hide_code=True)79def _(mo):80 mo.md(r"""81 ## Maximum Likelihood Estimation82 83 The core idea of MLE is to find the parameter values $\hat{\theta}$ that maximize the likelihood function:84 85 $$\hat{\theta} = \underset{\theta}{\operatorname{argmax}} \, L(\theta)$$86 87 The notation $\hat{\theta}$ represents our best estimate of the true parameters based on the observed data.88 89 ### Working with Log-Likelihood90 91 In practice, we usually work with the **log-likelihood** instead of the likelihood directly. Since logarithm is a monotonically increasing function, the maximum of $L(\theta)$ occurs at the same value of $\theta$ as the maximum of $\log L(\theta)$.92 93 Taking the logarithm transforms our product into a sum, which is much easier to work with:94 95 $$LL(\theta) = \log L(\theta) = \log \prod_{i=1}^n f(X_i=x_i|\Theta = \theta) = \sum_{i=1}^n \log f(X_i = x_i|\Theta = \theta)$$96 97 /// warning98 Working with products of many small probabilities can lead to numerical underflow. Taking the logarithm converts these products to sums, which is numerically more stable.99 ///100 """)101 return102 103 104@app.cell(hide_code=True)105def _(mo):106 mo.md(r"""107 ### Finding the Maximum108 109 To find the values of $\theta$ that maximize the log-likelihood, we typically:110 111 1. Take the derivative of $LL(\theta)$ with respect to each parameter112 2. Set each derivative equal to zero113 3. Solve for the parameters114 115 Let's see this approach in action with some common distributions.116 """)117 return118 119 120@app.cell(hide_code=True)121def _(mo):122 mo.md(r"""123 ## MLE for Bernoulli Distribution124 125 > _Note:_ The following derivation is included as reference material. The credit for this mathematical formulation belongs to ["Probability for Computer Scientists"](https://chrispiech.github.io/probabilityForComputerScientists/en/part5/mle/) by Chris Piech.126 127 Let's start with a simple example: estimating the parameter $p$ of a Bernoulli distribution.128 129 ### The Model130 131 A Bernoulli distribution has a single parameter $p$ which represents the probability of success (getting a value of 1). Its probability mass function (PMF) can be written as:132 133 $$f(x|p) = p^x(1-p)^{1-x}, \quad x \in \{0, 1\}$$134 135 This elegant formula works because:136 137 - When $x = 1$: $f(1|p) = p^1(1-p)^0 = p$138 - When $x = 0$: $f(0|p) = p^0(1-p)^1 = 1-p$139 140 ### Deriving the MLE141 142 Given $n$ independent Bernoulli trials $X_1, X_2, \ldots, X_n$, we want to find the value of $p$ that maximizes the likelihood of our observed data.143 144 Step 1: Write the likelihood function145 $$L(p) = \prod_{i=1}^n p^{x_i}(1-p)^{1-x_i}$$146 147 Step 2: Take the logarithm to get the log-likelihood148 $$\begin{align*}149 LL(p) &= \sum_{i=1}^n \log(p^{x_i}(1-p)^{1-x_i}) \\150 &= \sum_{i=1}^n \left[x_i \log(p) + (1-x_i)\log(1-p)\right] \\151 &= \left(\sum_{i=1}^n x_i\right) \log(p) + \left(n - \sum_{i=1}^n x_i\right) \log(1-p) \\152 &= Y\log(p) + (n-Y)\log(1-p)153 \end{align*}$$154 155 where $Y = \sum_{i=1}^n x_i$ is the total number of successes.156 157 Step 3: Find the value of $p$ that maximizes $LL(p)$ by setting the derivative to zero158 $$\begin{align*}159 \frac{d\,LL(p)}{dp} &= \frac{Y}{p} - \frac{n-Y}{1-p} = 0 \\160 \frac{Y}{p} &= \frac{n-Y}{1-p} \\161 Y(1-p) &= p(n-Y) \\162 Y - Yp &= pn - pY \\163 Y &= pn \\164 \hat{p} &= \frac{Y}{n} = \frac{\sum_{i=1}^n x_i}{n}165 \end{align*}$$166 167 /// tip168 The MLE for the parameter $p$ in a Bernoulli distribution is simply the **sample mean** - the proportion of successes in our data!169 ///170 """)171 return172 173 174@app.cell(hide_code=True)175def _(controls):176 controls.center()177 return178 179 180@app.cell(hide_code=True)181def _(generate_button, mo, np, plt, sample_size_slider, true_p_slider):182 # generate bernoulli samples when button is clicked183 bernoulli_button_value = generate_button.value184 185 # get parameter values186 bernoulli_true_p = true_p_slider.value187 bernoulli_n = sample_size_slider.value188 189 # generate data190 bernoulli_data = np.random.binomial(1, bernoulli_true_p, size=bernoulli_n)191 bernoulli_Y = np.sum(bernoulli_data)192 bernoulli_p_hat = bernoulli_Y / bernoulli_n193 194 # create visualization195 bernoulli_fig, (bernoulli_ax1, bernoulli_ax2) = plt.subplots(1, 2, figsize=(12, 5))196 197 # plot data histogram198 bernoulli_ax1.hist(bernoulli_data, bins=[-0.5, 0.5, 1.5], rwidth=0.8, color='lightblue')199 bernoulli_ax1.set_xticks([0, 1])200 bernoulli_ax1.set_xticklabels(['Failure (0)', 'Success (1)'])201 bernoulli_ax1.set_title(f'Bernoulli Data: {bernoulli_n} samples')202 bernoulli_ax1.set_ylabel('Count')203 bernoulli_y_counts = [bernoulli_n - bernoulli_Y, bernoulli_Y]204 for bernoulli_idx, bernoulli_count in enumerate(bernoulli_y_counts):205 bernoulli_ax1.text(bernoulli_idx, bernoulli_count/2, f"{bernoulli_count}", 206 ha='center', va='center', 207 color='white' if bernoulli_idx == 0 else 'black', 208 fontweight='bold')209 210 # calculate log-likelihood function211 bernoulli_p_values = np.linspace(0.01, 0.99, 100)212 bernoulli_ll_values = np.zeros_like(bernoulli_p_values)213 214 for bernoulli_i, bernoulli_p in enumerate(bernoulli_p_values):215 bernoulli_ll_values[bernoulli_i] = bernoulli_Y * np.log(bernoulli_p) + (bernoulli_n - bernoulli_Y) * np.log(1 - bernoulli_p)216 217 # plot log-likelihood218 bernoulli_ax2.plot(bernoulli_p_values, bernoulli_ll_values, 'b-', linewidth=2)219 bernoulli_ax2.axvline(x=bernoulli_p_hat, color='r', linestyle='--', label=f'MLE: $\\hat{{p}} = {bernoulli_p_hat:.3f}$')220 bernoulli_ax2.axvline(x=bernoulli_true_p, color='g', linestyle='--', label=f'True: $p = {bernoulli_true_p:.3f}$')221 bernoulli_ax2.set_xlabel('$p$ (probability of success)')222 bernoulli_ax2.set_ylabel('Log-Likelihood')223 bernoulli_ax2.set_title('Log-Likelihood Function')224 bernoulli_ax2.legend()225 226 plt.tight_layout()227 plt.gca()228 229 # Create markdown to explain the results230 bernoulli_explanation = mo.md(231 f"""232 ### Bernoulli MLE Results233 234 **True parameter**: $p = {bernoulli_true_p:.3f}$ 235 **Sample statistics**: {bernoulli_Y} successes out of {bernoulli_n} trials 236 **MLE estimate**: $\\hat{{p}} = \\frac{{{bernoulli_Y}}}{{{bernoulli_n}}} = {bernoulli_p_hat:.3f}$237 238 The plot on the right shows the log-likelihood function $LL(p) = Y\\log(p) + (n-Y)\\log(1-p)$. 239 The red dashed line marks the maximum likelihood estimate $\\hat{{p}}$, and the green dashed line 240 shows the true parameter value.241 242 /// note243 Try increasing the sample size to see how the MLE estimate gets closer to the true parameter value!244 ///245 """246 )247 248 # Display plot and explanation together249 mo.vstack([250 bernoulli_fig,251 bernoulli_explanation252 ])253 return254 255 256@app.cell(hide_code=True)257def _(mo):258 mo.md(r"""259 ## MLE for Normal Distribution260 261 Next, let's look at a more complex example: estimating the parameters $\mu$ and $\sigma^2$ of a Normal distribution.262 263 ### The Model264 265 A Normal (Gaussian) distribution has two parameters:266 - $\mu$: the mean267 - $\sigma^2$: the variance268 269 Its probability density function (PDF) is:270 271 $$f(x|\mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)$$272 273 ### Deriving the MLE274 275 Given $n$ independent samples $X_1, X_2, \ldots, X_n$ from a Normal distribution, we want to find the values of $\mu$ and $\sigma^2$ that maximize the likelihood of our observed data.276 277 Step 1: Write the likelihood function278 $$L(\mu, \sigma^2) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x_i - \mu)^2}{2\sigma^2}\right)$$279 280 Step 2: Take the logarithm to get the log-likelihood281 $$\begin{align*}282 LL(\mu, \sigma^2) &= \log\prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x_i - \mu)^2}{2\sigma^2}\right) \\283 &= \sum_{i=1}^n \log\left[\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x_i - \mu)^2}{2\sigma^2}\right)\right] \\284 &= \sum_{i=1}^n \left[-\frac{1}{2}\log(2\pi\sigma^2) - \frac{(x_i - \mu)^2}{2\sigma^2}\right] \\285 &= -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n (x_i - \mu)^2286 \end{align*}$$287 288 Step 3: Find the values of $\mu$ and $\sigma^2$ that maximize $LL(\mu, \sigma^2)$ by setting the partial derivatives to zero.289 290 For $\mu$:291 $$\begin{align*}292 \frac{\partial LL(\mu, \sigma^2)}{\partial \mu} &= \frac{1}{\sigma^2}\sum_{i=1}^n (x_i - \mu) = 0 \\293 \sum_{i=1}^n (x_i - \mu) &= 0 \\294 \sum_{i=1}^n x_i &= n\mu \\295 \hat{\mu} &= \frac{1}{n}\sum_{i=1}^n x_i296 \end{align*}$$297 298 For $\sigma^2$:299 $$\begin{align*}300 \frac{\partial LL(\mu, \sigma^2)}{\partial \sigma^2} &= -\frac{n}{2\sigma^2} + \frac{1}{2(\sigma^2)^2}\sum_{i=1}^n (x_i - \mu)^2 = 0 \\301 \frac{n}{2\sigma^2} &= \frac{1}{2(\sigma^2)^2}\sum_{i=1}^n (x_i - \mu)^2 \\302 n\sigma^2 &= \sum_{i=1}^n (x_i - \mu)^2 \\303 \hat{\sigma}^2 &= \frac{1}{n}\sum_{i=1}^n (x_i - \hat{\mu})^2304 \end{align*}$$305 306 /// tip307 The MLE for a Normal distribution gives us:308 309 - $\hat{\mu}$ = sample mean310 - $\hat{\sigma}^2$ = sample variance (using $n$ in the denominator, not $n-1$)311 ///312 """)313 return314 315 316@app.cell(hide_code=True)317def _(normal_controls):318 normal_controls.center()319 return320 321 322@app.cell(hide_code=True)323def _(324 mo,325 normal_generate_button,326 normal_sample_size_slider,327 np,328 plt,329 true_mu_slider,330 true_sigma_slider,331):332 # generate normal samples when button is clicked333 normal_button_value = normal_generate_button.value334 335 # get parameter values336 normal_true_mu = true_mu_slider.value337 normal_true_sigma = true_sigma_slider.value338 normal_true_var = normal_true_sigma**2339 normal_n = normal_sample_size_slider.value340 341 # generate random data342 normal_data = np.random.normal(normal_true_mu, normal_true_sigma, size=normal_n)343 344 # calculate mle estimates345 normal_mu_hat = np.mean(normal_data)346 normal_sigma2_hat = np.mean((normal_data - normal_mu_hat)**2) # mle variance using n347 normal_sigma_hat = np.sqrt(normal_sigma2_hat)348 349 # create visualization350 normal_fig, (normal_ax1, normal_ax2) = plt.subplots(1, 2, figsize=(12, 5))351 352 # plot histogram and density curves353 normal_bins = np.linspace(min(normal_data) - 1, max(normal_data) + 1, 30)354 normal_ax1.hist(normal_data, bins=normal_bins, density=True, alpha=0.6, color='lightblue', label='Data Histogram')355 356 # plot range for density curves357 normal_x = np.linspace(min(normal_data) - 2*normal_true_sigma, max(normal_data) + 2*normal_true_sigma, 1000)358 359 # plot true and mle densities360 normal_true_pdf = (1/(normal_true_sigma * np.sqrt(2*np.pi))) * np.exp(-0.5 * ((normal_x - normal_true_mu)/normal_true_sigma)**2)361 normal_ax1.plot(normal_x, normal_true_pdf, 'g-', linewidth=2, label=f'True: N({normal_true_mu:.2f}, {normal_true_var:.2f})')362 363 normal_mle_pdf = (1/(normal_sigma_hat * np.sqrt(2*np.pi))) * np.exp(-0.5 * ((normal_x - normal_mu_hat)/normal_sigma_hat)**2)364 normal_ax1.plot(normal_x, normal_mle_pdf, 'r--', linewidth=2, label=f'MLE: N({normal_mu_hat:.2f}, {normal_sigma2_hat:.2f})')365 366 normal_ax1.set_xlabel('x')367 normal_ax1.set_ylabel('Density')368 normal_ax1.set_title(f'Normal Distribution: {normal_n} samples')369 normal_ax1.legend()370 371 # create contour plot of log-likelihood372 normal_mu_range = np.linspace(normal_mu_hat - 2, normal_mu_hat + 2, 100)373 normal_sigma_range = np.linspace(max(0.1, normal_sigma_hat - 1), normal_sigma_hat + 1, 100)374 375 normal_mu_grid, normal_sigma_grid = np.meshgrid(normal_mu_range, normal_sigma_range)376 normal_ll_grid = np.zeros_like(normal_mu_grid)377 378 # calculate log-likelihood for each grid point379 for normal_i in range(normal_mu_grid.shape[0]):380 for normal_j in range(normal_mu_grid.shape[1]):381 normal_mu = normal_mu_grid[normal_i, normal_j]382 normal_sigma = normal_sigma_grid[normal_i, normal_j]383 normal_ll = -normal_n/2 * np.log(2*np.pi*normal_sigma**2) - np.sum((normal_data - normal_mu)**2)/(2*normal_sigma**2)384 normal_ll_grid[normal_i, normal_j] = normal_ll385 386 # plot log-likelihood contour387 normal_contour = normal_ax2.contourf(normal_mu_grid, normal_sigma_grid, normal_ll_grid, levels=50, cmap='viridis')388 normal_ax2.set_xlabel('μ (mean)')389 normal_ax2.set_ylabel('σ (standard deviation)')390 normal_ax2.set_title('Log-Likelihood Contour')391 392 # mark mle and true params393 normal_ax2.plot(normal_mu_hat, normal_sigma_hat, 'rx', markersize=10, label='MLE Estimate')394 normal_ax2.plot(normal_true_mu, normal_true_sigma, 'g*', markersize=10, label='True Parameters')395 normal_ax2.legend()396 397 plt.colorbar(normal_contour, ax=normal_ax2, label='Log-Likelihood')398 plt.tight_layout()399 plt.gca()400 401 # relevant markdown for the results402 normal_explanation = mo.md(403 rf"""404 ### Normal MLE Results405 406 **True parameters**: $\mu = {normal_true_mu:.3f}$, $\sigma^2 = {normal_true_var:.3f}$ 407 **MLE estimates**: $\hat{{\mu}} = {normal_mu_hat:.3f}$, $\hat{{\sigma}}^2 = {normal_sigma2_hat:.3f}$408 409 The left plot shows the data histogram with the true Normal distribution (green) and the MLE-estimated distribution (red dashed).410 411 The right plot shows the log-likelihood function as a contour map in the $(\mu, \sigma)$ parameter space. The maximum likelihood estimates are marked with a red X, while the true parameters are marked with a green star.412 413 /// note414 Notice how the log-likelihood contour is more stretched along the σ axis than the μ axis. This indicates that we typically estimate the mean with greater precision than the standard deviation.415 ///416 417 /// tip418 Increase the sample size to see how the MLE estimates converge to the true parameter values!419 ///420 """421 )422 423 # plot and explanation together424 mo.vstack([425 normal_fig,426 normal_explanation427 ])428 return429 430 431@app.cell(hide_code=True)432def _(mo):433 mo.md(r"""434 ## MLE for Linear Regression435 436 Now let's look at a more practical example: using MLE to derive linear regression.437 438 ### The Model439 440 Consider a model where:441 - We have pairs of observations $(X_1, Y_1), (X_2, Y_2), \ldots, (X_n, Y_n)$442 - The relationship between $X$ and $Y$ follows: $Y = \theta X + Z$443 - $Z \sim N(0, \sigma^2)$ is random noise444 - Our goal is to estimate the parameter $\theta$445 446 This means that for a given $X_i$, the conditional distribution of $Y_i$ is:447 448 $$Y_i | X_i \sim N(\theta X_i, \sigma^2)$$449 450 ### Deriving the MLE451 452 Step 1: Write the likelihood function for each data point $(X_i, Y_i)$453 $$f(Y_i | X_i, \theta) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(Y_i - \theta X_i)^2}{2\sigma^2}\right)$$454 455 Step 2: Write the likelihood for all data456 $$\begin{align*}457 L(\theta) &= \prod_{i=1}^n f(Y_i, X_i | \theta) \\458 &= \prod_{i=1}^n f(Y_i | X_i, \theta) \cdot f(X_i)459 \end{align*}$$460 461 Since $f(X_i)$ doesn't depend on $\theta$, we can simplify:462 $$L(\theta) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(Y_i - \theta X_i)^2}{2\sigma^2}\right) \cdot f(X_i)$$463 464 Step 3: Take the logarithm to get the log-likelihood465 $$\begin{align*}466 LL(\theta) &= \log \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(Y_i - \theta X_i)^2}{2\sigma^2}\right) \cdot f(X_i) \\467 &= \sum_{i=1}^n \log\left[\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(Y_i - \theta X_i)^2}{2\sigma^2}\right)\right] + \sum_{i=1}^n \log f(X_i) \\468 &= -\frac{n}{2} \log(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^n (Y_i - \theta X_i)^2 + \sum_{i=1}^n \log f(X_i)469 \end{align*}$$470 471 Step 4: Since we only care about maximizing with respect to $\theta$, we can drop terms that don't contain $\theta$:472 $$\hat{\theta} = \underset{\theta}{\operatorname{argmax}} \left[ -\frac{1}{2\sigma^2} \sum_{i=1}^n (Y_i - \theta X_i)^2 \right]$$473 474 This is equivalent to:475 $$\hat{\theta} = \underset{\theta}{\operatorname{argmin}} \sum_{i=1}^n (Y_i - \theta X_i)^2$$476 477 Step 5: Find the value of $\theta$ that minimizes the sum of squared errors by setting the derivative to zero:478 $$\begin{align*}479 \frac{d}{d\theta} \sum_{i=1}^n (Y_i - \theta X_i)^2 &= 0 \\480 \sum_{i=1}^n -2X_i(Y_i - \theta X_i) &= 0 \\481 \sum_{i=1}^n X_i Y_i - \theta X_i^2 &= 0 \\482 \sum_{i=1}^n X_i Y_i &= \theta \sum_{i=1}^n X_i^2 \\483 \hat{\theta} &= \frac{\sum_{i=1}^n X_i Y_i}{\sum_{i=1}^n X_i^2}484 \end{align*}$$485 486 /// tip487 **Key Insight**: MLE for this simple linear model gives us the least squares estimator! This is an important connection between MLE and regression.488 ///489 """)490 return491 492 493@app.cell(hide_code=True)494def _(linear_controls):495 linear_controls.center()496 return497 498 499@app.cell(hide_code=True)500def _(501 linear_generate_button,502 linear_sample_size_slider,503 mo,504 noise_sigma_slider,505 np,506 plt,507 true_theta_slider,508):509 # linear model data calc when button is clicked510 linear_button_value = linear_generate_button.value511 512 # get parameter values513 linear_true_theta = true_theta_slider.value514 linear_noise_sigma = noise_sigma_slider.value515 linear_n = linear_sample_size_slider.value516 517 # generate x data (uniformly between -3 and 3)518 linear_X = np.random.uniform(-3, 3, size=linear_n)519 520 # generate y data according to the model y = θx + z521 linear_Z = np.random.normal(0, linear_noise_sigma, size=linear_n)522 linear_Y = linear_true_theta * linear_X + linear_Z523 524 # calculate mle estimate525 linear_theta_hat = np.sum(linear_X * linear_Y) / np.sum(linear_X**2)526 527 # calculate sse for different theta values528 linear_theta_range = np.linspace(linear_true_theta - 1.5, linear_true_theta + 1.5, 100)529 linear_sse_values = np.zeros_like(linear_theta_range)530 531 for linear_i, linear_theta in enumerate(linear_theta_range):532 linear_y_pred = linear_theta * linear_X533 linear_sse_values[linear_i] = np.sum((linear_Y - linear_y_pred)**2)534 535 # convert sse to log-likelihood (ignoring constant terms)536 linear_ll_values = -linear_sse_values / (2 * linear_noise_sigma**2)537 538 # create visualization539 linear_fig, (linear_ax1, linear_ax2) = plt.subplots(1, 2, figsize=(12, 5))540 541 # plot scatter plot with regression lines542 linear_ax1.scatter(linear_X, linear_Y, color='blue', alpha=0.6, label='Data points')543 544 # plot range for regression lines545 linear_x_line = np.linspace(-3, 3, 100)546 547 # plot true and mle regression lines548 linear_ax1.plot(linear_x_line, linear_true_theta * linear_x_line, 'g-', linewidth=2, label=f'True: Y = {linear_true_theta:.2f}X')549 linear_ax1.plot(linear_x_line, linear_theta_hat * linear_x_line, 'r--', linewidth=2, label=f'MLE: Y = {linear_theta_hat:.2f}X')550 551 linear_ax1.set_xlabel('X')552 linear_ax1.set_ylabel('Y')553 linear_ax1.set_title(f'Linear Regression: {linear_n} data points')554 linear_ax1.grid(True, alpha=0.3)555 linear_ax1.legend()556 557 # plot log-likelihood function558 linear_ax2.plot(linear_theta_range, linear_ll_values, 'b-', linewidth=2)559 linear_ax2.axvline(x=linear_theta_hat, color='r', linestyle='--', label=f'MLE: θ = {linear_theta_hat:.3f}')560 linear_ax2.axvline(x=linear_true_theta, color='g', linestyle='--', label=f'True: θ = {linear_true_theta:.3f}')561 linear_ax2.set_xlabel('θ (slope parameter)')562 linear_ax2.set_ylabel('Log-Likelihood')563 linear_ax2.set_title('Log-Likelihood Function')564 linear_ax2.grid(True, alpha=0.3)565 linear_ax2.legend()566 567 plt.tight_layout()568 plt.gca()569 570 # relevant markdown to explain results571 linear_explanation = mo.md(572 f"""573 ### Linear Regression MLE Results574 575 **True parameter**: $\\theta = {linear_true_theta:.3f}$ 576 **MLE estimate**: $\\hat{{\\theta}} = {linear_theta_hat:.3f}$577 578 The left plot shows the scatter plot of data points with the true regression line (green) and the MLE-estimated regression line (red dashed).579 580 The right plot shows the log-likelihood function for different values of $\\theta$. The maximum likelihood estimate is marked with a red dashed line, and the true parameter is marked with a green dashed line.581 582 /// note583 The MLE estimate $\\hat{{\\theta}} = \\frac{{\\sum_{{i=1}}^n X_i Y_i}}{{\\sum_{{i=1}}^n X_i^2}}$ minimizes the sum of squared errors between the predicted and actual Y values.584 ///585 586 /// tip587 Try increasing the noise level to see how it affects the precision of the estimate!588 ///589 """590 )591 592 # show plot and explanation593 mo.vstack([594 linear_fig,595 linear_explanation596 ])597 return598 599 600@app.cell(hide_code=True)601def _(mo):602 mo.md(r"""603 ## Interactive Concept: Density/Mass Functions vs. Likelihood604 605 To better understand the distinction between likelihood and density/mass functions, let's create an interactive visualization. This concept is crucial for understanding why MLE works.606 """)607 return608 609 610@app.cell(hide_code=True)611def _(concept_controls):612 concept_controls.center()613 return614 615 616@app.cell(hide_code=True)617def _(concept_dist_type, mo, np, perspective_selector, plt, stats):618 # current distribution type619 concept_dist_type_value = concept_dist_type.value620 621 # view mode from dropdown622 concept_view_mode = "likelihood" if perspective_selector.value == "Likelihood Perspective" else "probability"623 624 # visualization based on distribution type625 concept_fig, concept_ax = plt.subplots(figsize=(10, 6))626 627 if concept_dist_type_value == "Normal":628 if concept_view_mode == "probability":629 # density function perspective: fixed params, varying data630 concept_mu = 0 # fixed parameter631 concept_sigma = 1 # fixed parameter632 633 # generate x values for the pdf634 concept_x = np.linspace(-4, 4, 1000)635 636 # plot pdf637 concept_pdf = stats.norm.pdf(concept_x, concept_mu, concept_sigma)638 concept_ax.plot(concept_x, concept_pdf, 'b-', linewidth=2, label='PDF: N(0, 1)')639 640 # highlight specific data values641 concept_data_points = [-2, -1, 0, 1, 2]642 concept_colors = ['#FF9999', '#FFCC99', '#99FF99', '#99CCFF', '#CC99FF']643 644 for concept_i, concept_data in enumerate(concept_data_points):645 concept_prob = stats.norm.pdf(concept_data, concept_mu, concept_sigma)646 concept_ax.plot([concept_data, concept_data], [0, concept_prob], concept_colors[concept_i], linewidth=2)647 concept_ax.scatter(concept_data, concept_prob, color=concept_colors[concept_i], s=50, 648 label=f'PDF at x={concept_data}: {concept_prob:.3f}')649 650 concept_ax.set_xlabel('Data (x)')651 concept_ax.set_ylabel('Probability Density')652 concept_ax.set_title('Density Function Perspective: Fixed Parameters (μ=0, σ=1), Different Data Points')653 654 else: # likelihood perspective655 # likelihood perspective: fixed data, varying parameters656 concept_data_point = 1.5 # fixed observed data657 658 # different possible parameter values (means)659 concept_mus = [-1, 0, 1, 2, 3]660 concept_sigma = 1661 662 # generate x values for multiple pdfs663 concept_x = np.linspace(-4, 6, 1000)664 665 concept_colors = ['#FF9999', '#FFCC99', '#99FF99', '#99CCFF', '#CC99FF']666 667 for concept_i, concept_mu in enumerate(concept_mus):668 concept_pdf = stats.norm.pdf(concept_x, concept_mu, concept_sigma)669 concept_ax.plot(concept_x, concept_pdf, color=concept_colors[concept_i], linewidth=2, alpha=0.7,670 label=f'N({concept_mu}, 1)')671 672 # mark the likelihood of the data point for this param673 concept_likelihood = stats.norm.pdf(concept_data_point, concept_mu, concept_sigma)674 concept_ax.plot([concept_data_point, concept_data_point], [0, concept_likelihood], concept_colors[concept_i], linewidth=2)675 concept_ax.scatter(concept_data_point, concept_likelihood, color=concept_colors[concept_i], s=50, 676 label=f'L(μ={concept_mu}|X=1.5) = {concept_likelihood:.3f}')677 678 # add vertical line at the observed data point679 concept_ax.axvline(x=concept_data_point, color='black', linestyle='--', 680 label=f'Observed data: X=1.5')681 682 concept_ax.set_xlabel('Data (x)')683 concept_ax.set_ylabel('Probability Density / Likelihood')684 concept_ax.set_title('Likelihood Perspective: Fixed Data Point (X=1.5), Different Parameter Values')685 686 elif concept_dist_type_value == "Bernoulli":687 if concept_view_mode == "probability":688 # probability perspective: fixed parameter, two possible data values689 concept_p = 0.3 # fixed parameter690 691 # bar chart for p(x=0) and p(x=1)692 concept_ax.bar([0, 1], [1-concept_p, concept_p], width=0.4, color=['#99CCFF', '#FF9999'], 693 alpha=0.7, label=f'PMF: Bernoulli({concept_p})')694 695 # text showing probabilities696 concept_ax.text(0, (1-concept_p)/2, f'P(X=0|p={concept_p}) = {1-concept_p:.3f}', ha='center', va='center', fontweight='bold')697 concept_ax.text(1, concept_p/2, f'P(X=1|p={concept_p}) = {concept_p:.3f}', ha='center', va='center', fontweight='bold')698 699 concept_ax.set_xlabel('Data (x)')700 concept_ax.set_ylabel('Probability')701 concept_ax.set_xticks([0, 1])702 concept_ax.set_xticklabels(['X=0', 'X=1'])703 concept_ax.set_ylim(0, 1)704 concept_ax.set_title('Probability Perspective: Fixed Parameter (p=0.3), Different Data Values')705 706 else: # likelihood perspective707 # likelihood perspective: fixed data, varying parameter708 concept_data_point = 1 # fixed observed data (success)709 710 # different possible parameter values711 concept_p_values = np.linspace(0.01, 0.99, 100)712 713 # calculate likelihood for each p value714 if concept_data_point == 1:715 # for x=1, likelihood is p716 concept_likelihood = concept_p_values717 concept_ax.plot(concept_p_values, concept_likelihood, 'b-', linewidth=2, 718 label=f'L(p|X=1) = p')719 720 # highlight specific values721 concept_highlight_ps = [0.2, 0.5, 0.8]722 concept_colors = ['#FF9999', '#99FF99', '#99CCFF']723 724 for concept_i, concept_p in enumerate(concept_highlight_ps):725 concept_ax.plot([concept_p, concept_p], [0, concept_p], concept_colors[concept_i], linewidth=2)726 concept_ax.scatter(concept_p, concept_p, color=concept_colors[concept_i], s=50, 727 label=f'L(p={concept_p}|X=1) = {concept_p:.3f}')728 729 concept_ax.set_title('Likelihood Perspective: Fixed Data Point (X=1), Different Parameter Values')730 731 else: # x=0732 # for x = 0, likelihood is (1-p)733 concept_likelihood = 1 - concept_p_values734 concept_ax.plot(concept_p_values, concept_likelihood, 'r-', linewidth=2, 735 label=f'L(p|X=0) = (1-p)')736 737 # highlight some specific values738 concept_highlight_ps = [0.2, 0.5, 0.8]739 concept_colors = ['#FF9999', '#99FF99', '#99CCFF']740 741 for concept_i, concept_p in enumerate(concept_highlight_ps):742 concept_ax.plot([concept_p, concept_p], [0, 1-concept_p], concept_colors[concept_i], linewidth=2)743 concept_ax.scatter(concept_p, 1-concept_p, color=concept_colors[concept_i], s=50, 744 label=f'L(p={concept_p}|X=0) = {1-concept_p:.3f}')745 746 concept_ax.set_title('Likelihood Perspective: Fixed Data Point (X=0), Different Parameter Values')747 748 concept_ax.set_xlabel('Parameter (p)')749 concept_ax.set_ylabel('Likelihood')750 concept_ax.set_xlim(0, 1)751 concept_ax.set_ylim(0, 1)752 753 elif concept_dist_type_value == "Poisson":754 if concept_view_mode == "probability":755 # probability perspective: fixed parameter, different data values756 concept_lam = 2.5 # fixed parameter757 758 # pmf for different x values plot759 concept_x_values = np.arange(0, 10)760 concept_pmf_values = stats.poisson.pmf(concept_x_values, concept_lam)761 762 concept_ax.bar(concept_x_values, concept_pmf_values, width=0.4, color='#99CCFF', 763 alpha=0.7, label=f'PMF: Poisson({concept_lam})')764 765 # highlight a few specific values766 concept_highlight_xs = [1, 2, 3, 4]767 concept_colors = ['#FF9999', '#99FF99', '#FFCC99', '#CC99FF']768 769 for concept_i, concept_x in enumerate(concept_highlight_xs):770 concept_prob = stats.poisson.pmf(concept_x, concept_lam)771 concept_ax.scatter(concept_x, concept_prob, color=concept_colors[concept_i], s=50, 772 label=f'P(X={concept_x}|λ={concept_lam}) = {concept_prob:.3f}')773 774 concept_ax.set_xlabel('Data (x)')775 concept_ax.set_ylabel('Probability')776 concept_ax.set_xticks(concept_x_values)777 concept_ax.set_title('Probability Perspective: Fixed Parameter (λ=2.5), Different Data Values')778 779 else: # likelihood perspective780 # likelihood perspective: fixed data, varying parameter781 concept_data_point = 4 # fixed observed data782 783 # different possible param values784 concept_lambda_values = np.linspace(0.1, 8, 100)785 786 # calc likelihood for each lambda value787 concept_likelihood = stats.poisson.pmf(concept_data_point, concept_lambda_values)788 789 concept_ax.plot(concept_lambda_values, concept_likelihood, 'b-', linewidth=2, 790 label=f'L(λ|X={concept_data_point})')791 792 # highlight some specific values793 concept_highlight_lambdas = [1, 2, 4, 6]794 concept_colors = ['#FF9999', '#99FF99', '#99CCFF', '#FFCC99']795 796 for concept_i, concept_lam in enumerate(concept_highlight_lambdas):797 concept_like_val = stats.poisson.pmf(concept_data_point, concept_lam)798 concept_ax.plot([concept_lam, concept_lam], [0, concept_like_val], concept_colors[concept_i], linewidth=2)799 concept_ax.scatter(concept_lam, concept_like_val, color=concept_colors[concept_i], s=50, 800 label=f'L(λ={concept_lam}|X={concept_data_point}) = {concept_like_val:.3f}')801 802 concept_ax.set_xlabel('Parameter (λ)')803 concept_ax.set_ylabel('Likelihood')804 concept_ax.set_title(f'Likelihood Perspective: Fixed Data Point (X={concept_data_point}), Different Parameter Values')805 806 concept_ax.legend(loc='best', fontsize=9)807 concept_ax.grid(True, alpha=0.3)808 plt.tight_layout()809 plt.gca()810 811 # relevant explanation based on view mode812 if concept_view_mode == "probability":813 concept_explanation = mo.md(814 f"""815 ### Density/Mass Function Perspective816 817 In the **density/mass function perspective**, the parameters of the distribution are **fixed and known**, and we evaluate the function at **different possible data values**.818 819 For the {concept_dist_type_value} distribution, we've fixed the parameter{'s' if concept_dist_type_value == 'Normal' else ''} and shown the {'density' if concept_dist_type_value == 'Normal' else 'probability mass'} function evaluated at different data points.820 821 This is the typical perspective when:822 823 - We know the true parameters of a distribution824 - We want to evaluate the {'density' if concept_dist_type_value == 'Normal' else 'probability mass'} at different observations825 - We make predictions based on our model826 827 **Mathematical notation**: $f(x | \theta)$828 """829 )830 else: # likelihood perspective831 concept_explanation = mo.md(832 f"""833 ### Likelihood Perspective834 835 In the **likelihood perspective**, the observed data is **fixed and known**, and we calculate how likely different parameter values are to have generated that data.836 837 For the {concept_dist_type_value} distribution, we've fixed the observed data point{'s' if concept_dist_type_value == 'Normal' else ''} and shown the likelihood of different parameter values.838 839 This is the perspective used in MLE:840 841 - We have observed data842 - We don't know the true parameters843 - We want to find parameters that best explain our observations844 845 **Mathematical notation**: $L(\theta | X = x)$846 847 /// tip848 The value of $\\theta$ that maximizes this likelihood function is the MLE estimate $\\hat{{\\theta}}$!849 ///850 """851 )852 853 # Display plot and explanation together854 mo.vstack([855 concept_fig,856 concept_explanation857 ])858 return859 860 861@app.cell(hide_code=True)862def _(mo):863 mo.md(r"""864 ## 🤔 Test Your Understanding865 866 Which of the following statements about Maximum Likelihood Estimation are correct? Click each statement to check your answer.867 868 /// details | Probability and likelihood have different interpretations: probability measures the chance of data given parameters, while likelihood measures how likely parameters are given data.869 ✅ **Correct!**870 871 Probability measures how likely it is to observe particular data when we know the parameters. Likelihood measures how likely particular parameter values are, given observed data.872 873 Mathematically, probability is $P(X=x|\theta)$ while likelihood is $L(\theta|X=x)$.874 ///875 876 /// details | We use log-likelihood instead of likelihood because it's mathematically simpler and numerically more stable.877 ✅ **Correct!**878 879 We work with log-likelihood for several reasons:880 1. It converts products into sums, which is easier to work with mathematically881 2. It avoids numerical underflow when multiplying many small probabilities882 3. Logarithm is a monotonically increasing function, so the maximum of the likelihood occurs at the same parameter values as the maximum of the log-likelihood883 ///884 885 /// details | For a Bernoulli distribution, the MLE for parameter p is the sample mean of the observations.886 ✅ **Correct!**887 888 For a Bernoulli distribution with parameter $p$, given $n$ independent samples $X_1, X_2, \ldots, X_n$, the MLE estimator is:889 890 $$\hat{p} = \frac{\sum_{i=1}^n X_i}{n}$$891 892 This is simply the sample mean, or the proportion of successes (1s) in the data.893 ///894 895 /// details | For a Normal distribution, MLE gives unbiased estimates for both mean and variance parameters.896 ❌ **Incorrect.**897 898 While the MLE for the mean ($\hat{\mu} = \frac{1}{n}\sum_{i=1}^n X_i$) is unbiased, the MLE for variance:899 900 $$\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (X_i - \hat{\mu})^2$$901 902 is a biased estimator. It uses $n$ in the denominator rather than $n-1$ used in the unbiased estimator.903 ///904 905 /// details | MLE estimators are always unbiased regardless of the distribution.906 ❌ **Incorrect.**907 908 MLE is not always unbiased, though it often is asymptotically unbiased (meaning the bias approaches zero as the sample size increases).909 910 A notable example is the MLE estimator for the variance of a Normal distribution:911 $$\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (X_i - \hat{\mu})^2$$912 913 This estimator is biased, which is why we often use the unbiased estimator:914 $$s^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i - \hat{\mu})^2$$915 916 Despite occasional bias, MLE estimators have many desirable properties, including consistency and asymptotic efficiency.917 ///918 """)919 return920 921 922@app.cell(hide_code=True)923def _(mo):924 mo.md(r"""925 ## Summary926 927 Maximum Likelihood Estimation really is one of those elegant ideas that sits at the core of modern statistics. When you get down to it, MLE is just about finding the most plausible explanation for the data we've observed. It's like being a detective - you have some clues (your data), and you're trying to piece together the most likely story (your parameters) that explains them.928 929 We've seen how this works with different distributions. For the Bernoulli, it simply gives us the sample proportion. For the Normal, it gives us the sample mean and a slightly biased estimate of variance. And for linear regression, it provides a mathematical justification for the least squares method that everyone learns in basic stats classes.930 931 What makes MLE so useful in practice is that it tends to give us estimates with good properties. As you collect more data, the estimates generally get closer to the true values (consistency) and do so efficiently. That's why MLE is everywhere in statistics and machine learning - from simple regression models to complex neural networks.932 933 The most important takeaway? Next time you're fitting a model to data, remember that you're not just following a recipe - you're finding the parameters that make your observed data most likely to have occurred. That's the essence of statistical inference.934 """)935 return936 937 938@app.cell(hide_code=True)939def _(mo):940 mo.md(r"""941 ## Further Reading942 943 If you're curious to dive deeper into this topic, check out "Statistical Inference" by Casella and Berger - it's the classic text that many statisticians learned from. For a more machine learning angle, Bishop's "Pattern Recognition and Machine Learning" shows how MLE connects to more advanced topics like EM algorithms and Bayesian methods.944 945 Beyond the basics we've covered, you might explore Bayesian estimation (which incorporates prior knowledge), Fisher Information (which tells us how precisely we can estimate parameters), or the EM algorithm (for when we have missing data or latent variables). Each of these builds on the foundation of likelihood that we've established here.946 """)947 return948 949 950@app.cell(hide_code=True)951def _(mo):952 mo.md(r"""953 ## Appendix (helper functions and imports)954 """)955 return956 957 958@app.cell959def _():960 import marimo as mo961 return (mo,)962 963 964@app.cell965def _():966 import numpy as np967 import matplotlib.pyplot as plt968 from scipy import stats969 import plotly.graph_objects as go970 import polars as pl971 from matplotlib import cm972 973 # Set a consistent random seed for reproducibility974 np.random.seed(42)975 976 # Set a nice style for matplotlib977 plt.style.use('seaborn-v0_8-darkgrid')978 return np, plt, stats979 980 981@app.cell(hide_code=True)982def _(mo):983 # Create interactive elements984 true_p_slider = mo.ui.slider(985 start =0.01, 986 stop =0.99, 987 value=0.3, 988 step=0.01, 989 label="True probability (p)"990 )991 992 sample_size_slider = mo.ui.slider(993 start =10, 994 stop =1000, 995 value=100, 996 step=10, 997 label="Sample size (n)"998 )999 1000 generate_button = mo.ui.button(label="Generate New Sample", kind="success")1001 1002 controls = mo.vstack([1003 mo.vstack([true_p_slider, sample_size_slider]), 1004 generate_button1005 ], justify="space-between")1006 return controls, generate_button, sample_size_slider, true_p_slider1007 1008 1009@app.cell(hide_code=True)1010def _(mo):1011 # Create interactive elements for Normal distribution1012 true_mu_slider = mo.ui.slider(1013 start =-5, 1014 stop =5, 1015 value=0, 1016 step=0.1, 1017 label="True mean (μ)"1018 )1019 1020 true_sigma_slider = mo.ui.slider(1021 start =0.5, 1022 stop =3, 1023 value=1, 1024 step=0.1, 1025 label="True standard deviation (σ)"1026 )1027 1028 normal_sample_size_slider = mo.ui.slider(1029 start =10, 1030 stop =500, 1031 value=50, 1032 step=10, 1033 label="Sample size (n)"1034 )1035 1036 normal_generate_button = mo.ui.button(label="Generate New Sample", kind="warn")1037 1038 normal_controls = mo.hstack([1039 mo.vstack([true_mu_slider, true_sigma_slider, normal_sample_size_slider]), 1040 normal_generate_button1041 ], justify="space-between")1042 return (1043 normal_controls,1044 normal_generate_button,1045 normal_sample_size_slider,1046 true_mu_slider,1047 true_sigma_slider,1048 )1049 1050 1051@app.cell(hide_code=True)1052def _(mo):1053 # Create interactive elements for linear regression1054 true_theta_slider = mo.ui.slider(1055 start =-2, 1056 stop =2, 1057 value=0.5, 1058 step=0.1, 1059 label="True slope (θ)"1060 )1061 1062 noise_sigma_slider = mo.ui.slider(1063 start =0.1, 1064 stop =2, 1065 value=0.5, 1066 step=0.1, 1067 label="Noise level (σ)"1068 )1069 1070 linear_sample_size_slider = mo.ui.slider(1071 start =10, 1072 stop =200, 1073 value=50, 1074 step=10, 1075 label="Sample size (n)"1076 )1077 1078 linear_generate_button = mo.ui.button(label="Generate New Sample", kind="warn")1079 1080 linear_controls = mo.hstack([1081 mo.vstack([true_theta_slider, noise_sigma_slider, linear_sample_size_slider]), 1082 linear_generate_button1083 ], justify="space-between")1084 return (1085 linear_controls,1086 linear_generate_button,1087 linear_sample_size_slider,1088 noise_sigma_slider,1089 true_theta_slider,1090 )1091 1092 1093@app.cell(hide_code=True)1094def _(mo):1095 # Interactive elements for likelihood vs probability demo1096 concept_dist_type = mo.ui.dropdown(1097 options=["Normal", "Bernoulli", "Poisson"],1098 value="Normal",1099 label="Distribution"1100 )1101 1102 # Replace buttons with a simple dropdown selector1103 perspective_selector = mo.ui.dropdown(1104 options=["Probability Perspective", "Likelihood Perspective"],1105 value="Probability Perspective",1106 label="View"1107 )1108 1109 concept_controls = mo.vstack([1110 mo.hstack([concept_dist_type, perspective_selector])1111 ])1112 return concept_controls, concept_dist_type, perspective_selector1113 1114 1115if __name__ == "__main__":1116 app.run()1117