compositionalSpline {robCompositions}

Compositional spline

Description

This code implements the compositional smoothing splines grounded on the theory of Bayes spaces.

Usage

compositionalSpline(
  t,
  clrf,
  knots,
  w,
  order,
  der,
  alpha,
  spline.plot = FALSE,
  basis.plot = FALSE
)

Arguments

t

class midpoints

clrf

clr transformed values at class midpoints, i.e., fcenLR(f(t))

knots

sequence of knots

w

weights

order

order of the spline (i.e., degree + 1)

der

lth derivation

alpha

smoothing parameter

spline.plot

if TRUE, the resulting spline is plotted

basis.plot

if TRUE, the ZB-spline basis system is plotted

Details

The compositional splines enable to construct a spline basis in the centred logratio (clr) space of density functions (ZB-spline basis) and consequently also in the original space of densities (CB-spline basis).The resulting compositional splines in the clr space as well as the ZB-spline basis satisfy the zero integral constraint. This enables to work with compositional splines consistently in the framework of the Bayes space methodology.

Augmented knot sequence is obtained from the original knots by adding #(order-1) multiple endpoints.

Value

J

value of the functional J

ZB_coef

ZB-spline basis coeffcients

CV

score of cross-validation

GCV

score of generalized cross-validation

Author(s)

J. Machalova [e-mail], R. Talska [e-mail]

References

Machalova, J., Talska, R., Hron, K. Gaba, A. Compositional splines for representation of density functions. Comput Stat (2020). https://doi.org/10.1007/s00180-020-01042-7

Examples

# Example (Iris data):
SepalLengthCm <- iris$Sepal.Length
Species <- iris$Species
iris1 <- SepalLengthCm[iris$Species==levels(iris$Species)[1]]
h1 <- hist(iris1, plot = FALSE)
midx1 <- h1$mids
midy1 <- matrix(h1$density, nrow=1, ncol = length(h1$density), byrow=TRUE)
clrf  <- cenLR(rbind(midy1,midy1))$x.clr[1,]
knots <- seq(min(h1$breaks),max(h1$breaks),l=5)
order <- 4
der <- 2
alpha <- 0.99

sol1 <- compositionalSpline(t = midx1, clrf = clrf, knots = knots, 
  w = rep(1,length(midx1)), order = order, der = der, 
  alpha = alpha, spline.plot = TRUE)
sol1$GCV
ZB_coef <- sol1$ZB_coef
t <- seq(min(knots),max(knots),l=500)
t_step <- diff(t[1:2])
ZB_base <- ZBsplineBasis(t=t,knots,order)$ZBsplineBasis
sol1.t <- ZB_base%*%ZB_coef
sol2.t <- fcenLRinv(t,t_step,sol1.t)
h2 = hist(iris1,prob=TRUE,las=1)
points(midx1,midy1,pch=16)
lines(t,sol2.t,col="darkred",lwd=2)
# Example (normal distrubution):
# generate n values from normal distribution
set.seed(1)
n = 1000; mean = 0; sd = 1.5
raw_data = rnorm(n,mean,sd)
  
# number of classes according to Sturges rule
n.class = round(1+1.43*log(n),0)
  
# Interval midpoints
parnition = seq(-5,5,length=(n.class+1))
t.mid = c(); for (i in 1:n.class){t.mid[i]=(parnition[i+1]+parnition[i])/2}
  
counts = table(cut(raw_data,parnition))
prob = counts/sum(counts)                # probabilities
dens.raw = prob/diff(parnition)          # raw density data
clrf =  cenLR(rbind(dens.raw,dens.raw))$x.clr[1,]  # raw clr density data
  
# set the input parameters for smoothing 
knots = seq(min(parnition),max(parnition),l=5)
w = rep(1,length(clrf))
order = 4
der = 2
alpha = 0.5
spline = compositionalSpline(t = t.mid, clrf = clrf, knots = knots, 
  w = w, order = order, der = der, alpha = alpha, 
  spline.plot=TRUE, basis.plot=FALSE)
  
# ZB-spline coefficients
ZB_coef = spline$ZB_coef
  
# ZB-spline basis evaluated on the grid "t.fine"
t.fine = seq(min(knots),max(knots),l=1000)
ZB_base = ZBsplineBasis(t=t.fine,knots,order)$ZBsplineBasis
  
# Compositional spline in the clr space (evaluated on the grid t.fine)
comp.spline.clr = ZB_base%*%ZB_coef
  
# Compositional spline in the Bayes space (evaluated on the grid t.fine)
comp.spline = fcenLRinv(t.fine,diff(t.fine)[1:2],comp.spline.clr)
  
# Unit-integral representation of truncated true normal density function 
dens.true = dnorm(t.fine, mean, sd)/trapzc(diff(t.fine)[1:2],dnorm(t.fine, mean, sd))
  
# Plot of compositional spline together with raw density data
matplot(t.fine,comp.spline,type="l",
    lty=1, las=1, col="darkblue", xlab="t", 
    ylab="density",lwd=2,cex.axis=1.2,cex.lab=1.2,ylim=c(0,0.28))
matpoints(t.mid,dens.raw,pch = 8, col="darkblue", cex=1.3)
  
# Add true normal density function
matlines(t.fine,dens.true,col="darkred",lwd=2)