Fitting Gaussian to data: Density-Estimation vs Regression

Question

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.

Questions

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":

2. What is the intuition interpretation of those two approaches?
3. 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).

library(ggplot2)
library(dplyr)
library(tidyr)

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")

  plot(g1)
  plot(g2)

}

Plots

Answer 1

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.

Answer 2

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.