library("gamlss2")
## packages, code, and data
library("distributions3")
data("cars", package = "datasets")
## fit heteroscedastic normal GAMLSS model
## stopping distance (ft) explained by speed (mph)
<- gamlss2(dist ~ s(speed) | s(speed), data = cars, family = NO) m
GAMLSS-RS iteration 1: Global Deviance = 407.3647 eps = 0.125474
GAMLSS-RS iteration 2: Global Deviance = 405.766 eps = 0.003924
GAMLSS-RS iteration 3: Global Deviance = 405.7473 eps = 0.000045
GAMLSS-RS iteration 4: Global Deviance = 405.7473 eps = 0.000000
## obtain predicted distributions for three levels of speed
<- prodist(m, newdata = data.frame(speed = c(10, 20, 30)))
d print(d)
1 2
"GDF NO(mu = 23.02, sigma = 10.09)" "GDF NO(mu = 59.14, sigma = 18.49)"
3
"GDF NO(mu = 96.54, sigma = 33.89)"
## obtain quantiles (works the same for any distribution object 'd' !)
quantile(d, 0.5)
1 2 3
23.02058 59.13735 96.53656
quantile(d, c(0.05, 0.5, 0.95), elementwise = FALSE)
q_0.05 q_0.5 q_0.95
1 6.428687 23.02058 39.61248
2 28.721600 59.13735 89.55309
3 40.795319 96.53656 152.27779
quantile(d, c(0.05, 0.5, 0.95), elementwise = TRUE)
1 2 3
6.428687 59.137347 152.277794
## visualization
plot(dist ~ speed, data = cars)
<- data.frame(speed = 0:240/4)
nd $dist <- prodist(m, newdata = nd)
nd$fit <- quantile(nd$dist, c(0.05, 0.5, 0.95))
ndmatplot(nd$speed, nd$fit, type = "l", lty = 1, col = "slategray", add = TRUE)
## moments
mean(d)
1 2 3
23.02058 59.13735 96.53656
variance(d)
1 2 3
101.7507 341.9341 1148.4146
## simulate random numbers
random(d, 5)
r_1 r_2 r_3 r_4 r_5
1 28.76613 19.67624 21.26357 -9.481159 4.512701
2 56.42252 50.97136 68.48574 55.733796 35.136203
3 111.99632 114.53278 68.36612 137.467845 102.742254
## density and distribution
pdf(d, 50 * -2:2)
d_-100 d_-50 d_0 d_50 d_100
1 1.992414e-34 1.652166e-13 0.0029253352 0.001105976 8.928201e-15
2 1.783797e-18 5.885926e-10 0.0001297196 0.019094897 1.877372e-03
3 5.850207e-10 1.024848e-06 0.0002035755 0.004585323 1.171096e-02
cdf(d, 50 * -2:2)
p_-100 p_-50 p_0 p_50 p_100
1 1.637069e-34 2.260584e-13 0.0112397196 0.99625942 1.0000000
2 3.783013e-18 1.795168e-09 0.0006917068 0.31060410 0.9864409
3 3.324736e-09 7.657479e-06 0.0021951055 0.08483964 0.5407018
## Poisson example
data("FIFA2018", package = "distributions3")
<- gamlss2(goals ~ s(difference), data = FIFA2018, family = PO) m2
GAMLSS-RS iteration 1: Global Deviance = 355.3922 eps = 0.045332
GAMLSS-RS iteration 2: Global Deviance = 355.3922 eps = 0.000000
<- prodist(m2, newdata = data.frame(difference = 0))
d2 print(d2)
1
"GDF PO(mu = 1.237)"
quantile(d2, c(0.05, 0.5, 0.95))
[1] 0 1 3
## note that log_pdf() can replicate logLik() value
sum(log_pdf(prodist(m2), FIFA2018$goals))
[1] -177.6961
logLik(m2)
'log Lik.' -177.6961 (df=2.005144)