We validate that the fast deviance evaluator matches the generic
gamlss.dist density route (up to floating-point tolerance)
across common families.
check_fast_vs_generic() on fitted
gamlss models and read the returned diagnostics.library(gamlss)
library(SelectBoost.gamlss)
set.seed(2025)
n <- 4000
families <- list(
NO = gamlss.dist::NO(),
PO = gamlss.dist::PO(),
LOGNO = gamlss.dist::LOGNO(),
GA = gamlss.dist::GA(),
IG = gamlss.dist::IG(),
LO = gamlss.dist::LO(),
LOGITNO = gamlss.dist::LOGITNO(),
GEOM = gamlss.dist::GEOM(),
BE = gamlss.dist::BE()
)
gen_fun <- list(
NO = function(n) gamlss.dist::rNO(n, mu = 0, sigma = 1),
PO = function(n) gamlss.dist::rPO(n, mu = 2.5),
LOGNO = function(n) gamlss.dist::rLOGNO(n, mu = 0, sigma = 0.6),
GA = function(n) gamlss.dist::rGA(n, mu = 2, sigma = 0.5),
IG = function(n) gamlss.dist::rIG(n, mu = 2, sigma = 0.5),
LO = function(n) gamlss.dist::rLO(n, mu = 0, sigma = 1),
LOGITNO = function(n) gamlss.dist::rLOGITNO(n, mu = 0, sigma = 1),
GEOM = function(n) gamlss.dist::rGEOM(n, mu = 3),
BE = function(n) gamlss.dist::rBE(n, mu = 0.4, sigma = 0.2)
)
summ <- lapply(names(families), function(fname) {
fam <- families[[fname]]; gen <- gen_fun[[fname]]
y <- gen(n); dat <- data.frame(y = y)
fit <- try(gamlss::gamlss(y ~ 1, data = dat, family = fam), silent = TRUE)
if (inherits(fit, "try-error")) return(NULL)
chk <- check_fast_vs_generic(fit, dat, tol = 1e-6)
data.frame(family = fname, ll_fast = chk$ll_fast, ll_generic = chk$ll_generic,
abs_diff = chk$abs_diff, pass = chk$pass)
})
#> GAMLSS-RS iteration 1: Global Deviance = 11375.33
#> GAMLSS-RS iteration 2: Global Deviance = 11375.33
#> GAMLSS-RS iteration 1: Global Deviance = 14782.83
#> GAMLSS-RS iteration 2: Global Deviance = 14782.83
#> GAMLSS-RS iteration 1: Global Deviance = 7267.968
#> GAMLSS-RS iteration 2: Global Deviance = 7267.968
#> GAMLSS-RS iteration 1: Global Deviance = 10699.55
#> GAMLSS-RS iteration 2: Global Deviance = 10699.55
#> GAMLSS-RS iteration 1: Global Deviance = 11694.15
#> GAMLSS-RS iteration 2: Global Deviance = 11694.15
#> GAMLSS-RS iteration 1: Global Deviance = 15982.17
#> GAMLSS-RS iteration 2: Global Deviance = 15982.16
#> GAMLSS-RS iteration 3: Global Deviance = 15982.16
#> GAMLSS-RS iteration 1: Global Deviance = 18087.96
#> GAMLSS-RS iteration 1: Global Deviance = -7154.868
#> GAMLSS-RS iteration 2: Global Deviance = -7422.771
#> GAMLSS-RS iteration 3: Global Deviance = -7422.823
#> GAMLSS-RS iteration 4: Global Deviance = -7422.823
summ <- do.call(rbind, Filter(Negate(is.null), summ))
summ
#> family ll_fast ll_generic abs_diff pass
#> 1 NO -5688.032 -5688.032 0 TRUE
#> 2 PO -10568.077 -10568.077 0 TRUE
#> 3 LOGNO -4398.621 -4398.621 0 TRUE
#> 4 GA -8009.049 -8009.049 0 TRUE
#> 5 IG -8144.880 -8144.880 0 TRUE
#> 6 LO -7991.389 -7991.389 0 TRUE
#> 7 GEOM -9326.121 -9326.121 0 TRUE
#> 8 BE 1678.217 1678.217 0 TRUEBelow we attempt an automatic sweep across many families; those without installed densities or generators will be skipped gracefully.
fams <- c("NO","PO","LOGNO","GA","IG","LO","LOGITNO","GEOM","BE",
"NBI","NBII","LOGLOG","DEL","ZAGA","ZIP","ZINBI","DPO","GPO",
"ZAIG","ZALG","BCT","BCPE","ZIPF","ZIPFmu","SICHEL","GLG",
"BETA4","RS","WEI","GIG")
n <- 3000
res <- lapply(fams, function(fam) {
y <- try(.gen_family(fam, n), silent = TRUE)
if (inherits(y, "try-error") || is.null(y)) return(NULL)
dat <- data.frame(y = y)
fam_obj <- try(getFromNamespace(fam, "gamlss.dist")(), silent = TRUE)
if (inherits(fam_obj, "try-error")) return(NULL)
fit <- try(gamlss::gamlss(y ~ 1, data = dat, family = fam_obj), silent = TRUE)
if (inherits(fit, "try-error")) return(NULL)
chk <- check_fast_vs_generic(fit, dat, tol = 1e-6)
data.frame(family = fam, abs_diff = chk$abs_diff, pass = chk$pass)
})
#> GAMLSS-RS iteration 1: Global Deviance = 8615.228
#> GAMLSS-RS iteration 2: Global Deviance = 8615.228
#> GAMLSS-RS iteration 1: Global Deviance = 10988.5
#> GAMLSS-RS iteration 2: Global Deviance = 10988.5
#> GAMLSS-RS iteration 1: Global Deviance = 5327.758
#> GAMLSS-RS iteration 2: Global Deviance = 5327.758
#> GAMLSS-RS iteration 1: Global Deviance = 7943.186
#> GAMLSS-RS iteration 2: Global Deviance = 7943.186
#> GAMLSS-RS iteration 1: Global Deviance = 8701.265
#> GAMLSS-RS iteration 2: Global Deviance = 8701.265
#> GAMLSS-RS iteration 1: Global Deviance = 11876.4
#> GAMLSS-RS iteration 2: Global Deviance = 11876.4
#> GAMLSS-RS iteration 3: Global Deviance = 11876.4
#> GAMLSS-RS iteration 1: Global Deviance = 13655.01
#> GAMLSS-RS iteration 1: Global Deviance = -5375.158
#> GAMLSS-RS iteration 2: Global Deviance = -5578.084
#> GAMLSS-RS iteration 3: Global Deviance = -5578.124
#> GAMLSS-RS iteration 4: Global Deviance = -5578.124
#> GAMLSS-RS iteration 1: Global Deviance = 11283.26
#> GAMLSS-RS iteration 2: Global Deviance = 11283.26
#> GAMLSS-RS iteration 1: Global Deviance = 11023
#> GAMLSS-RS iteration 2: Global Deviance = 11023
#> GAMLSS-RS iteration 3: Global Deviance = 11023
#> GAMLSS-RS iteration 1: Global Deviance = 10670
#> GAMLSS-RS iteration 2: Global Deviance = 10669.65
#> GAMLSS-RS iteration 3: Global Deviance = 10669.34
#> GAMLSS-RS iteration 4: Global Deviance = 10669.07
#> GAMLSS-RS iteration 5: Global Deviance = 10668.83
#> GAMLSS-RS iteration 6: Global Deviance = 10668.61
#> GAMLSS-RS iteration 7: Global Deviance = 10668.42
#> GAMLSS-RS iteration 8: Global Deviance = 10668.24
#> GAMLSS-RS iteration 9: Global Deviance = 10668.09
#> GAMLSS-RS iteration 10: Global Deviance = 10667.95
#> GAMLSS-RS iteration 11: Global Deviance = 10667.82
#> GAMLSS-RS iteration 12: Global Deviance = 10667.71
#> GAMLSS-RS iteration 13: Global Deviance = 10667.6
#> GAMLSS-RS iteration 14: Global Deviance = 10667.51
#> GAMLSS-RS iteration 15: Global Deviance = 10667.43
#> GAMLSS-RS iteration 16: Global Deviance = 10667.35
#> GAMLSS-RS iteration 17: Global Deviance = 10667.28
#> GAMLSS-RS iteration 18: Global Deviance = 10667.22
#> GAMLSS-RS iteration 19: Global Deviance = 10667.17
#> GAMLSS-RS iteration 20: Global Deviance = 10667.12
#> Warning in RS(): Algorithm RS has not yet converged
#> GAMLSS-RS iteration 1: Global Deviance = 9320.173
#> GAMLSS-RS iteration 2: Global Deviance = 9320.173
#> GAMLSS-RS iteration 1: Global Deviance = 10312.69
#> GAMLSS-RS iteration 2: Global Deviance = 10311.55
#> GAMLSS-RS iteration 3: Global Deviance = 10311.48
#> GAMLSS-RS iteration 4: Global Deviance = 10311.44
#> GAMLSS-RS iteration 5: Global Deviance = 10311.41
#> GAMLSS-RS iteration 6: Global Deviance = 10311.4
#> GAMLSS-RS iteration 7: Global Deviance = 10311.39
#> GAMLSS-RS iteration 8: Global Deviance = 10311.38
#> GAMLSS-RS iteration 9: Global Deviance = 10311.38
#> GAMLSS-RS iteration 10: Global Deviance = 10311.38
#> GAMLSS-RS iteration 11: Global Deviance = 10311.38
#> GAMLSS-RS iteration 12: Global Deviance = 10311.38
#> GAMLSS-RS iteration 1: Global Deviance = 11566.28
#> GAMLSS-RS iteration 2: Global Deviance = 11566.18
#> GAMLSS-RS iteration 3: Global Deviance = 11566.18
#> GAMLSS-RS iteration 1: Global Deviance = 12274.44
#> GAMLSS-RS iteration 2: Global Deviance = 12274.44
#> GAMLSS-RS iteration 1: Global Deviance = 9841.625
#> GAMLSS-RS iteration 2: Global Deviance = 9841.625
#> GAMLSS-RS iteration 1: Global Deviance = 12231.44
#> GAMLSS-RS iteration 2: Global Deviance = 12229.24
#> GAMLSS-RS iteration 3: Global Deviance = 12228.23
#> GAMLSS-RS iteration 4: Global Deviance = 12227.59
#> GAMLSS-RS iteration 5: Global Deviance = 12227.14
#> GAMLSS-RS iteration 6: Global Deviance = 12226.81
#> GAMLSS-RS iteration 7: Global Deviance = 12226.55
#> GAMLSS-RS iteration 8: Global Deviance = 12226.34
#> GAMLSS-RS iteration 9: Global Deviance = 12226.17
#> GAMLSS-RS iteration 10: Global Deviance = 12226.02
#> GAMLSS-RS iteration 11: Global Deviance = 12225.9
#> GAMLSS-RS iteration 12: Global Deviance = 12225.79
#> GAMLSS-RS iteration 13: Global Deviance = 12225.7
#> GAMLSS-RS iteration 14: Global Deviance = 12225.62
#> GAMLSS-RS iteration 15: Global Deviance = 12225.56
#> GAMLSS-RS iteration 16: Global Deviance = 12225.5
#> GAMLSS-RS iteration 17: Global Deviance = 12225.44
#> GAMLSS-RS iteration 18: Global Deviance = 12225.4
#> GAMLSS-RS iteration 19: Global Deviance = 12225.35
#> GAMLSS-RS iteration 20: Global Deviance = 12225.32
#> Warning in RS(): Algorithm RS has not yet converged
#> GAMLSS-RS iteration 1: Global Deviance = 10091.85
#> GAMLSS-RS iteration 2: Global Deviance = 10091.37
#> GAMLSS-RS iteration 3: Global Deviance = 10091.36
#> GAMLSS-RS iteration 4: Global Deviance = 10091.36
#> GAMLSS-RS iteration 1: Global Deviance = 9563.572
#> GAMLSS-RS iteration 2: Global Deviance = 9505.614
#> GAMLSS-RS iteration 3: Global Deviance = 9501.911
#> GAMLSS-RS iteration 4: Global Deviance = 9501.197
#> GAMLSS-RS iteration 5: Global Deviance = 9501.037
#> GAMLSS-RS iteration 6: Global Deviance = 9500.992
#> GAMLSS-RS iteration 7: Global Deviance = 9500.977
#> GAMLSS-RS iteration 8: Global Deviance = 9500.971
#> GAMLSS-RS iteration 9: Global Deviance = 9500.968
#> GAMLSS-RS iteration 10: Global Deviance = 9500.967
#> GAMLSS-RS iteration 11: Global Deviance = 9500.966
res <- do.call(rbind, Filter(Negate(is.null), res))
res
#> family abs_diff pass
#> 1 NO 0 TRUE
#> 2 PO 0 TRUE
#> 3 LOGNO 0 TRUE
#> 4 GA 0 TRUE
#> 5 IG 0 TRUE
#> 6 LO 0 TRUE
#> 7 GEOM 0 TRUE
#> 8 BE 0 TRUE
#> 9 NBI 0 TRUE
#> 10 NBII 0 TRUE
#> 11 DEL 0 TRUE
#> 12 ZAGA 0 TRUE
#> 13 ZINBI 0 TRUE
#> 14 DPO 0 TRUE
#> 15 GPO 0 TRUE
#> 16 ZAIG 0 TRUE
#> 17 SICHEL 0 TRUE
#> 18 WEI 0 TRUE
#> 19 GIG 0 TRUEfams <- c("NO","PO","LOGNO","GA","IG","LO","LOGITNO","GEOM","BE",
"NBI","NBII","LOGLOG","DEL","ZAGA","ZIP","ZINBI","DPO","GPO",
"ZAIG","ZALG","BCT","BCPE","ZIPF","ZIPFmu","SICHEL","GLG",
"BETA4","RS","WEI","GIG")
n <- 3000
tols <- SelectBoost.gamlss:::.family_tolerance()
tol_default <- tols[['.default']]
rows <- lapply(fams, function(fam) {
# try generator
y <- try(.gen_family(fam, n), silent = TRUE)
if (inherits(y, "try-error") || is.null(y)) {
return(data.frame(family=fam, status="skip", reason="generator missing/failed", abs_diff=NA_real_, pass=NA))
}
dat <- data.frame(y = y)
# try family object
fam_obj <- try(getFromNamespace(fam, "gamlss.dist")(), silent = TRUE)
if (inherits(fam_obj, "try-error")) {
return(data.frame(family=fam, status="skip", reason="family constructor missing", abs_diff=NA_real_, pass=NA))
}
# try fit
fit <- try(gamlss::gamlss(y ~ 1, data = dat, family = fam_obj), silent = TRUE)
if (inherits(fit, "try-error")) {
return(data.frame(family=fam, status="skip", reason="fit failed", abs_diff=NA_real_, pass=NA))
}
# evaluate
tol_fam <- tols[[fam]] %||% tol_default
chk <- check_fast_vs_generic(fit, dat, tol = tol_fam)
data.frame(family=fam, status="ok", reason="", abs_diff=chk$abs_diff, pass=chk$pass)
})
#> GAMLSS-RS iteration 1: Global Deviance = 8484.742
#> GAMLSS-RS iteration 2: Global Deviance = 8484.742
#> GAMLSS-RS iteration 1: Global Deviance = 10872.98
#> GAMLSS-RS iteration 2: Global Deviance = 10872.98
#> GAMLSS-RS iteration 1: Global Deviance = 5341.986
#> GAMLSS-RS iteration 2: Global Deviance = 5341.986
#> GAMLSS-RS iteration 1: Global Deviance = 8015.778
#> GAMLSS-RS iteration 2: Global Deviance = 8015.778
#> GAMLSS-RS iteration 1: Global Deviance = 8692.182
#> GAMLSS-RS iteration 2: Global Deviance = 8692.182
#> GAMLSS-RS iteration 1: Global Deviance = 11909.62
#> GAMLSS-RS iteration 2: Global Deviance = 11909.55
#> GAMLSS-RS iteration 3: Global Deviance = 11909.55
#> GAMLSS-RS iteration 1: Global Deviance = 13515
#> GAMLSS-RS iteration 1: Global Deviance = -5382.726
#> GAMLSS-RS iteration 2: Global Deviance = -5583.397
#> GAMLSS-RS iteration 3: Global Deviance = -5583.435
#> GAMLSS-RS iteration 4: Global Deviance = -5583.435
#> GAMLSS-RS iteration 1: Global Deviance = 11308.87
#> GAMLSS-RS iteration 2: Global Deviance = 11308.87
#> GAMLSS-RS iteration 1: Global Deviance = 11075.33
#> GAMLSS-RS iteration 2: Global Deviance = 11075.33
#> GAMLSS-RS iteration 3: Global Deviance = 11075.33
#> GAMLSS-RS iteration 1: Global Deviance = 10735.22
#> GAMLSS-RS iteration 2: Global Deviance = 10735.04
#> GAMLSS-RS iteration 3: Global Deviance = 10734.89
#> GAMLSS-RS iteration 4: Global Deviance = 10734.76
#> GAMLSS-RS iteration 5: Global Deviance = 10734.64
#> GAMLSS-RS iteration 6: Global Deviance = 10734.54
#> GAMLSS-RS iteration 7: Global Deviance = 10734.46
#> GAMLSS-RS iteration 8: Global Deviance = 10734.38
#> GAMLSS-RS iteration 9: Global Deviance = 10734.32
#> GAMLSS-RS iteration 10: Global Deviance = 10734.27
#> GAMLSS-RS iteration 11: Global Deviance = 10734.22
#> GAMLSS-RS iteration 12: Global Deviance = 10734.19
#> GAMLSS-RS iteration 13: Global Deviance = 10734.15
#> GAMLSS-RS iteration 14: Global Deviance = 10734.13
#> GAMLSS-RS iteration 15: Global Deviance = 10734.11
#> GAMLSS-RS iteration 16: Global Deviance = 10734.09
#> GAMLSS-RS iteration 17: Global Deviance = 10734.08
#> GAMLSS-RS iteration 18: Global Deviance = 10734.08
#> GAMLSS-RS iteration 19: Global Deviance = 10734.07
#> GAMLSS-RS iteration 20: Global Deviance = 10734.07
#> GAMLSS-RS iteration 1: Global Deviance = 9155.055
#> GAMLSS-RS iteration 2: Global Deviance = 9155.055
#> GAMLSS-RS iteration 1: Global Deviance = 10422.64
#> GAMLSS-RS iteration 2: Global Deviance = 10421.36
#> GAMLSS-RS iteration 3: Global Deviance = 10421.2
#> GAMLSS-RS iteration 4: Global Deviance = 10421.1
#> GAMLSS-RS iteration 5: Global Deviance = 10421.05
#> GAMLSS-RS iteration 6: Global Deviance = 10421.01
#> GAMLSS-RS iteration 7: Global Deviance = 10420.99
#> GAMLSS-RS iteration 8: Global Deviance = 10420.98
#> GAMLSS-RS iteration 9: Global Deviance = 10420.97
#> GAMLSS-RS iteration 10: Global Deviance = 10420.97
#> GAMLSS-RS iteration 11: Global Deviance = 10420.97
#> GAMLSS-RS iteration 12: Global Deviance = 10420.96
#> GAMLSS-RS iteration 13: Global Deviance = 10420.96
#> GAMLSS-RS iteration 1: Global Deviance = 11345.71
#> GAMLSS-RS iteration 2: Global Deviance = 11345.63
#> GAMLSS-RS iteration 3: Global Deviance = 11345.63
#> GAMLSS-RS iteration 1: Global Deviance = 12225.98
#> GAMLSS-RS iteration 2: Global Deviance = 12225.98
#> GAMLSS-RS iteration 1: Global Deviance = 9806.759
#> GAMLSS-RS iteration 2: Global Deviance = 9806.759
#> GAMLSS-RS iteration 1: Global Deviance = 12050.13
#> GAMLSS-RS iteration 2: Global Deviance = 12049.42
#> GAMLSS-RS iteration 3: Global Deviance = 12049.14
#> GAMLSS-RS iteration 4: Global Deviance = 12048.99
#> GAMLSS-RS iteration 5: Global Deviance = 12048.9
#> GAMLSS-RS iteration 6: Global Deviance = 12048.84
#> GAMLSS-RS iteration 7: Global Deviance = 12048.8
#> GAMLSS-RS iteration 8: Global Deviance = 12048.77
#> GAMLSS-RS iteration 9: Global Deviance = 12048.75
#> GAMLSS-RS iteration 10: Global Deviance = 12048.74
#> GAMLSS-RS iteration 11: Global Deviance = 12048.72
#> GAMLSS-RS iteration 12: Global Deviance = 12048.71
#> GAMLSS-RS iteration 13: Global Deviance = 12048.71
#> GAMLSS-RS iteration 14: Global Deviance = 12048.7
#> GAMLSS-RS iteration 15: Global Deviance = 12048.7
#> GAMLSS-RS iteration 16: Global Deviance = 12048.69
#> GAMLSS-RS iteration 17: Global Deviance = 12048.69
#> GAMLSS-RS iteration 18: Global Deviance = 12048.69
#> GAMLSS-RS iteration 19: Global Deviance = 12048.69
#> GAMLSS-RS iteration 20: Global Deviance = 12048.68
#> Warning in RS(): Algorithm RS has not yet converged
#> GAMLSS-RS iteration 1: Global Deviance = 10251.6
#> GAMLSS-RS iteration 2: Global Deviance = 10251.24
#> GAMLSS-RS iteration 3: Global Deviance = 10251.24
#> GAMLSS-RS iteration 4: Global Deviance = 10251.24
#> GAMLSS-RS iteration 1: Global Deviance = 9539.068
#> GAMLSS-RS iteration 2: Global Deviance = 9504.096
#> GAMLSS-RS iteration 3: Global Deviance = 9502.578
#> GAMLSS-RS iteration 4: Global Deviance = 9502.384
#> GAMLSS-RS iteration 5: Global Deviance = 9502.352
#> GAMLSS-RS iteration 6: Global Deviance = 9502.346
#> GAMLSS-RS iteration 7: Global Deviance = 9502.344
#> GAMLSS-RS iteration 8: Global Deviance = 9502.343
res <- do.call(rbind, rows)
res
#> family status reason abs_diff pass
#> 1 NO ok 0 TRUE
#> 2 PO ok 0 TRUE
#> 3 LOGNO ok 0 TRUE
#> 4 GA ok 0 TRUE
#> 5 IG ok 0 TRUE
#> 6 LO ok 0 TRUE
#> 7 LOGITNO skip fit failed NA NA
#> 8 GEOM ok 0 TRUE
#> 9 BE ok 0 TRUE
#> 10 NBI ok 0 TRUE
#> 11 NBII ok 0 TRUE
#> 12 LOGLOG skip generator missing/failed NA NA
#> 13 DEL ok 0 TRUE
#> 14 ZAGA ok 0 TRUE
#> 15 ZIP skip generator missing/failed NA NA
#> 16 ZINBI ok 0 TRUE
#> 17 DPO ok 0 TRUE
#> 18 GPO ok 0 TRUE
#> 19 ZAIG ok 0 TRUE
#> 20 ZALG skip generator missing/failed NA NA
#> 21 BCT skip generator missing/failed NA NA
#> 22 BCPE skip generator missing/failed NA NA
#> 23 ZIPF skip generator missing/failed NA NA
#> 24 ZIPFmu skip generator missing/failed NA NA
#> 25 SICHEL ok 0 TRUE
#> 26 GLG skip generator missing/failed NA NA
#> 27 BETA4 skip generator missing/failed NA NA
#> 28 RS skip generator missing/failed NA NA
#> 29 WEI ok 0 TRUE
#> 30 GIG ok 0 TRUE