How to formalize or calculate parameters that maximize the likelihood without numerical method?

Question

I am trying to understand how can we manually calculate the parameter with highest likelihoods. I don't want to use Maximum likelihood or any numerical method and want to understand the proccess. My explanation might be wrong, kindly let me know if I haven't explained it right.

If I have heights observations $X= [45,55]$ and have first parameter values of $\mu = [40,60]$ and $\sigma = [1,2]$

Likelihood can be given as $P(X|\theta) = N(X;\theta)$ in this case.

If we assume the observations i.i.d, what I understand is that we will have:

$$L(\theta|X) = \prod_{i=1}^{N} N(x_i;\theta) $$

Do I need to loop (product/summation) for each combination of $\mu$ and $sigma$ and argmax? How can I formulate this?

What I understand is, I will will take all combinations of these parameters and evaluate for X and see which give maximum likelihood, but I am not sure about two values of X, should I take average of that?

Other thing that confuses me is regarding my other question, if I need to go through each combination of $\mu$ and $\sigma$ than isn't it same that we do for marginal likelihood $P(X)$ of Bayes rule, so aren't nominator(likelihood*prior) and marginal likelihood (denominator) in Bayes equally difficult to compute?

Update: I know how to derivate and know that it is effecient to estimate this with MLE, I just wanted make sense of difficulty as I need to compare it with marginal likelihood/evidence and see why they both arent same. I have just written the code to show what i mean by manual/brute-force approach here.

import numpy as np
import scipy.stats as stats
import math
arr = []
max_l = max_m = max_s = 0

for m in [40,60]:
  for s in [1,2]:
    if stats.norm.pdf(55,m,s)>max_l:
      print('max is ',m,s)
      max_m = m
      max_s = s 
      max_l = stats.norm.pdf(55,m,s)

print(f"mu:{max_m},sigma:{max_s},max_likilihood:{max_l}")

This produces

mu:60,sigma:2,max_likilihood:0.00876415024678427

I am not sure, how to integrate different $X= [45,55]$ in this loop. Maybe multiple the two evaulation inside second loop stats.norm.pdf(55,m,s)*stats.norm.pdf(45,m,s)?

Answer 1

You’re asking how to estimate the maximum likelihood parameters of a normal distribution from data. That is—you want the parameters with highest likelihood, given the data.

How you find these is up to you; there are several optimization procedures you could turn to. You could start at some initial guess of $\mu$ and $\sigma$, then climb the gradient of the log-likelihood to iteratively improve these values. You could even use the Hessian, to similar effect. (These work because the Gaussian likelihood is convex.) But you said you don’t want that. Besides, there’s a faster, exact solution. So let’s try something else.

It turns out, there are closed-form expressions for the optimal parameters.

$$\hat{\mu} = \overline{x} \equiv \frac{1}{n}\sum_{i=1}^n x_i$$

$$\hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \overline{x})^2$$

The best mu is just the sample mean. The best sigma is just the standard deviation.

---

(I did trick you a bit here. There’s no way to get out of using calculus—because the maximum is where the gradient of the likelihood is zero. But you can trust, for a Gaussian distribution, that someone’s already done that work. If you trust this, you can use their closed-form solution.)

---

Edit: Okay, you've altered your question pretty substantially—it's a separate question. I'm going to answer the updated one below.

You're no longer seeking the overall maximum likelihood parameters. You are trying to select from a finite set of possible parameter assignments the one that has highest likelihood.

Your loop is entirely reasonable. As your likelihood equation (for $L(\theta \mid X)$) shows, you should multiply the probabilities of each observation. You can do it the way you showed. Another slightly cleaner way is this Python code:

import numpy as np
import scipy.stats as stats
import math
max_l = max_m = max_s = 0
observations = [45, 55]  # Let's keep these together.

for m in [40, 60]:
  for s in [1, 2]:
    # Compute the product of all observations' probabilities.
    likelihood = stats.norm.pdf(observations, loc=m, scale=s).prod()
    if likelihood > max_l:
      print('max is ', m, s)
      max_m = m
      max_s = s 
      max_l = likelihood

print(f"mu:{max_m}, sigma:{max_s}, max_likelihood:{max_l}")

If you have many samples, you'll quickly find yourself at the limits of computers' floating-point representation of numbers. You should instead do all of your computations in log-space.