n <- 10 x <- rnorm(n) f <- function(m, s) sum(log(0.5 * dnorm(x) + 0.5 * dnorm(x, m, s))) vf <- Vectorize(f) m <- sort(c(x, seq(-3, 3, len = 100))) s <- seq(0.0000001, 3, len = 100) ms <- expand.grid(m = m, s = s) ll <- matrix(vf(ms$m, ms$s), length(m)) library(rgl) clear3d() surface3d(unique(ms$m), unique(ms$s), ll / sqrt(n), col = "red") contour(unique(ms$m), unique(ms$s), ll, nlevels = 40)