My R Code for the MVUE post

library(tidyverse)

library(clipr)

sample.n <- 1e4

#Random Parameters

sim <- data.frame("mu" = runif(sample.n, min = 1, max = 4.5),

"sigma" = runif(sample.n, min = 0.5, max = 1.5),

"size" = sample(6:10, sample.n, replace = TRUE))

#Consistant Paramaters

sim <- data.frame("mu" = rep(2.4, sample.n),

"sigma" = rep(1, sample.n),

"size" = rep(1e4, sample.n))

sim <- sim %>%

rowwise() %>%

mutate(samples = list(round(rlnorm(size, mu, sigma),2)),

true_mean = exp(mu + sigma^2 /2),

sample_mean = mean(unlist(samples)),

sample_mean_diff = sample_mean - true_mean,

MLE_mean = exp(mean(log(unlist(samples))) + sd(log(unlist(samples)))^2 / 2),

MLE_mean_diff = MLE_mean - true_mean,

MVUE_mean = exp(mean(log(unlist(samples)))) *

(1 + (size - 1) * (sd(log(unlist(samples)))^2 /2)/size +

(size - 1)^3 * (sd(log(unlist(samples)))^2 /2)^2 / (factorial(2) * size^2 * (size+1)) +

(size - 1)^5 * (sd(log(unlist(samples)))^2 /2)^3 / (factorial(3) * size^3 * (size+1) * (size+3)) +

(size - 1)^7 * (sd(log(unlist(samples)))^2 /2)^4 / (factorial(4) * size^4 * (size+1) * (size+3) * (size+5)) +

(size - 1)^9 * (sd(log(unlist(samples)))^2 /2)^5 / (factorial(5) * size^5 * (size+1) * (size+3) * (size+5) * (size+7))),

MVUE_mean_diff = MVUE_mean - true_mean,

MVUE_MLE_diff = MLE_mean - MVUE_mean,

MVUE_MLE_pro = MLE_mean / MVUE_mean

)

#BIAS

(sample_mean_BIAS <- mean(sim$sample_mean_diff))

(MLE_mean_BIAS <- mean(sim$MLE_mean_diff))

(MVUE_mean_BIAS <- mean(sim$MVUE_mean_diff))

#Variance

sample_mean_VAR <- var(sim$sample_mean)

MLE_mean_VAR <- var(sim$MLE_mean)

MVUE_mean_VAR <- var(sim$MVUE_mean)

#MSE

sample_mean_MSE <- var(sim$sample_mean) + sample_mean_BIAS^2

MLE_mean_MSE <- var(sim$MLE_mean) + MLE_mean_BIAS^2

MVUE_mean_MSE <- var(sim$MVUE_mean) + MVUE_mean_BIAS^2

#to copy and paste into IHstat

data <- unlist(sim$samples[1])

write_clip(data)

ggplot() +

geom_density(data = sim, aes(x = MLE_mean, colour = "MLE")) +

geom_density(data = sim, aes(x = MVUE_mean, colour = "MVUE")) +

geom_density(data = sim, aes(x = sample_mean, colour = "Sample")) +

xlab("Estimate Value") +

ggtitle("Distribution of Mean Estimators", "10,000 Samples") +

geom_vline(xintercept = exp(2.4 + 1^2 /2))