# Setup/Functions ---------------------------------------------------------

# Loss Function 

loss_fn <- function(x, data, L_p) sum(abs(data - x)^L_p)

# Mode functions to get the (first) mode for discrete or continuous data

mode_discrete <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

mode_continuous <- function(x) {
  lim.inf=min(x)-1; lim.sup=max(x)+1
  s<-density(x,from=lim.inf,to=lim.sup,bw=0.2)
  n<-length(s$y)
  v1<-s$y[1:(n-2)];
  v2<-s$y[2:(n-1)];
  v3<-s$y[3:n]
  ix<-1+which((v1<v2)&(v2>v3))
  s$x[which(s$y==max(s$y))]
}

# Wrapper for optim (numerical optimisation) to simplify code

numerical_solution <- function(L_p, data, fn=loss_fn, init=1) {
  optim(init, 
        fn=loss_fn, data=data, L_p=L_p, 
        method = "Brent", lower = 0, upper = 100)$par
}

n <- 1e6 # sample size (number of simulated values)

# Poisson Data ------------------------------------------------------------

lambda <- 6 # rate parameter

Poisson_data <- rpois(n, lambda)

mean(Poisson_data)
numerical_solution(2, Poisson_data)

median(Poisson_data)
numerical_solution(1, Poisson_data)

mode_discrete(Poisson_data)
numerical_solution(1e-6, Poisson_data)

# Normal (Gaussian) Data --------------------------------------------------

mu <- 12
std_dev <- 3

Gaussian_data <- rnorm(n, mu, std_dev)

mean(Gaussian_data)
numerical_solution(2, Gaussian_data)

median(Gaussian_data)
numerical_solution(1, Gaussian_data)

mode_continuous(Gaussian_data)
numerical_solution(1e-6, Gaussian_data)

# Negative Binomial -------------------------------------------------------

library(stats)

size <- 30
p <- 0.6

neg_binom_data <- rnbinom(n, size = size, prob = p)

mean(neg_binom_data)
numerical_solution(2, neg_binom_data)

median(neg_binom_data)
numerical_solution(1, neg_binom_data)

mode_discrete(neg_binom_data)
numerical_solution(1e-6, neg_binom_data)

# Exponential Distribution ------------------------------------------------

exp_data <- rexp(n, rate = lambda) # simulate data

mean(exp_data)
numerical_solution(2, exp_data)

median(exp_data)
numerical_solution(1, exp_data)

mode_continuous(exp_data)
numerical_solution(1e-6, exp_data)

# Plot of Data ------------------------------------------------------------

par(mfrow = c(2, 2))
hist(Poisson_data, breaks = 20, main = "Poisson", xlab = NULL)
hist(Gaussian_data, breaks = 100, main = "Gaussian", xlab = NULL)
hist(neg_binom_data, breaks = 100, main = "Negative Binomial", xlab = NULL)
hist(exp_data, breaks = 100, main = "Exponential", xlab = NULL)