What is the name of the density estimation method where all possible pairs are used to create a Normal mixture distribution? I just thought  of a neat (not necessarily good) way of creating one dimensional density estimates and my question is:
Does this density estimation method have a name? If not, is it a special case of some other method in the literature?
Here is the method: We have a vector $X = [x_1,x_2,...,x_n]$ which we assume is drawn from some unknown distribution we would like to estimate. A way of doing this is to take all possible pairs of values in $X$ and for each pair $[x_i,x_j]_{i \neq j}$ fit a Normal distribution using maximum likelihood. The resulting density estimate is then the mixture distribution that consists of all the resulting Normals, where each Normal is given equal weight.
The figure below illustrates using this method on the vector $[-1.3,0.15,0.73,1.4]$. Here the circles are the datapoints, the coloured Normals are the maximum likelihood distributions estimated using each possible pair and the thick black line shows the resulting density estimate (that is, the mixture distribution).

By the way, it is really easy to implement a method in R that draws a sample from the resulting mixture distribution:
# Generating some "data"
x <- rnorm(30)

# Drawing from the density estimate using the method described above.
density_estimate_sample <- replicate(9999, {
  pair <- sample(x, size = 2)
  rnorm(1, mean(pair), sd(pair))
})

# Plotting the density estimate compared with 
# the "data" and the "true" density.
hist(x ,xlim=c(-5, 5), main='The "data"')
hist(density_estimate_sample, xlim=c(-5, 5), main='Estimated density')
hist(rnorm(9999), xlim=c(-5, 5), main='The "true" density')


 A: This is an intriguing idea, because the estimator of the standard deviation appears to be less sensitive to outliers than the usual root-mean-square approaches.  However, I doubt this estimator has been published.  There are three reasons why: it is computationally inefficient, it is biased, and even when the bias is corrected, it is statistically inefficient (but only a little).  These can be seen with a little preliminary analysis, so let's do that first and then draw the conclusions.
Analysis
The ML estimators of the mean $\mu$ and standard deviation $\sigma$ based on data $(x_i, x_j)$ are
$$\hat\mu(x_i,x_j) = \frac{x_i+x_j}{2}$$
and
$$\hat\sigma(x_i,x_j) = \frac{|x_i-x_j|}{2}.$$
Therefore the method described in the question is
$$\hat\mu(x_1, x_2, \ldots, x_n) = \frac{2}{n(n-1)} \sum_{i\gt j} \frac{x_i+x_j}{2} = \frac{1}{n}\sum_{i=1}^nx_i,$$
which is the usual estimator of the mean, and
$$\hat\sigma(x_1, x_2, \ldots, x_n) = \frac{2}{n(n-1)}\sum_{i\gt j}\frac{|x_i-x_j|}{2} = \frac{1}{n(n-1)}\sum_{i,j}|x_i-x_j|.$$
The expected value of this estimator is readily found by exploiting the exchangeability of the data, which implies $E = \mathbb{E}(|x_i-x_j|)$ is independent of $i$ and $j$.  Whence
$$\mathbb{E}(\hat\sigma(x_1, x_2, \ldots, x_n)) = \frac{1}{n(n-1)}\sum_{i,j}\mathbb{E}(|x_i-x_j|) = E.$$
But since $x_i$ and $x_j$ are independent Normal variates, their difference is a zero-mean Normal with variance $2\sigma^2$.  Its absolute value therefore is $\sqrt{2}\sigma$ times a $\chi(1)$ distribution, whose mean is $\sqrt{2/\pi}$.  Consequently
$$E = \frac{2}{\sqrt{\pi}} \sigma.$$
The coefficient $2/\sqrt{\pi} \approx 1.128$ is the bias in this estimator.
In the same way, but with considerably more work, one could compute the variance of $\hat\sigma$, but--as we will see--there's unlikely to be much interest in this, so I will just estimate it with a quick simulation.
Conclusions


*

*The estimator is biased.  $\hat\sigma$ has a substantial constant bias of about +13%.  This could be corrected.  In this example with a sample size of $n=20,000$ both the biased and bias-corrected estimators are plotted over the histogram.  The 13% error is apparent.


*It is computationally inefficient.  Because the sum of absolute values, $\sum_{i,j}|x_i-x_j|$, has no algebraic simplification, its calculation requires $O(n^2)$ effort instead of the $O(n)$ effort for almost any other estimator.  This scales badly, making it prohibitively expensive once $n$ exceeds $10,000$ or so.  For instance, computing the previous figure required 45 seconds of CPU time and 8 GB RAM in R.  (On other platforms the RAM requirements would be much smaller, perhaps at a slight cost in computation time.)

*It is statistically inefficient.  To give it the best showing, let's consider the unbiased version and compare it to the unbiased version of either the least squares or maximum likelihood estimator
$$\hat\sigma_{OLS} = \sqrt{\left(\frac{1}{n-1} \sum_{i=1}^n \left(x_i - \hat\mu\right)^2\right)} \frac{(n-1)\Gamma((n-1)/2)}{2\Gamma(n/2)}.$$
The R code below demonstrates that the unbiased version of the estimator in the question is surprisingly efficient: across a range of sample sizes from $n=3$ to $n=300$ its variance is usually about 1% to 2% greater than the variance of $\hat\sigma_{OLS}$.  This means you should plan on paying an extra 1% to 2% more for samples in order to achieve any given level of precision in estimating $\sigma$.
Afterward
The form of $\hat\sigma$ is reminiscent of the robust and resistant Theil-Sen estimator--but instead of using the medians of the absolute differences, it uses their means.  If the objective is to have an estimator that is resistant to outlying values or one that is robust to departures from the Normality assumption, then using the median would be more advisable.  

Code
sigma <- function(x) sum(abs(outer(x, x, '-'))) / (2*choose(length(x), 2))
#
# sigma is biased.
#
y <- rnorm(1e3) # Don't exceed 2E4 or so!
mu.hat <- mean(y)
sigma.hat <- sigma(y)

hist(y, freq=FALSE,
     main="Biased (dotted red) and Unbiased (solid blue) Versions of the Estimator",
     xlab=paste("Sample size of", length(y)))
curve(dnorm(x, mu.hat, sigma.hat), col="Red", lwd=2, lty=3, add=TRUE)
curve(dnorm(x, mu.hat, sqrt(pi/4)*sigma.hat), col="Blue", lwd=2, add=TRUE)
#
# The variance of sigma is too large.
#
N <- 1e4
n <- 10
y <- matrix(rnorm(n*N), nrow=n)
sigma.hat <- apply(y, 2, sigma) * sqrt(pi/4)
sigma.ols <- apply(y, 2, sd) / (sqrt(2/(n-1)) * exp(lgamma(n/2)-lgamma((n-1)/2)))

message("Mean of unbiased estimator is ", format(mean(sigma.hat), digits=4))
message("Mean of unbiased OLS estimator is ", format(mean(sigma.ols), digits=4))
message("Variance of unbiased estimator is ", format(var(sigma.hat), digits=4))
message("Variance of unbiased OLS estimator is ", format(var(sigma.ols), digits=4))
message("Efficiency is ", format(var(sigma.ols) / var(sigma.hat), digits=4))

