The family one distribution, seven models

Every model below is GKw with parameters held fixed. Follow an edge to fix a parameter and drop to a simpler model; the number in the badge is how many free parameters remain.

Uniform is GKw(1, 1, 1, 0, 1) — the neutral value of every parameter. Kw and Beta are not nested in one another, and neither are EKw/Mc or BKw/KKw.

What each parameter does

49 functions, one naming rule

prefix + family code. Row = your model, column = the job. blue distribution API   red likelihood API.

Mind the underscore. Only the Beta row is irregular: d/p/q/r keep the trailing _ so they do not mask stats::dbeta — the likelihood trio does not. llbeta, not llbeta_.

Distribution functions d · p · q · r

dgkw(x, alpha, beta, gamma, delta, lambda, log = FALSE)
pgkw(q, ...,  lower.tail = TRUE, log.p = FALSE)
qgkw(p, ...,  lower.tail = TRUE, log.p = FALSE)
rgkw(n, ...)

A sub-family drops the parameters it fixes: dkw(x, alpha, beta, log = FALSE).

The base-R contract these keep

• Recycling. x and every parameter are recycled to a common length: dkw(c(.3,.6), c(1,2,3,4), 2) gives 4 values.
• Attributes. dim, dimnames and names of the first argument survive when the lengths agree, as in stats::dbeta.
• Outside [0,1] the density is 0; at the two endpoints it is the limit (box, left).
• An invalid parameter is an error — α ≤ 0, δ < 0, NA all stop(). They do not return NaN.
• Defaults are the neutral values, so dgkw(x) is the uniform density.

Shapes the family reaches

Real output. The boundary rule (column 1) predicts each end.

What it will not do. At most one interior mode — a second peak is always a divergence at 0 or 1, as in the bottom row. Two interior humps need a mixture, not a bigger GKw.

Pick a starting family

Fit up, not down. Add a parameter only when a likelihood-ratio test (page 2) says it earns its place.

Fit by maximum likelihood

Objective, gradient and Hessian drop straight into stats::optim.

x <- rekw(2000, 2, 3, 1.5)        # data on (0,1)

start <- gkwgetstartvalues(x, family = "ekw")

fit <- optim(start,
             fn     = llekw,   # −logL
             gr     = grekw,   # −score
             data   = x,       # passed through
             method = "BFGS",
             hessian = TRUE)

est <- fit$par                    # MLE
logLik <- -fit$value              # note the sign
se  <- sqrt(diag(solve(fit$hessian)))
ci  <- est + outer(se, c(-1.96, 1.96))

fn is already the negative log-likelihood, so optim's default minimisation is the right direction and fit$hessian is the observed information — invert it directly, no sign flip.

For the exact information matrix, use the analytical Hessian:

se <- sqrt(diag(solve( hsekw(est, x) )))

They agree to ~5e−7 relative: optim's is a finite difference.

Starting values

gkwgetstartvalues(x, family = "gkw", n_starts = 5)

Method of moments by Nelder–Mead, returning a named vector already in the par order that family's ll* expects. family is one of "gkw" "bkw" "kkw" "ekw" "mc" "kw" "beta" (case-insensitive).

• Deterministic. set.seed() has no effect and it never touches .Random.seed. Widen the search with n_starts; there is no seed.
• Four fixed starting points are always used, so n_starts < 4 behaves as 4; cost grows linearly and the objective is non-increasing in n_starts.
• Results are clipped to a box — α, β ∈ (0.1, 50), γ ∈ (0.1, 10) but (0.1, 50) for "beta", δ ∈ (0.01, 10), λ ∈ (0.1, 20) — so a value sitting exactly on a bound is a hint that the box bit, not an estimate.
• Values outside (0,1) are truncated with a warning. Treat that warning as a scale problem in the data — a 0–100 percentage — not as noise.

Return values under failure

An infinite objective gives the optimiser no direction — it sits there reporting convergence. Screen first: stopifnot(all(x > 0 & x < 1)). An exact 0 or 1 is the usual culprit; squeeze it in with (x*(n-1) + 0.5)/n.

Likelihood functions ll · gr · hs

llgkw(par, data)   → scalar   −ℓ(θ)
grgkw(par, data)   → vector   −∇ℓ(θ)
hsgkw(par, data)   → matrix   −∇²ℓ(θ)

All three carry a minus sign: minimise ll*, and hs*(est, x) is the observed information. par is a plain numeric vector — the names are ignored, only the order matters.

Why analytical derivatives

One pass over the data instead of a finite-difference stencil per parameter. Measured at n = 20 000, GKw, against numDeriv:

Agreement with numDeriv to ~9e−6 absolute — the difference is the finite-difference truncation error, not the analytical form.

Fit all seven, then rank

fam <- c("gkw","bkw","kkw","ekw","mc","kw","beta")
ll  <- list(llgkw, llbkw, llkkw, llekw, llmc, llkw, llbeta)
gr  <- list(grgkw, grbkw, grkkw, grekw, grmc, grkw, grbeta)

tab <- Map(function(f, l, g) {
  s <- gkwgetstartvalues(x, family = f)
  o <- optim(s, l, g, data = x, method = "BFGS")
  c(k = length(s), AIC = 2*o$value + 2*length(s))
}, fam, ll, gr)

tab <- do.call(rbind, tab); tab[order(tab[, "AIC"]), ]

The three lists run in the same order as page 1's matrix, and gkwgetstartvalues() returns each par already in the order its family expects — nothing is aligned by hand.

Compare nested models

Comparable by likelihood-ratio test only if one is the other with parameters fixed; df counts them.

lrt <- 2 * (fit_red$value - fit_full$value)
pchisq(lrt, df, lower.tail = FALSE)

† δ = 0 is on the boundary, so χ²(df) is conservative: the null is 0.5·χ²(df−1) + 0.5·χ²(df).

Not nested: kw vs beta, ekw vs mc, bkw vs kkw. Compare those with AIC/BIC only.

Information criteria

k <- length(fit$par);  ll <- -fit$value
AIC <- -2*ll + 2*k
BIC <- -2*ll + k*log(length(x))

Lower is better; a gap under ~2 is not evidence. All seven are fitted to the same x, so they compare directly.

Check the fit

# PIT residuals — uniform if the model holds
u <- pekw(x, est[1], est[2], est[3])
ks.test(u, "punif")

# Q-Q against the fitted quantiles
q <- qekw(ppoints(length(x)), est[1], est[2], est[3])
qqplot(q, sort(x)); abline(0, 1, col = 2)

Estimating from the same x makes the KS p-value optimistic.