Nothing
cub_spline_est <- function(xtab, ytab, x, kn = ceiling((length(xtab))^(1/4)),
method = "u", all.dea = FALSE, control = list("tm_limit" = 700)){
# verification
stopifnot(method%in%c("u","m","mc"), length(xtab)==length(ytab))
# initialisation
cv = switch(method, m = 0, mc = 1, u = -1)
spl_ord <- 4 # cubic spline
# support = [l_end,r_end]
l_end <- min(x)
r_end <- max(x)
if(method=="u"){
dea_fdh_index <- 1:(length(xtab))
}else
{dea_fdh_index <- which(dea(xtab, ytab, RTS = cv, ORIENTATION = "out")$eff==1)
}
# change the number ok knots to select all the DEA enveloppe points
lgrid <- (r_end-l_end)/(kn+1)
if(!all.dea){
knots <- c(l_end, quantile(xtab[dea_fdh_index], prob = 1:kn/(kn+1), type = 7), r_end)
}else
{knots <- c(l_end, sort(xtab[dea_fdh_index]), r_end)
}
rind <- which(diff(knots)<0.001)+1 # remove redundant knots
if(length(rind>0))
knots <- knots[-rind]
kn <- length(knots)-2 # corrected number of internal knots
knots <- c(-3*lgrid+l_end, -2*lgrid+l_end, -lgrid+l_end, knots,
r_end+lgrid, r_end+2*lgrid, r_end+3*lgrid)
x_cal <- seq(l_end, r_end, length.out = 5001)
der <- rep(0, length(x_cal))
bb <- splines::spline.des(knots, ord = spl_ord, derivs = der, x = x_cal, outer.ok = TRUE)
opt_coef <- apply(bb$design, 2, sum)/(length(x_cal)-1)
# design and constraint matrix
desM <- splines::spline.des(knots, ord = spl_ord, x = xtab, outer.ok = TRUE)$design
nbasis <- kn + spl_ord
monoC <- cbind(-diag(nbasis-1), 0) + cbind(0, diag(nbasis-1))
if(method=="u"){
mat <- desM
rhs <- ytab
dir <- rep(">=", length(xtab))
}
if(cv==0){
# formulation of linear programming
mat <- rbind(desM, monoC)
rhs <- c(ytab, rep(0, nbasis-1))
dir <- c(rep(">=", length(xtab)), rep(">=", nbasis-1))
}
if(cv==1){
# concavity constraint
nknots <- knots[c(spl_ord:(length(knots)-(spl_ord-1)))]
der_n <- rep(2, length(nknots))
convC <- spline.des(knots, ord = spl_ord, derivs = der_n, x = nknots, outer.ok = TRUE)$design
# formulation of linear programming
mat <- rbind(desM, monoC, convC)
rhs <- c(ytab, rep(0, nbasis-1), rep(0, length(nknots)))
dir <- c(rep(">=", length(xtab)), rep(">=", nbasis-1), rep("<=", length(nknots)))
}
bounds <- list(lower = list(ind = 1:length(opt_coef), val = rep(-Inf, length(opt_coef))),
upper = list(ind = 1:length(opt_coef), val = rep(Inf, length(opt_coef))))
Sol <- Rglpk_solve_LP(opt_coef, mat, dir, rhs, bounds, types = NULL, max = FALSE, control = control)
if(Sol$status!=0){
# warning("It seems no optimal solution has been found by Glpk")
OPT <- rep(NA, length(opt_coef))
}else
OPT <- Sol$solution
xgrid <- x
der <- rep(0, length(xgrid))
fitt <- spline.des(knots, ord = spl_ord, derivs = der, x = xgrid, outer.ok = TRUE)$design %*% OPT
return(fitt)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.