R Under development (unstable) (2024-11-18 r87347 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> library("mlt")
Loading required package: basefun
Loading required package: variables
> library("lattice")
> options(digits = 4)
> 
> data("faithful")
> 
> aic <- numeric(20)
> 
> for (o in 2:(length(aic) + 1)) {
+ Bs <- Bernstein_basis(numeric_var("waiting", support = range(faithful$waiting) + c(-5, 5)),
+                       order = o, ui = "incre")
+ m <- ctm(Bs)
+ mod <- mlt(m, data = faithful)
+ yp <- mkgrid(mod, 50)[["waiting"]]
+ 
+ aic[o - 1] <- AIC(mod)
+ 
+ pd <- data.frame(waiting = yp)
+ pd$p <- predict(mod, q = yp, type = "distribution")
+ 
+ plot(p ~ waiting, data = pd,
+      col = "red", pch = 21, main = paste("order", o, "aic", aic[o - 1]))
+ lines(ecdf(faithful$waiting))
+ 
+ }
> 
> plot(aic)
> 
> o <- which.min(aic) + 1
> Bs <- Bernstein_basis(numeric_var("waiting", support = range(faithful$waiting) + c(-5, 5)),
+                       order = o, ui = "incre")
> m <- ctm(Bs)
> mod <- mlt(m, data = faithful)
> 
> abline(h = AIC(mod))
> 
> pd$d <- predict(mod, q = yp, type = "density")
> 
> plot(d ~ waiting, data = pd, type = "l", col = "red", lwd = 3)
> lines(density(faithful$waiting))
> lines(rug(faithful$waiting))
> abline(h = 0)
> 
> p <- 1:99 / 100
> q <- predict(mod, p = p, K = 100, type = "quantile")
> 
> plot(p, q)
> lines(p, quantile(faithful$waiting, p))
> 
> Bs <- Bernstein_basis(numeric_var("waiting", support = range(faithful$waiting) + c(-5, 5)),
+                       order = o, ui = "incre")
> m <- ctm(Bs)
> mod <- mlt(m, data = faithful)
> 
> # H1 <- mod$optim(coef(mod), hessian = TRUE)$hessian
> H2 <- mod$hessian(coef(mod), weights(mod))
> 
> X <- model.matrix(m, faithful)
> Xprime <- model.matrix(m, faithful, deriv = c(waiting = 1))
> w <- drop((Xprime %*% coef(mod))^2)
> H3 <- crossprod(X) + crossprod(Xprime * w, Xprime)
> max(abs(H3 - H2))
[1] 9.886
> 
> cov2cor(vcov(mod))
             Bs1(waiting) Bs2(waiting) Bs3(waiting) Bs4(waiting) Bs5(waiting)
Bs1(waiting)       1.0000      -0.9653       0.9032      -0.8310       0.7514
Bs2(waiting)      -0.9653       1.0000      -0.9793       0.9310      -0.8628
Bs3(waiting)       0.9032      -0.9793       1.0000      -0.9835       0.9383
Bs4(waiting)      -0.8310       0.9310      -0.9835       1.0000      -0.9837
Bs5(waiting)       0.7514      -0.8628       0.9383      -0.9837       1.0000
Bs6(waiting)      -0.6650       0.7788      -0.8678       0.9358      -0.9824
Bs7(waiting)       0.5742      -0.6830       0.7767      -0.8583       0.9286
Bs8(waiting)      -0.4756       0.5730      -0.6626       0.7479      -0.8313
Bs9(waiting)       0.3573      -0.4341       0.5080      -0.5821       0.6607
             Bs6(waiting) Bs7(waiting) Bs8(waiting) Bs9(waiting)
Bs1(waiting)      -0.6650       0.5742      -0.4756       0.3573
Bs2(waiting)       0.7788      -0.6830       0.5730      -0.4341
Bs3(waiting)      -0.8678       0.7767      -0.6626       0.5080
Bs4(waiting)       0.9358      -0.8583       0.7479      -0.5821
Bs5(waiting)      -0.9824       0.9286      -0.8313       0.6607
Bs6(waiting)       1.0000      -0.9795       0.9080      -0.7427
Bs7(waiting)      -0.9795       1.0000      -0.9700       0.8292
Bs8(waiting)       0.9080      -0.9700       1.0000      -0.9214
Bs9(waiting)      -0.7427       0.8292      -0.9214       1.0000
> 
> if (FALSE) {
+ library("multcomp") ### since 1.0-3
+ 
+ mp <- parm(coef(mod), vcov(mod))
+ y <- mkgrid(mod, 30)$waiting
+ g <- glht(mp, linfct = model.matrix(mod$model,
+     data = data.frame(waiting = y)))
+ 
+ mc <- confint(g)
+ umc <- confint(g, calpha = qnorm(.975))
+ p <- mod$model$todistr$p
+ plot(y, p(mc$confint[, "Estimate"]), type = "l")
+ lines(y, p(mc$confint[, "lwr"]))
+ lines(y, p(mc$confint[, "upr"]))
+ lines(y, p(umc$confint[, "lwr"]))
+ lines(y, p(umc$confint[, "upr"]))
+ 
+ library("survival")
+ cm <- coxph(Surv(waiting, rep(TRUE, nrow(faithful))) ~ 1, data = faithful)
+ plot(survfit(cm))
+ lines(y, 1 - p(mc$confint[, "Estimate"]), col = "red")
+ lines(y, 1 - p(mc$confint[, "lwr"]), col = "red")
+ lines(y, 1 - p(mc$confint[, "upr"]), col = "red")
+ }
> 
> proc.time()
   user  system elapsed 
   4.98    0.79    5.76