R Under development (unstable) (2023-11-23 r85618 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 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(SEL)
Loading required package: splines
Loading required package: quadprog
Loading required package: lattice
> ### Check calculation of the entropy term
> foo <- function(x, fit){
+   predict(fit, newdata = x)^2
+ }
> ### Case 1
> x <- c(2,5,7,9)
> y <- c(0.05, 0.2, 0.8, 0.9)
> fit1 <- SEL(x, y, gamma = 1, bounds = c(0,10))
> as.numeric(fit1$coefs%*%fit1$Omega%*%fit1$coefs)
[1] 0.1235645
> integrate(foo, 0, 10, fit=fit1)
0.1235647 with absolute error < 8.3e-05
> 
> ### Case 2
> x <- c(-0.9, -0.5, -0.1)
> y <- c(0.1, 0.8, 0.9)
> fit2 <- SEL(x, y, Delta = 0.03, bounds = c(-1,0), N = 10)
> as.numeric(fit2$coefs%*%fit2$Omega%*%fit2$coefs)
[1] 1.302755
> integrate(foo, -1, 0, fit=fit2)
1.302755 with absolute error < 9.1e-05
> 
> ### Case 3
> x <- seq(1e5, 2e5, length = 10)
> y <- c(seq(0.1, 0.4, length = 5), seq(0.65,0.9,length=5))
> fit3 <- SEL(x, y, bounds = c(5e4, 2.5e5), N=0, d=20, gamma = 2,
+             constr = "unimodal", mode = 1.5e5)
> as.numeric(fit3$coefs%*%fit3$Omega%*%fit3$coefs)
[1] 7.455394e-06
> integrate(foo, 5e4, 2.5e5, fit=fit3)
7.455394e-06 with absolute error < 1.8e-06
> 
> ### Case 4
> pistar <- function(x) dbeta(x, 1.5, 1.5)
> x <- c(0.1, 0.5, 0.9)
> y <- c(0.1, 0.8, 0.9)
> fit2 <- SEL(x, y, Delta = 0.03, bounds = c(0,1), N = 10, pistar=pistar)
> fit2
Fitted SEL object

Bounds:  0 1 
Inner Knots: 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 
Degree:  4 / Dimension:  15 
Shape Constraint: none 
gamma:  0.03168 / Delta 0.03 
Brier divergence to pistar:  0.3465014 
> sq <- seq(0, 1, length=2001)
> mean((predict(fit2, newdata=sq)-dbeta(sq, 1.5, 1.5))^2)
[1] 0.346443
> 
> proc.time()
   user  system elapsed 
   0.25    0.09    0.34