Recently I discussed the following topic with a friend. The setting is that we have a one-dimensional set of data. (In the example it was points of students, we would be grading.) The goal was to make a density estimation, but not to use anything "fancy" like kernel density estimation, but just use a Gaussian as the estimation. (Sure that makes the big assumption that the data is Gaussian, but that is not the point here.) We discussed two ways:

  1. Make a density estimation using an unsupervised learning method, e.g. using EM-algorithm. In this case the claim is that simply calculating the mean of the data and the standard derivation is already giving one the parameters to get the Gaussian parametrized the right way.
  2. Add up the number of occurrences for each value and then use a supervised learning regression with the Gaussian as function, powered by an optimisation algorithm.

Through discussion we found out that the two clearly have different outcomes - though probably similar. For case 2 we optimise the parameters of the Gaussian such, that the sum of the distances from the Gaussian to the occurrences is minimized (along the y-axis if you will). For case 1 we optimise the parameters along the x-axis if you will.


  1. Preliminary question: Is it correct that the EM-algorithm will have the same result as just calculating mean and standard deviation from the data?

Assuming the answer to the preliminary question is "yes, the results are the same":

  1. What is the intuition interpretation of those two approaches?
  2. None of them can be wrong in itself, but one of them could be wrong in the sense of wrong-usage. Meaning: I want to do something and I use the wrong method for it, because of misunderstood interpretation of what happened. So in that sense: Is one of them wrong in a way?

Example Code

I managed to express myself in R code. As one can see from the plots, the result is definitely not the same for any dataset. Only if the data is Gaussian and large the results get to be similar. But that means little to me ($2^x$ and $exp(x)$ converge both to $infinity$ for $x -> infinity$ and they have little in common otherwise).


std <- 10

datasets <- list(
  data.frame( x = round(rnorm(1000, sd = std))),
  data.frame( x = round(rnorm(1000, sd = std))) %>% filter(x > 0)

for(data1 in datasets){

  data1$densest <- dnorm(data1$x, mean = mean(data1$x), sd = std, log = FALSE)
  data2 <- data1 %>%
    group_by(x) %>%
    summarise(count = n())

  f <- function(x, m, sd, k) {
    k * exp(-0.5 * ((x - m)/sd)^2) # 1/sqrt(2*pi*sd^2) *

  cost <- function(par) {
    rhat <- f(data2$x, par[1], par[2], par[3])
    sum((data2$count - rhat)^2)

  o <- optim(c(0, std, 10), cost, method="BFGS", control=list(reltol=1e-9))
  data1$regr <- f(data1$x, o$par[1], o$par[2], o$par[3])
  data1$regrNormalized <- dnorm(data1$x, mean = o$par[1], sd = abs(o$par[2]), log = FALSE)

  g1 <- ggplot(data=data1, aes(x=x)) +
    geom_point(data=data2, aes(x=x, y=count), alpha=0.2)+
    geom_line(aes(y=densest), color="green") +
    geom_point(data = data.frame(mean = mean(data1$x)), aes(x=mean, y=0), color="green") +
    geom_line(aes(y=regr), color="blue") +
    geom_point(data = data.frame(mean = o$par[1]), aes(x=mean, y=0), color="blue")

  g2 <- ggplot(data=data1, aes(x=x)) +
    geom_line(aes(y=densest), color="green") +
    geom_point(data = data.frame(mean = mean(data1$x)), aes(x=mean, y=0), color="green") +
    geom_line(aes(y=regrNormalized), color="blue") +
    geom_point(data = data.frame(mean = o$par[1]), aes(x=mean, y=0), color="blue")




enter image description here enter image description here enter image description here enter image description here

  • 1
    $\begingroup$ What's wrong with the maximum-likelihood estimators of the mean and variance of your Gaussian given the data? ie sample mean and mean square deviation from the sample mean? Closed form solutions, no messing. Or are you really interested in comparing outcomes of different methods for different (less tractable) distributions? $\endgroup$
    – Spacedman
    Dec 2, 2016 at 9:20
  • $\begingroup$ @Spacedman: Yes I am. It is an actual data science question, not a practical, applied, works-so-I-am-fine-with-it business analyst question. (Apologies, for my cynicism - too many data non-scientists around theses days. You know what I mean ;-).) This is about the intuitive interpretation of the mathematically different approaches. It's a math question. $\endgroup$
    – Make42
    Dec 2, 2016 at 18:49
  • $\begingroup$ If I understand method (2) correctly, you are binning the data first into counts - this loses information. That should be enough for intuition to see that you could get different results - or at least you are estimating different things. $\endgroup$
    – Spacedman
    Dec 4, 2016 at 14:26
  • $\begingroup$ @Spacedman: No, it does not lose information. The set {1,1,1,1,2,2,2,3,4,4,4,4,4,5,5} equals the map {1->4, 2->3, 3->1, 4->5, 5->2} in information. The reason those are different because of different loss functions (where the loss function in 1 is implicit, just as the loss function for K-Means is also just implicit - as an example.) $\endgroup$
    – Make42
    Dec 5, 2016 at 15:16
  • $\begingroup$ ah, I was thrown by your Gaussian assumption again. A set of integers can't be from a Gaussian because continuous. But you don't care about that. If you could include some code for what you did it might be possible to better see what you are doing and if the asymptotic results (for large amounts of data) converge - because its still a bit unclear. $\endgroup$
    – Spacedman
    Dec 5, 2016 at 15:24

2 Answers 2


The difference comes about because the sample mean and sd estimators are assuming the data is a random sample from the distribution on the full range of the distribution, in this case -Inf to +Inf. In your second case above there's no data below zero. The sample mean estimator can't see that you've only got half a Gaussian there, it will always be a positive number.

On the other hand, the curve-fitting approach of the least-squares fitting procedure (which is presumably possible via EM, or any of a dozen other algorithms) can see the data as part of a curve that has a Gaussian density shape. If the data is just curving upwards, it will fit to the left side of a Gaussian. If its curving down it will fit to the right side. If its flat it will either fit it to the peak with a very large variance, or a tail and might even go degenerate and divide by zero.

So you are doing parametric estimation of a distribution, or curve-fitting, which are categorically different things.


Assuming that your data is really Gaussian, the best way to fit to a Gaussian is simply to calculate the mean and variance of your data (sample mean and sample variance) and to use those parameters Gaussian fit. I.e. fit the moments. No need for any complex E-M estimation or supervised learning.

The only possible weakness of fitting using estimated moments is that if your sample is small, you might be overly sensitive to outlier points, so you have to be careful about that.

  • $\begingroup$ I don't think you answered any of my questions. I edited the post to make sure, you know what the questions are: It is clear that the different Gaussians mean different things, because they use different loss functions. It is not about "right" and "wrong", but about using the right one for the right goal. To know what methods to use for what goal, we need an intuitive understanding of what the methods mean. $\endgroup$
    – Make42
    Dec 2, 2016 at 18:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.