Do we also need to model a probability distribution for the decoder of a VAE? I'm working on understanding VAEs, mostly through video lectures of Stanford cs231n, in particular lecture 13 tackles on this topic and I think I have a good theoretical grasp.

However, when looking at actual code of implementations, such as this code  from this blog of VAEs I see some differences which I can't quite understand.

Please take a look at this VAE architecture visualization from the class, specifically the decoder part. From the way it is presented here I understand that the decoder network outputs mean and covariance for the data distribution. To get an actual output (i.e. image) we need to sample from the distribution that is parametrized by mean and covariance - the outputs of the decoder.

Now if you look at the code from the Keras blog VAE implementation, you will see that there is no such thing. A decoder takes in a sample from latent space and directly maps its input (sampled z) to an output (e.g. image), not to parameters of a distribution from which an output is to be sampled. 

Am I missing something or does this implementation not correspond to the one presented in the lecture? I've been trying to make sense of it for quite some time now but still can't seem to understand the discrepancy.
Thanks @nbro for pointing this out.

The pictorial architecture in the slides uses the Gaussian loss, which when coupled with Maximum Likelihood Estimation gives the squared error loss (not to remove any tractability issues). The main reason we do the encoder Gaussian trick is to force the latent variable $z$ to be normal so that we can apply $KL$ $Divergence$ to optimise an otherwise intractable integral. You can get a better intuition and reasoning in this video. 

The pictorial architecture is basically taking the Gaussian loss so that the final loss becomes the squared error loss effectively. Also the loss term used in your blog link is exactly the same loss term used in the original paper, but the blog is using CE loss (it is the more common loss used for classification). I am not sure how they are using the CE loss though, since it is only valid for $0$ and $1$ values and AFAIK MNIST data-set has grayscale images. 

I am not exactly sure how they implement the randomness of the Gaussian loss in the decoder structure, but in the simplest of cases they just take the MSE

Check out this blog on VAE's (where they have taken the mean $\Sigma$ which they have abbreviated as mean, I have not checked their implementation detail to know what they exactly mean by that) and also this answer on Data Science on implementation of VAE's4 (both of which gives a more general form of loss). Also for the exact Mathematics check out Appendix C of the original paper.