Why is everything based on likelihoods even though likelihoods are so small?

The key lies not in the absolute size of the likelihood values but in their relative comparison and the mathematical principles underlying likelihood-based methods. The smallness of the likelihood is expected when dealing with continuous distributions and a product of many probabilities because you're essentially multiplying a lot of numbers that are less than 1.

The utility of likelihoods comes from their comparative nature, not their absolute values. When we compare likelihoods across different sets of parameters, we're looking for which parameters make the observed data "most likely" relative to other parameter sets, rather than looking for a likelihood that suggests the data is likely in an absolute sense.

The scale of likelihood values is often less important than how these values change relative to changes in parameters. This is why in many statistical methods, such as MLE, we're interested in finding the parameters that maximize the likelihood function, as these are considered the best estimates given the data.

Because likelihood values can be extremely small, in practice, statisticians often work with the log of the likelihood. This transformation turns products into sums, making the values more manageable and the optimization problems easier to solve, while preserving the location of the maximum.

set.seed(123)
random_numbers <- rnorm(50, mean = 5, sd = 5)

# Function to calculate log likelihood of a normal distribution
log_likelihood <- function(data, mean, sd) {
  sum(dnorm(data, mean, sd, log = TRUE))
}

# Calculating log likelihood for the correct parameters
log_likelihood_correct <- log_likelihood(random_numbers, 5, 5)
print(log_likelihood_correct)
[1] -147.4507

# Calculating log likelihood for incorrect parameters
log_likelihood_incorrect <- log_likelihood(random_numbers, 6, 6)
print(log_likelihood_incorrect)
[1] -150.5959

# Comparison
print(log_likelihood_correct > log_likelihood_incorrect)
[1] TRUE

First, as others have mentioned, we usually work with the logarithm of the likelihood function, for various mathematical and computational reasons.

Second, since the likelihood function depends on the data, it is convenient to transform it to a function with standardized maxima (see Pickles 1986).

set.seed(123)
random_numbers <- rnorm(50, mean = 5, sd = 5)

max_likelihood <- prod(dnorm(random_numbers, mean = 5, sd = 5))

nonmax_likelihood <- rep(0,1000)
j <- 1 

for (k in seq(0,10,length.out=1000)) {
  nonmax_likelihood[j] <- prod(dnorm(random_numbers, mean=k, sd=5))
  j <- j+1
}

par(mfrow = c(1, 2))

plot(seq(0,10,length.out=1000),nonmax_likelihood/max_likelihood, 
     xlab="Mean", ylab="Relative likelihood")

plot(seq(0,10,length.out=1000),log(nonmax_likelihood) - log(max_likelihood),
     xlab="Mean", ylab="Relative log-likelihood")

I can think of two things that might help you.

First, likelihoods are defined only to a proportionality factor and their utility comes from their use in a ratio and while they are proportional to the relevant probability, they are not probabilities. That means that if you are uncomfortable with the values in the range of 10−6510−65 then you could simply multiply them all by 10651065 without changing the ratios. Of course, there is no need to do as the ratio effectively does it for you. The likelihood ratio for the two distributions is about 25 times in favour of the 5,5 distribution over the 6,6 distribution. That would typically be thought of as being fairly strong (but not overwhelmingly strong) support by the data (and the statistical model) for the 5,5 distribution over the 6,6 distribution.

Second, I usually find a plot of the likelihood as a function of a parameter to be helpful. You have set up the system with two parameters that are effectively 'of interest' and so the relevant likelihood function would be three dimensional and thus awkward. (Those dimensions being the population mean, the standard deviation, and the likelihood values.) It would be easier for you to fix one of those parameters and explore the likelihoods as a function of the other. My justification for looking at the full likelihood function rather than a singular ratio of two selected points in parameter space is that it contains more information and it allows the data to speak with less distortion.