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.