Loss given Activation Function and Probability Model
At the output we can have an activation that takes a weighted linear sum from the previous layer and transforms it into a new range. Here are a few activations you could use:

$x = x$, the identity activation which is what is usually used at the end for regression models

$\sigma(x) = \frac{1}{1+e^{-x}}$ which maps to $(0, 1)$

$exp(x) = e^x$ which maps to $(0,\infty)$

These activations are differentiable which means backpropagation will keep working. Now that we can manipulate our output to turn into ranges that we want, we can use them to use them for probability distributions. In the $\sigma$ case we can use them directly as the probability at event $Y$ as opposed to not $Y$, so $P(Y|X)$. We can also use multiple output nodes to parameterize a probability distribution. For example let's try to approximate a Gaussian distribution conditioned on our input features $x$. Instead of regressing directly on $y$ we can regress $\mu$ and $\sigma$ that define the Gaussian distribution. For $\mu$ we can use the identity activation because $\mu$ is not constrained, $\sigma$ needs to be strictly positive so we could use $e^x$ for that.

However, in our training data, we don't have $\mu$ and $\sigma$, we only have a $y$, a sample from this probability distribution. What we do know is that given a specific $\mu$ and $\sigma$ what the likelihood is of that $y$. A value of 100 is extremely unlikely from a normal distribution with $\mu = 20$ and $\sigma=20$, but a value of $20$ has a much higher density. The density of a Gaussian is $f(y|\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(y-\mu)^2}{2\sigma^2}}$. This means we can check the likelihood of a $y$ given our current estimate of $\mu$ and $\sigma$ given $x$ and try and make the likelihood higher. Or minimize the negative log-likelihood, which has less numerical issues as opposed to normal likelihood.