Nothing
#' @importFrom ggridges geom_density_ridges
NULL
#' A function for summarize HTS
#' @description A summarizer for HTS
#' @param x a HTS
#' @return A matrix with basic statistics.
#' @examples
#' summaryHTS(subsetHTS(RetHTS, from = 1, to = 10))
#' @export
summaryHTS <- function(x) {
df <- data.frame(
Tstamp = numeric(),
T_period_start = numeric(),
T_period_end = numeric(),
mean = numeric(),
std = numeric(),
skew = numeric(),
kurt = numeric(),
min = numeric(),
Q1 = numeric(),
MED = numeric(),
Q3 = numeric(),
Max = numeric()
)
nr <- length(x@data)
for (i in 1:nr) {
tmpo <- x@data[[i]]
newrow <- data.frame(
Tstamp = tmpo@tstamp,
T_period_start = tmpo@period$start,
T_period_end = tmpo@period$end,
mean = tmpo@m,
std = tmpo@s,
skew = skewH(tmpo),
kurt = kurtH(tmpo),
min = tmpo@x[1],
Q1 = compQ(tmpo, 0.25),
MED = compQ(tmpo, 0.5),
Q3 = compQ(tmpo, 0.75),
Max = tmpo@x[length(tmpo@x)]
)
df <- rbind(df, newrow)
}
return(df)
}
# OK assigned
plot.gg <- function(x, type = "HISTO", col = "green", border = "black") {
xlabel <- paste("m= ", format(x@m, nsmall = 4), " std= ", format(x@s, nsmall = 4))
if (type == "HISTO") {
lowers <- x@x[1:(length(x@x) - 1)]
uppers <- x@x[2:length(x@x)]
maxdens <- 10 / (max(uppers) - lowers[1])
ampl <- uppers - lowers
# ampl[which(ampl==0)]=1
dens <- (x@p[2:length(x@p)] - x@p[1:(length(x@p) - 1)]) / ampl
dens[which(dens > maxdens)] <- maxdens
df <- data.frame(xm = lowers, xM = uppers, ym = dens * 0, yM = dens)
p <- with(df, ggplot(df, aes(xmin = xm, xmax = xM, ymin = ym, ymax = yM)) +
geom_rect(alpha = 0.7, fill = col, colour = border) +
geom_vline(xintercept = x@m, colour = col, linetype = "dashed", size = 0.5) +
xlab(xlabel) +
ylab("density") +
ggtitle("Histogram"))
}
if (type == "CDF") {
xs <- x@x
ps <- x@p
df <- data.frame(x = xs, cdf = ps)
p <- with(df, ggplot(df, aes(x = x, y = cdf)) +
geom_line(colour = border) +
xlab("x") +
ylab("p") +
ggtitle("Cumulative Distribution Function"))
}
if (type == "QF") {
xs <- x@x
ps <- x@p
df <- data.frame(x = xs, cdf = ps)
p <- with(df, ggplot(df, aes(x = cdf, y = x)) +
geom_line(colour = border) +
xlab("p") +
ylab("x") +
ggtitle("Quantile Function"))
}
if (type == "DENS") {
# generate 200 random points according to the QF
rn <- 200
xn <- c(rep(0, rn + 1))
random_no <- c(0:rn) / rn
for (i in 0:rn) {
xn[i + 1] <- compQ(x, random_no[i + 1])
}
df <- data.frame(x = xn)
p <- with(df, ggplot(df, aes(x = x)) +
geom_density(alpha = 0.7, fill = col, colour = border) +
xlab(xlabel) +
ylab("density") +
ggtitle("Density plot (KDE)"))
p <- p + geom_vline(xintercept = x@m, colour = "black", linetype = "dashed", size = 0.5)
}
if (type == "HBOXPLOT") {
qua <- c(0, 0.25, 0.5, 0.75, 1)
df <- data.frame(ymin = 0, lower = 0, middle = 0, upper = 0, ymax = 0)
for (i in 1:5) {
df[[i]] <- compQ(x, qua[i])
}
df$x <- factor(0)
p <- with(df, ggplot(df, aes_all(names(df))) +
geom_boxplot(stat = "identity", fill = col, size = 1, colour = border) +
ylab(xlabel) +
xlab("") +
ggtitle("Horizontal Boxplot") +
coord_flip())
# xn=c(0,0,0,0,0)
# for (i in 1:5){
# xn[i]=compQ(x,qua[i])
# }
# df=data.frame(x=xn)
# p=with(df,ggplot(df, aes(x=factor(0), ymin = x[1], lower = x[2], middle =x[3], upper =x[4], ymax =x[5]))+
# geom_boxplot(stat = "identity", fill=col, size=1, colour=border)+
# ylab(xlabel) +xlab("")+
# ggtitle("Horizontal Boxplot")+coord_flip()
# )
}
if (type == "VBOXPLOT") {
qua <- c(0, 0.25, 0.5, 0.75, 1)
df <- data.frame(ymin = 0, lower = 0, middle = 0, upper = 0, ymax = 0)
for (i in 1:5) {
df[[i]] <- compQ(x, qua[i])
}
df$x <- factor(0)
p <- with(df, ggplot(df, aes_all(names(df))) +
geom_boxplot(stat = "identity", fill = col, size = 1, colour = border) +
ylab("x") +
xlab(xlabel) +
ggtitle("Vertical Boxplot"))
# p=with(df,ggplot(df, aes(x=factor(0), ymin = x[1], lower = x[2], middle =x[3], upper =x[4], ymax =x[5]))+
# geom_boxplot(stat = "identity", fill=col, size=1, colour=border)+
# ylab("x") +xlab(xlabel)+
# ggtitle("Vertical Boxplot")
# )
}
return(p)
}
WH.joy <- function(DATA, var) {
n <- get.MatH.nrows(DATA)
n_names <- get.MatH.rownames(DATA)
vp <- c(0:1000) / 1000
vec1 <- numeric()
labs <- character()
for (i in 1:n) {
tmp <- COMP_Q_VECT(DATA@M[i, var][[1]]@x, DATA@M[i, var][[1]]@p, vp)
vec1 <- c(vec1, tmp)
labs <- c(labs, rep(n_names[i], length(tmp)))
}
vname <- rep(get.MatH.varnames(DATA)[var], length(labs))
DF <- data.frame(x = vec1, y = as.factor(labs), z = vname)
return(DF)
}
WH.joymult <- function(DATA, list_of_vars) { ## AS SUBSTUTE OF PLOT MATH DENS
DF <- WH.joy(DATA, list_of_vars[1])
if (length(list_of_vars) > 1) {
for (i in 2:length(list_of_vars)) {
tmp <- WH.joy(DATA, list_of_vars[i])
DF <- rbind(DF, tmp)
}
}
DF$y <- factor(DF$y, levels = rev(get.MatH.rownames(DATA)))
# levels(DF$y)=rev(levels(DF$y))
return(DF)
}
WH.joy.HIST <- function(DATA, var, tr = 1) {
n <- get.MatH.nrows(DATA)
n_names <- get.MatH.rownames(DATA)
vec1 <- numeric()
vec2 <- numeric()
vec3 <- numeric()
labs <- character()
for (i in 1:n) {
dom <- DATA@M[i, var][[1]]@x
widths <- diff(dom)
widths[which(widths == 0)] <- 1e-16
p <- diff(DATA@M[i, var][[1]]@p)
densit <- p / widths
tmpx <- numeric() # sort(c(dom,dom))
tmpy <- numeric()
for (j in 1:length(densit)) {
tmpx <- c(tmpx, dom[j], dom[j], dom[j + 1], dom[j + 1])
tmpy <- c(tmpy, 0, densit[j], densit[j], 0)
}
# tmpy=c(tmpy,0)
vec1 <- c(vec1, tmpx)
vec2 <- c(vec2, tmpy)
vec3 <- c(vec3, rep(i, length(tmpx)))
labs <- c(labs, rep(n_names[i], length(tmpx)))
}
vname <- rep(get.MatH.varnames(DATA)[var], length(labs))
if (tr == 1) {
treshold <- 2 # quantile(vec2, probs = 0.8)
} else {
treshold <- 1.5 * diff(quantile(vec2, probs = c(0.25, 0.75))) + quantile(vec2, probs = c(0.75))
}
vec2[vec2 > treshold] <- treshold
DF <- data.frame(X = vec1, Y = vec3, H = vec2, y = labs, z = as.factor(vname))
return(DF)
}
WH.joymult.HISTO <- function(DATA, list_of_vars, tr = 1) {
DF <- WH.joy.HIST(DATA, list_of_vars[1], tr = tr)
if (length(list_of_vars) > 1) {
for (i in 2:length(list_of_vars)) {
tmp <- WH.joy.HIST(DATA, list_of_vars[i], tr = tr)
DF <- rbind(DF, tmp)
}
}
DF$y <- factor(DF$y, levels = rev(get.MatH.rownames(DATA)))
return(DF)
}
# OK assigned
plot.M <- function(x, type = "HISTO", border = "black", angL = 330) {
varsno <- ncol(x@M)
indno <- nrow(x@M)
if (varsno == 1 && indno == 1) {
tmpo <- x@M[1, 1][[1]]
p <- plot.gg(tmpo, type = type, col = "green", border = border)
}
else {
if (type == "HISTO") {
df <- WH.joymult.HISTO(DATA = x, list_of_vars = c(1:get.MatH.ncols(x)), tr = 1)
CC <- get.MatH.ncols(x)
p <- with(
df,
ggplot(df, aes(x = X, y = y, height = H, group = y, fill = z)) +
geom_density_ridges(stat = "identity", scale = 1.1, alpha = 0.7, color = border) +
facet_wrap(~z, scales = "free_x", ncol = CC) +
theme(
strip.text.x = element_text(face = "bold"),
axis.text.x = element_text(size = 6),
axis.text.y = element_text(face = "bold", vjust = 0),
legend.position = "none",
axis.title.x = element_blank(), axis.title.y = element_blank()
)
)
}
if (type == "DENS") {
df <- WH.joymult(DATA = x, list_of_vars = c(1:get.MatH.ncols(x)))
CC <- get.MatH.ncols(x)
p <- with(
df,
ggplot(df, aes(x = x, y = y, fill = z)) +
geom_density_ridges(scale = 1.3, rel_min_height = 0.01, alpha = 0.7) +
facet_wrap(~z, scales = "free_x", ncol = CC) +
theme(
strip.text.x = element_text(face = "bold"),
axis.text.x = element_text(size = 6),
axis.text.y = element_text(face = "bold", vjust = 0),
legend.position = "none",
axis.title.x = element_blank(), axis.title.y = element_blank()
)
)
}
if (type == "BOXPLOT") {
df <- data.frame(
ymin = numeric(),
yQ1 = numeric(),
yME = numeric(),
yQ3 = numeric(),
ymax = numeric(),
indiv = character(), vars = character(), stringsAsFactors = FALSE
)
for (i in 1:indno) {
for (j in 1:varsno) {
qua <- c(0, 0.25, 0.5, 0.75, 1)
xn <- c(0, 0, 0, 0, 0)
tmpo <- x@M[i, j][[1]]
for (k in 1:5) {
xn[k] <- compQ(tmpo, qua[k])
}
newrow <- data.frame(
ymin = xn[1],
yQ1 = xn[2],
yME = xn[3],
yQ3 = xn[4],
ymax = xn[5],
indiv = rownames(x@M)[i],
vars = colnames(x@M)[j],
stringsAsFactors = FALSE
)
df <- rbind(df, newrow)
}
}
df$indiv <- as.factor(df$indiv)
df$indiv <- factor(df$indiv, levels = rownames(x@M))
df$vars <- as.factor(df$vars)
df$vars <- factor(df$vars, levels = colnames(x@M))
p <- with(df, ggplot(df, aes(x = factor(0), ymin = ymin, lower = yQ1, middle = yME, upper = yQ3, ymax = ymax, fill = vars)) +
geom_boxplot(stat = "identity", size = 0.5, colour = border) +
facet_grid(vars ~ indiv, scales = "free_y") +
scale_x_discrete(breaks = NULL) +
theme(
legend.position = "none",
axis.title.x = element_blank(),
strip.text.x = element_text(size = 12, face = "bold", angle = 30),
strip.text.y = element_text(size = 12, face = "bold")
))
# print(p)
}
}
return(p)
}
# OK assigned
plot.HTS.1v <- function(x, type = "BOXPLOT", border = "black", maxno.perplot = 15) {
numofp <- length(x@data)
numoflines <- ceiling(numofp / maxno.perplot)
listofP <- list()
df1 <- summaryHTS(x)
clint <- c(1:numofp) / numofp
TT <- max(df1$Tstamp) - min(df1$Tstamp)
Tini <- min(df1$Tstamp)
for (plo in 1:numoflines) {
selected <- c(((plo - 1) * maxno.perplot + 1), min(plo * maxno.perplot, numofp))
y <- subsetHTS(x, selected[1], selected[2])
if (type == "BOXPLOT") {
df <- df1[(selected[1]:selected[2]), ]
p <- with(df, ggplot(df, aes(
x = as.factor(Tstamp), ymin = min, lower = Q1, middle = MED,
upper = Q3, ymax = Max, fill = as.factor(Tstamp)
)) +
geom_boxplot(stat = "identity") +
guides(fill = FALSE) +
theme(
legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(vjust = 0.5, size = 8)
))
# print(p)
listofP[[plo]] <- p
}
if (type == "VIOLIN") {
nr <- length(y@data)
df <- data.frame(x = numeric(), t = numeric(), grad = numeric())
df2 <- data.frame(x = numeric(), t = numeric(), xe = numeric(), te = numeric(), grad = numeric())
for (i in 1:nr) {
points <- 100
xn <- c(rep(0, points))
random_no <- c(0:points) / points
for (p in 1:points) {
xn[p] <- compQ(x@data[[i]], random_no[p])
}
newrow <- data.frame(x = xn, t = rep(y@data[[i]]@tstamp), grad = rep((y@data[[i]]@tstamp - Tini) / TT))
df <- rbind(df, newrow)
if (i == nr) {
newrow2 <- data.frame(
x = y@data[[i]]@m, t = rep(y@data[[i]]@tstamp),
xe = y@data[[i]]@m, te = rep(y@data[[i]]@tstamp),
grad = rep((y@data[[i]]@tstamp - Tini) / TT)
)
}
else {
newrow2 <- data.frame(
x = y@data[[i]]@m, t = rep(y@data[[i]]@tstamp),
xe = y@data[[i + 1]]@m, te = rep(y@data[[i + 1]]@tstamp),
grad = rep((y@data[[i]]@tstamp - Tini) / TT)
)
}
df2 <- rbind(df2, newrow2)
}
p <- with(df, ggplot(df, aes(factor(t), x, fill = grad)) +
geom_violin() +
geom_point(data = df2, aes(factor(t), x)) +
geom_line(data = df2, aes(x = factor(t), y = x, group = factor(0)), size = 0.5, alpha = 0.7, linetype = "dashed") +
scale_fill_gradient2(limits = c(0, 1), low = "red", mid = "white", high = "green", midpoint = 0.5) +
theme(
plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"),
legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 330, vjust = 0.5, size = 10)
))
listofP[[plo]] <- p
}
}
# return(listofP)
multiplot(listofP)
}
multiplot <- function(plotlist = NULL, ..., file, cols = 1, layout = NULL) {
# require(grid)
# require(ggplot2)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots <- length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots / cols)),
ncol = cols, nrow = ceiling(numPlots / cols)
)
}
if (numPlots == 1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(
layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col
))
}
}
}
#' A function for comparing observed vs predicted histograms
#' @description This function allows the representation of observed vs predicted histograms. It can
#' be used as a tool for interpreting preditive methods (for exampe, the regression of histogrma data)
#' @param PRED a \code{MatH} object with one column, the predicted data
#' @param OBS a \code{MatH} object with one column, the observed data
#' @param type a string. "HISTO" (default), if ones want to compare histograms\cr
#' "CDF", if ones want to compare cumulative distribution functions;\cr
#' "DENS" if ones want to compare approximated densities (using KDE);\cr
#' @param ncolu number of columns in which is arranged the plot, default is 2. If you have a lot of data consider
#' to choose higher values.
#' @return A plot with compared histogram-valued data.
#' @examples
#' ## do a regression
#' pars <- WH.regression.two.components(BLOOD, Yvar = 1, Xvars = c(2:3))
#' ## predict data
#' PRED <- WH.regression.two.components.predict(data = BLOOD[, 2:3], parameters = pars)
#' ## define observed data
#' \dontrun{
#' OBS <- BLOOD[, 1]
#' plotPredVsObs(PRED, OBS, "HISTO")
#' plotPredVsObs(PRED, OBS, "CDF")
#' plotPredVsObs(PRED, OBS, "DENS")
#' }
#' @export
plotPredVsObs <- function(PRED, OBS, type = "HISTO", ncolu = 2) {
# require ("vioplot")
nobj <- get.MatH.nrows(PRED)
maxpercol <- ceiling(nobj / ncolu)
ListofP <- list()
M <- data.frame()
cc <- 0
for (ind in 1:nobj) {
cc <- cc + 1
if (cc > ncolu) cc <- 1
if (type == "HISTO") {
x <- OBS@M[ind, 1][[1]]
lowers <- x@x[1:(length(x@x) - 1)]
uppers <- x@x[2:length(x@x)]
maxdens <- 10 / (max(uppers) - lowers[1])
ampl <- uppers - lowers
# ampl[which(ampl==0)]=1
dens <- (x@p[2:length(x@p)] - x@p[1:(length(x@p) - 1)]) / ampl
dens[which(dens > maxdens)] <- maxdens
df <- data.frame(
xm = lowers, xM = uppers, ym = dens * 0, yM = dens,
Type = rep("OBS", length(lowers)),
NameO = rep(rownames(OBS@M)[ind], length(lowers)),
colu = rep(cc, length(lowers))
)
x <- PRED@M[ind, 1][[1]]
lowers <- x@x[1:(length(x@x) - 1)]
uppers <- x@x[2:length(x@x)]
maxdens <- 10 / (max(uppers) - lowers[1])
ampl <- uppers - lowers
# ampl[which(ampl==0)]=1
dens <- (x@p[2:length(x@p)] - x@p[1:(length(x@p) - 1)]) / ampl
dens[which(dens > maxdens)] <- maxdens
df2 <- data.frame(
xm = lowers, xM = uppers, ym = dens * 0, yM = dens,
Type = rep("PRED", length(lowers)),
NameO = rep(rownames(OBS@M)[ind], length(lowers)),
colu = rep(cc, length(lowers))
)
M <- rbind(M, df, df2)
}
if (type == "CDF") {
x <- PRED@M[ind, 1][[1]]
xs <- x@x
ps <- x@p
df <- data.frame(
x = xs, cdf = ps,
Type = rep("PRED", length(xs)),
NameO = rep(rownames(OBS@M)[ind], length(xs)),
colu = rep(cc, length(xs))
)
x <- OBS@M[ind, 1][[1]]
xs <- x@x
ps <- x@p
df2 <- data.frame(
x = xs, cdf = ps,
Type = rep("OBS", length(xs)),
NameO = rep(rownames(OBS@M)[ind], length(xs)),
colu = rep(cc, length(xs))
)
M <- rbind(M, df, df2)
}
if (type == "DENS") {
# generate 200 random points according to the QF
rn <- 200
xn1 <- c(rep(0, rn))
xn2 <- c(rep(0, rn))
random_no <- c(0:rn) / rn
x1 <- PRED@M[ind, 1][[1]]
x2 <- OBS@M[ind, 1][[1]]
for (i in 1:rn) {
xn1[i] <- compQ(x1, random_no[i])
xn2[i] <- compQ(x2, random_no[i])
}
df1 <- data.frame(
x = xn1,
Type = rep("PRED", length(xn1)),
NameO = rep(rownames(OBS@M)[ind], length(xn1)),
colu = rep(cc, length(xn1))
)
df2 <- data.frame(
x = xn2,
Type = rep("OBS", length(xn2)),
NameO = rep(rownames(OBS@M)[ind], length(xn2)),
colu = rep(cc, length(xn2))
)
M <- rbind(M, df1, df2)
}
}
ListofP <- list()
levels(M$Type) <- c("OBS", "PRED")
for (CC in 1:ncolu) {
if (type == "HISTO") {
p <- with(M, ggplot(subset(M, colu == CC), aes(xmin = xm, xmax = xM, ymin = ym, ymax = yM, fill = Type)) +
geom_rect(alpha = 0.5) +
facet_grid(NameO ~ .) +
theme(legend.position = "bottom", legend.title = element_blank()))
}
if (type == "CDF") {
p <- with(M, ggplot(subset(M, colu == CC), aes(x = x, y = cdf)) +
geom_line(aes(linetype = Type)) +
facet_grid(NameO ~ .) +
theme(
legend.position = "bottom",
axis.title.x = element_blank(),
legend.title = element_blank()
))
}
if (type == "DENS") {
p <- with(M, ggplot(subset(M, colu == CC), aes(x = x, fill = Type)) +
geom_density(alpha = 0.6, colour = "gray") +
facet_grid(NameO ~ .) +
theme(
legend.position = "bottom",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank()
))
}
ListofP <- c(ListofP, list(p))
}
multiplot(ListofP, cols = ncolu)
}
#' A function for plotting functions of errors
#' @description This function allows the representation of the difference between observed histograms and
#' the respective predicted ones. It can
#' be used as a tool for interpreting preditive methods (for exampe, the regression of histogrma data)
#' @param PRED a \code{MatH} object with one column, the predicted data
#' @param OBS a \code{MatH} object with one column, the observed data
#' @param type a string. "HISTO_QUA" (default), if ones want to compare histograms quantile differences\cr
#' "HISTO_DEN", if ones want to show the histogram densities differences;\cr
#' "DENS_KDE" if ones want to show the differences between approximated densities (using KDE);\cr
#' @param np number of points considered for density or quantile computation (default=200).
#' @return A plot with functions of differences between observed and predicted histograms, and a
#' Root Mean Squared value computing by using the L2 Wasserstein distance.
#' @examples
#' ## do a regression
#' pars <- WH.regression.two.components(BLOOD, Yvar = 1, Xvars = c(2:3))
#' ## predict data
#' PRED <- WH.regression.two.components.predict(data = BLOOD[, 2:3], parameters = pars)
#' ## define observed data
#' OBS <- BLOOD[, 1]
#' plot_errors(PRED, OBS, "HISTO_QUA")
#' plot_errors(PRED, OBS, "HISTO_DEN")
#' plot_errors(PRED, OBS, "DENS_KDE")
#' @export
plot_errors <- function(PRED, OBS, type = "HISTO_QUA", np = 200) {
RMS_W <- 0
nobj <- get.MatH.nrows(PRED)
dists <- rep(0, nobj)
# maxpercol=ceiling(nobj/ncolu)
ListofP <- list()
M <- data.frame()
# M2=data.frame()
# M3=data.frame()
minimum <- min(min((get.MatH.stats(OBS, stat = "min"))$mat), min((get.MatH.stats(PRED, stat = "min"))$mat))
maximum <- max(max((get.MatH.stats(OBS, stat = "max"))$mat), max((get.MatH.stats(PRED, stat = "max"))$mat))
vals <- minimum + (maximum - minimum) * c(0:np) / np
range <- maximum - minimum
NM <- WH.bind.col(OBS, PRED)
colnames(NM@M) <- c("OBS", "PRED")
NM <- registerMH(NM)
for (i in 1:nobj) {
xo <- NM@M[i, 1][[1]]
xp <- NM@M[i, 2][[1]]
dists[i] <- WassSqDistH(xo, xp)
RMS_W <- RMS_W + 1 / nobj * dists[i]
if (type == "HISTO_QUA") {
# tmp=register(xo,xp)
x <- rep(0, (np + 1))
p <- rep(0, (np + 1))
for (j in 0:np) {
p[j + 1] <- j / np
x[j + 1] <- compQ(xo, p[j + 1]) - compQ(xp, p[j + 1])
# x=xo@x-xp@x
# p=xo@p
}
df <- data.frame(quantile = x, p = p, NameO = rep(rownames(OBS@M)[i], length(x)))
M <- rbind(M, df)
}
if (type == "HISTO_DEN") {
fo <- rep(0, np)
fp <- rep(0, np)
mp <- 0.5 * (vals[2:np] + vals[1:(np - 1)])
ran <- (vals[2:np] - vals[1:(np - 1)])
for (v in 1:np) {
fo[v] <- compP(xo, vals[v])
fp[v] <- compP(xp, vals[v])
}
do <- (fo[2:np] - fo[1:(np - 1)]) / ran
dp <- (fp[2:np] - fp[1:(np - 1)]) / ran
diffden <- do - dp
df2 <- data.frame(x = mp, density = diffden, NameO = rep(rownames(OBS@M)[i], length(mp)))
M <- rbind(M, df2)
}
if (type == "DENS_KDE") {
# smooth
po <- rep(0, np)
pp <- rep(0, np)
for (v in 1:np) {
po[v] <- compQ(xo, (v / np))
pp[v] <- compQ(xp, (v / np))
}
do <- density(po, from = minimum, to = maximum)
dp <- density(pp, from = minimum, to = maximum)
diffdkde <- do$y - dp$y
df3 <- data.frame(x = do$x, density = diffdkde, NameO = rep(rownames(OBS@M)[i], length(do$x)))
M <- rbind(M, df3)
}
}
if (type == "HISTO_QUA") {
p <- with(M, ggplot(data = M, aes(x = p, y = quantile, group = NameO, colour = NameO)) +
geom_line() +
geom_point() +
ggtitle("Plot of quantile differences"))
}
if (type == "HISTO_DEN") {
p <- with(M, ggplot(data = M, aes(x = x, y = density, group = NameO, colour = NameO)) +
geom_line() +
ggtitle("Plot of density differences"))
}
if (type == "DENS_KDE") {
p <- with(M, ggplot(data = M, aes(x = x, y = density, group = NameO, colour = NameO)) +
geom_line() +
ggtitle("Plot of smoothed (KDE) \n density differences"))
}
print(p)
RMS_W <- sqrt(RMS_W)
print(RMS_W)
# print(RMS_W/range)
# print(sqrt(dists))
# summary(sqrt(dists))
}
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.