i already do for one single mle sample size which is n=50 for OBPL-distribution, but now i wanna implement the loop simulation to do replicates R=1000times on different random samples of various sizes of n=50, 300,1050,2050. and calc average MLEs,biases,and mean square errors(MSEs)are obtained at each sample size . but stuck . This is the single mle code
Load required libraries
library(stats)
Define the inverse regularized incomplete beta function
inv_beta <- function(u, a, b) {
qbeta(u, shape1 = a, shape2 = b) # Regularized beta function in R
}
Define the quantile function for OBP-Logistic distribution
quantile_obp_logistic <- function(u, a, b, mu, s) {
Compute K
K <- inv_beta(u, a, b) / (1 + inv_beta(u, a, b))
Compute x using the quantile function
x <- mu - s * log((1 - K) / K)
return(x)
}
Generate random samples for OBP-Logistic
generate_obp_logistic <- function(n, a, b, mu, s) {
u <- runif(n) # Generate uniform random variables
x <- sapply(u, quantile_obp_logistic, a = a, b = b, mu = mu, s = s)
return(x)
}
Set parameters for the OBP-Logistic distribution (Set 1)
a <- 0.5 # Shape parameter a
b <- 1.5 # Shape parameter b
mu <- 0.75 # Location parameter
s <- 0.5 # Scale parameter
n <- 50 # Number of samples
Generate random samples
set.seed(123) # For reproducibility
random_samples <- generate_obp_logistic(n, a, b, mu, s)
random_samples
Plot density
plot(density(random_samples),
main = "Density Plot of Random Samples from OBP-Logistic",
xlab = "x",
col = "blue",
lwd = 2)
Define the log-likelihood function for OBP-Logistic
log_likelihood <- function(params, data) {
a <- params[1]
b <- params[2]
mu <- params[3]
s <- params[4]
Ensure parameters are valid
if (a <= 0 || b <= 0 || s <= 0) {
return(-Inf)
}
Calculate log-likelihood
n <- length(data)
term1 <- n * log(gamma(a + b)) - n * log(s) - n * log(gamma(a)) - n * log(gamma(b))
term2 <- a * sum(log(1 / exp(-(data - mu) / s)))
term3 <- -(a + b) * sum(log(1 + 1 / exp(-(data - mu) / s)))
log_likelihood <- term1 + term2 + term3
return(-log_likelihood) # Negative for minimization
}
Perform MLE to fit the OBP-Logistic distribution
Initial parameter estimates
initial_params <- c(a = 0.5, b = 1.5, mu = mean(random_samples), s = sd(random_samples))
Optimization using optim
result <- optim(
par = initial_params,
fn = log_likelihood,
data = random_samples,
method = "L-BFGS-B",
lower = c(0.01, 0.01, -Inf, 0.01), # Lower bounds for parameters
upper = c(Inf, Inf, Inf, Inf) # Upper bounds for parameters
)
Extract estimated parameters
estimated_params <- result$par
names(estimated_params) <- c("a", "b", "mu", "s")
print("Estimated Parameters:")
print(estimated_params)
Summary of optimization
print("Optimization Result:")
print(result)