I have always seen statistical inferences being made on "clean and exact" probability distributions - for example, (multivariate) normal distributions. But can statistical inferences be made on "rough and irregular" probability distributions? I tried to illustrate my question using the R programming language:

Suppose we have the following data:

var_1 = rnorm(1000, 10, 5)
var_2 = rnorm(1000, 5, 10)
my_data = data.frame(var_1, var_2)

      var_1      var_2
1  7.197622  -4.957987
2  8.849113  -5.399550
3 17.793542   4.820198
4 10.352542   3.678249
5 10.646439 -20.493428
6 18.575325  15.405735

I can try and estimate the probability distribution for this data using two different methods: kernel density estimation and mixture models.

1) Kernel density estimates

Using R, we can estimate the multivariate kernel density estimate of this data and generate a probability distribution function:

# Load library

# Perform kernel density estimate (choose arbitrary bandwidth parameters)

H <- diag(c(1.25, 0.75))
kde <- ks::kde(x = my_data, H = H)

# Plot the results
image(kde$eval.points[[1]], kde$eval.points[[2]], kde$estimate,
      col = viridis::viridis(20), xlab = "x", ylab = "y")
points(kde$x) # The data is returned in $x

Enter image description here


kd <- with(my_data, MASS::kde2d(my_data$var_1, my_data$var_2, n = 50))

 plot_ly(x = kd$x, y = kd$y, z = kd$z) %>% add_surface()

Enter image description here

I am not sure if it is possible to evaluate how well this kernel density estimate (i.e., the choice of the bandwidth/smoothing parameters) fits the data (e.g., aic, bic, and likelihood) - in reality, we would fit multiple kernel density estimates to this same data using different choices of kernels, smoothing and bandwidth parameters. The final kernel density estimate would create a probability distribution function of the following form:

Enter image description here

2) Mixture models

Another popular approach is to fit (Gaussian) mixture models to this data:


# Fit mixture model with two components:
mod1 <- Mclust(my_data, G = 2)

Enter image description here

  summary(mod1, parameters = TRUE)

We can see the results of the mixture:

 summary(mod1, parameters = TRUE)
Gaussian finite mixture model fitted by EM algorithm

Mclust EEI (diagonal, equal volume and shape) model with two components:

 log-likelihood    n df       BIC      ICL
      -6747.084 1000  7 -13542.52 -14181.7

Clustering table:
  1   2
750 250

Mixing probabilities:
        1         2
0.6613209 0.3386791

          [,1]     [,2]
var_1 9.438299 11.33490
var_2 2.250080 11.62347

         var_1    var_2
var_1 23.75622  0.00000
var_2  0.00000 82.16364
         var_1    var_2
var_1 23.75622  0.00000
var_2  0.00000 82.16364

The log-likelihood of this fitted mixture model (assume no covariances) is -6747.084, and the general form of this model is: 0.66 * Bivariate_Normal (mean = (9.4, 2.2), variance = (23.75622, 82)) + 0.33 * Bivariate_Normal (mean = (11.3, 11.6), variance = (23.7, 82))

I am not sure if my understanding of multivariate theory is correct, but I have plotted (what I think is) the mixture model (i.e., the resulting probability distribution):


bivar <- bmbvpdf (
    mean.X1 = 0.66*9.4 + 0.33*11, mean.X2 = 0.66*2.2 + 0.33*11.6,

    sd.X1 = sqrt(0.66* 23 + 0.33 * 82), sd.X2 = sqrt(0.66* 23 + 0.33 * 82),


Enter image description here

If once the probability distribution function for some data is estimated using kernel density estimates or by mixture models, can we use these resulting probability distributions for statistical inference? Is it even possible to take the kernel density estimate and find the conditional probability distribution (e.g., Prob(x1 given x2 = c))?

  • $\begingroup$ Note: I was trying to find if there exist any evaluation metrics which can be used to check how well the kernel density estimate fits the data (e.g. bic, aic, likelihood). I found a useful link over here: cran.r-project.org/web/packages/kedd/vignettes/kedd.pdf . Has anyone ever heard of evaluation metrics for kernel density estimates? thanks! $\endgroup$
    – stats555
    Oct 17 at 1:06
  • 2
    $\begingroup$ what kind of statistical inference do you have mind? could you go into more detail? $\endgroup$
    – user551504
    Oct 17 at 1:42
  • 1
    $\begingroup$ Yes. In fact some tests are based on empirical distribution functions. Two-sample Kolmogorov-Smirnov test, for example. Contingency table tests, for categorical data for another. $\endgroup$
    – Alexis
    Oct 17 at 1:46
  • 2
    $\begingroup$ You could use a permutation test to see if two empirical distributions that seem to have similar shapes have different means, different variances, etc. You could use a nonparametric bootstrap confidence interval to estimate the mean of a single univariate empirical distribution. $\endgroup$
    – BruceET
    Oct 17 at 2:12
  • $\begingroup$ thank you everyone for your replies! $\endgroup$
    – stats555
    Oct 17 at 2:12

You have some very nice graphics in your question. I am not entirely sure I understand exactly what each illustrates in terms of your main question.

Ignoring that your two variables are sampled from normal distributions, I want to do an elementary exploration these two samples: (a) Are they distinguishable by standard tests. (b) If they are mixed, what is the mean of the mixture. A few of the methods and challenges of doing statistical inference on data from unknown population types are discussed briefly.


A Kolmogorov-Smirnov test finds a clear difference between the two samples. This test is notorious for low power, but we have large enough samples that it is useful

x1 = rnorm(1000,10,5)
x2 = rnorm(1000,5,10)
ks.test(x1, x2)

        Two-sample Kolmogorov-Smirnov test

data:  x1 and x2
D = 0.33, p-value < 2.2e-16
alternative hypothesis: two-sided

The test statistic of the K-S test is the largest vertical distance between the empirical CDFs of the two samples.

hdr = "ECDF Plots of Samples 1 (blue) and 2"
plot(ecdf(x2), col="brown", main=hdr)
lines(ecdf(x1), col="blue")

enter image description here

A Welch 2-sample t test shows a highly significant difference in means, and an F-test shows a highly significant difference in variances. If data are nearly normal the Welch test is valid (especially on account of the moderately large sample size). The F test is not at all robust against departures from normality.

t.test(x1, x2)$p.val
[1] 4.329103e-37

var.test(x1, x2)$p.val
[1] 1.202862e-103


Mix the two samples to obtain a single sample x of size $n=1000.$ Because the difference in locations of x1 and x2 is small compared to their variances, the mixture in unimodal.

x = sample(c(x1, x2), 1000)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-21.040   3.278   8.668   7.676  12.710  38.904 
hist(x, prob=T, col="skyblue2")

enter image description here

Confidence intervals of the mean (t interval) and median (Wilcoxon) of the mixture distribution. Both styles of confidence interval are questionable: the t interval $(7.17, 8.19)$ for the mean because the mixture distribution is not normal, and the nonparametric interval $(7.57, 8.54)$ for the median because the mixture distribution is not symmetrical.

[1] 7.165765 8.185229
[1] 0.95

wilcox.test(x, conf.int=T)$conf.int
[1] 7.568812 8.538802
[1] 0.95

One style of 95% nonparametric bootstrap CI $(7.19, 8.17)$ is shown below, it is similar to the t interval.

a.obs = mean(x); a.obs
[1] 7.675497
re.d = replicate(2000, mean(sample(x,1000,rep=T))-a.obs)
LU = quantile(re.d, c(.975,.025))
a.obs - LU
   97.5%     2.5% 
7.186136 8.166439 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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