# random number generation
set.seed(123)
n_sim <- 10000
n <- 50
mu <- 10
sigma <- 2
estimates <- replicate(n_sim, {
sample <- rnorm(n, mean = mu, sd = sigma)
mean(sample)
})
# Estimate the bias
bias <- mean(estimates) - mu
cat("Bias:", bias, "\n")
# Increase sample size
sample_sizes <- c(10, 50, 100, 500, 1000)
true_mean <- 10
sd <- 2
n_sim <- 1000
# Store means and SDs
mean_estimates <- numeric(length(sample_sizes))
sd_estimates <- numeric(length(sample_sizes))
for (i in seq_along(sample_sizes)) {
n <- sample_sizes[i]
estimates <- replicate(n_sim, mean(rnorm(n, mean = true_mean, sd = sd)))
mean_estimates[i] <- mean(estimates)
sd_estimates[i] <- sd(estimates)
}
# Draw Consistency plot
png("consistency_plot.png", width = 800, height = 600)
plot(sample_sizes, mean_estimates, type = "b", pch = 19, col = "blue",
ylim = c(true_mean - 1, true_mean + 1),
xlab = "Sample Size (n)", ylab = "Estimated Mean",
main = "Check for Consistency of Sample Mean")
abline(h = true_mean, col = "red", lty = 2, lwd = 2)
arrows(sample_sizes, mean_estimates - sd_estimates,
sample_sizes, mean_estimates + sd_estimates,
angle = 90, code = 3, length = 0.05, col = "gray")
legend("topright", legend = c("Estimated Mean", "True Mean = 10"),
col = c("blue", "red"), pch = c(19, NA), lty = c(1, 2), lwd = c(2, 2))
dev.off()
# checking efficiency
var_mean <- var(replicate(1000, mean(rnorm(50, mu, sigma))))
var_median <- var(replicate(1000, median(rnorm(50, mu, sigma))))
cat("Var(mean):", var_mean, " Var(median):", var_median, "\n")
## simulate random number sequentially(Using MLE method)
simulate_mle <- function(dist = "normal", true_param, n = 50, n_sim = 1000) {
estimates <- numeric(n_sim)
for (i in 1:n_sim) {
x <- switch(dist,
"normal" = rnorm(n, mean = true_param, sd = 1),
"binomial" = rbinom(n, size = 1, prob = true_param),
"exponential" = rexp(n, rate = true_param),
stop("Unsupported distribution"))
estimates[i] <- switch(dist,
"normal" = mean(x),
"binomial" = mean(x),
"exponential" = 1 / mean(x))
}
return(estimates)
}
evaluate_properties <- function(estimates, true_param) {
list(
bias = mean(estimates) - true_param,
variance = var(estimates),
mse = mean((estimates - true_param)^2)
)
}
# Settings of sequence random number
distributions <- c("normal", "binomial", "exponential")
true_params <- c(normal = 5, binomial = 0.6, exponential = 2)
n <- 50
n_sim <- 1000
for (dist in distributions) {
cat("\n---", toupper(dist), "Distribution ---\n")
estimates <- simulate_mle(dist = dist, true_param = true_params[[dist]], n = n, n_sim = n_sim)
props <- evaluate_properties(estimates, true_params[[dist]])
cat("Bias:", props$bias, "\n")
cat("Variance:", props$variance, "\n")
cat("MSE:", props$mse, "\n")
}
# Loop over distributions
for (dist in distributions) {
cat("\n---", toupper(dist), "Distribution ---\n")
estimates <- simulate_mle(dist = dist, true_param = true_params[[dist]], n = n, n_sim = n_sim)
props <- evaluate_properties(estimates, true_params[[dist]])
cat("Bias: ", props$bias, "\n")
cat("Variance: ", props$variance, "\n")
cat("MSE: ", props$mse, "\n")
# Save histogram to PDF
file_name <- paste0("MLE_Histogram_", dist, ".pdf")
pdf(file = file_name, width = 7, height = 5)
hist(estimates,
main = paste("MLE Estimates -", toupper(dist)),
xlab = "Estimated Value",
col = "skyblue", border = "black", breaks = 30)
abline(v = true_params[[dist]], col = "red", lwd = 2, lty = 2)
legend("topright", legend = paste("True =", true_params[[dist]]),
col = "red", lwd = 2, lty = 2)
dev.off()
cat("Histogram saved to", file_name, "\n")
}
