finding the closest fit of a normal distribution to a set of sample data points

max-log-lik-2-params.R

library(tidyverse)
library(scico)

set.seed(2)

mu.0 <- runif(1, -2, 2)
sigma.0 <- runif(1, 0.1, 1)
samples <- rnorm(300, mean = mu.0, sd = sigma.0)

list(
  muMax = 0,
  muMin = -2.0,
  muRes = 0.01,
  sigmaMax = 2.0,
  sigmaMin = 0.5,
  sigmaRes = 0.001
) -> opt

mu    <- seq(opt$muMin,    opt$muMax,    opt$muRes)
sigma <- seq(opt$sigmaMin, opt$sigmaMax, opt$sigmaRes)
params <- crossing(mu, sigma)

params %>%
  rowwise() %>%
  mutate(
    lL = sum(log(dnorm(samples, sd = sigma, mean = mu)))
  ) %>%
  ungroup() -> params

params %>%
  arrange(-lL) %>%
  top_n(1) %>%
  as.list() -> peak

params %>%
  ggplot(aes(x = mu, y = sigma, fill = lL)) +
  geom_tile() +
  geom_vline(xintercept = peak$mu, linetype = "dashed", color="deeppink2") +
  geom_hline(yintercept = peak$sigma, linetype = "dashed", color="deeppink2") +
  scale_x_continuous(breaks = c(opt$muMin, peak$mu, opt$muMax)) +
  scale_y_continuous(breaks = c(opt$sigmaMin, peak$sigma, opt$sigmaMax)) +
  scale_fill_scico(palette="lajolla") +
  theme_linedraw() -> plot

print(plot)
cat("mu.0:", mu.0, "\nsigma.0:", sigma.0)