# \donttest{
## Example 1: Basic Distribution Functions
library(gkwdist)
# Set parameters for GKw distribution
alpha <- 2.0
beta <- 3.0
gamma <- 1.5
delta <- 2.0
lambda <- 1.2
# Create sequence of x values
x <- seq(0.01, 0.99, length.out = 200)
# Compute density
dens <- dgkw(x, alpha, beta, gamma, delta, lambda)
# Compute CDF
cdf <- pgkw(x, alpha, beta, gamma, delta, lambda)
# Compute specific quantiles
probs <- c(0.05, 0.25, 0.5, 0.75, 0.95)
quants <- qgkw(probs, alpha, beta, gamma, delta, lambda)
print(round(quants, 4))
# Generate random sample
set.seed(123)
sample <- rgkw(1000, alpha, beta, gamma, delta, lambda)
# PDF
plot(x, dens,
type = "l", lwd = 2, col = "darkblue",
main = "GKw Probability Density Function",
xlab = "x", ylab = "f(x)", las = 1
)
grid(col = "gray80", lty = 2)
# CDF
plot(x, cdf,
type = "l", lwd = 2, col = "darkred",
main = "GKw Cumulative Distribution Function",
xlab = "x", ylab = "F(x)", las = 1
)
grid(col = "gray80", lty = 2)
# Histogram of random sample
hist(sample,
breaks = 30, probability = TRUE,
col = "lightblue", border = "white",
main = "Random Sample from GKw",
xlab = "x", ylab = "Density", las = 1
)
lines(x, dens, col = "darkblue", lwd = 2)
grid(col = "gray80", lty = 2)
# Q-Q plot
theoretical <- qgkw(ppoints(length(sample)), alpha, beta, gamma, delta, lambda)
empirical <- sort(sample)
plot(theoretical, empirical,
pch = 19, col = rgb(0, 0, 1, 0.3),
main = "Q-Q Plot",
xlab = "Theoretical Quantiles", ylab = "Sample Quantiles", las = 1
)
abline(0, 1, col = "red", lwd = 2, lty = 2)
grid(col = "gray80", lty = 2)
## Example 2: Comparing Distribution Families
# Create comparison plot
x <- seq(0.001, 0.999, length.out = 500)
# GKw (5 parameters) - most flexible
plot(x, dgkw(x, 2, 3, 1.5, 2, 1.2),
type = "l", lwd = 2, col = "black",
main = "GKw Distribution (5 params)",
xlab = "x", ylab = "Density", las = 1, ylim = c(0, 3)
)
grid(col = "gray80", lty = 2)
# BKw (4 parameters)
plot(x, dbkw(x, 2, 3, 1.5, 2),
type = "l", lwd = 2, col = "darkgreen",
main = "BKw Distribution (4 params)",
xlab = "x", ylab = "Density", las = 1, ylim = c(0, 3)
)
grid(col = "gray80", lty = 2)
# EKw (3 parameters)
plot(x, dekw(x, 2, 3, 1.5),
type = "l", lwd = 2, col = "darkred",
main = "EKw Distribution (3 params)",
xlab = "x", ylab = "Density", las = 1, ylim = c(0, 3)
)
grid(col = "gray80", lty = 2)
# Kw (2 parameters) - base distribution
plot(x, dkw(x, 2, 3),
type = "l", lwd = 2, col = "darkblue",
main = "Kw Distribution (2 params)",
xlab = "x", ylab = "Density", las = 1, ylim = c(0, 3)
)
grid(col = "gray80", lty = 2)
## Example 3: Maximum Likelihood Estimation
# Generate data from Kumaraswamy distribution
set.seed(2024)
n <- 2000
true_alpha <- 2.5
true_beta <- 3.5
data <- rkw(n, true_alpha, true_beta)
# Obtain starting values
start_vals <- gkwgetstartvalues(data, family = "kw", n_starts = 3)
# Maximum likelihood estimation with analytical gradient
fit <- optim(
par = start_vals,
fn = llkw, # Negative log-likelihood
gr = grkw, # Analytical gradient
data = data,
method = "BFGS",
hessian = TRUE,
control = list(maxit = 500)
)
# Extract results
mle <- fit$par
se <- sqrt(diag(solve(fit$hessian)))
# Construct confidence intervals (95%)
ci <- data.frame(
Parameter = c("alpha", "beta"),
True = c(true_alpha, true_beta),
MLE = mle,
SE = se,
Lower = mle - 1.96 * se,
Upper = mle + 1.96 * se
)
print(ci, digits = 4)
# Goodness-of-fit diagnostic
x_grid <- seq(0.001, 0.999, length.out = 200)
fitted_dens <- dkw(x_grid, mle[1], mle[2])
true_dens <- dkw(x_grid, true_alpha, true_beta)
hist(data,
breaks = 40, probability = TRUE,
col = "lightgray", border = "white",
main = "Kumaraswamy Distribution Fit",
xlab = "Data", ylab = "Density", las = 1
)
lines(x_grid, fitted_dens, col = "red", lwd = 2, lty = 1)
lines(x_grid, true_dens, col = "blue", lwd = 2, lty = 2)
legend("topright",
legend = c("Data", "Fitted", "True"),
col = c("gray", "red", "blue"),
lwd = c(8, 2, 2), lty = c(1, 1, 2),
bty = "n"
)
grid(col = "gray80", lty = 2)
## Example 4: Model Selection Using Information Criteria
# Generate data from Exponentiated Kumaraswamy
set.seed(456)
n <- 1500
data <- rekw(n, alpha = 2, beta = 3, lambda = 1.5)
# Define competing models
models <- list(
Beta = list(
ll = function(par) llbeta(par, data),
gr = function(par) grbeta(par, data),
start = gkwgetstartvalues(data, family = "beta", n_starts = 2),
k = 2
),
Kw = list(
ll = function(par) llkw(par, data),
gr = function(par) grkw(par, data),
start = gkwgetstartvalues(data, family = "kw", n_starts = 2),
k = 2
),
EKw = list(
ll = function(par) llekw(par, data),
gr = function(par) grekw(par, data),
start = gkwgetstartvalues(data, family = "ekw", n_starts = 2),
k = 3
),
Mc = list(
ll = function(par) llmc(par, data),
gr = function(par) grmc(par, data),
start = gkwgetstartvalues(data, family = "mc", n_starts = 2),
k = 3
)
)
# Fit all models
results <- lapply(names(models), function(name) {
m <- models[[name]]
fit <- optim(par = m$start, fn = m$ll, gr = m$gr, method = "BFGS")
loglik <- -fit$value
aic <- -2 * loglik + 2 * m$k
bic <- -2 * loglik + m$k * log(n)
data.frame(
Model = name,
k = m$k,
LogLik = round(loglik, 2),
AIC = round(aic, 2),
BIC = round(bic, 2),
stringsAsFactors = FALSE
)
})
# Combine and sort by AIC
comparison <- do.call(rbind, results)
comparison <- comparison[order(comparison$AIC), ]
rownames(comparison) <- NULL
print(comparison)
cat("\nBest model by AIC:", comparison$Model[1], "\n")
cat("Best model by BIC:", comparison$Model[which.min(comparison$BIC)], "\n")
# }
Run the code above in your browser using DataLab