## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4)

## ----setup--------------------------------------------------------------------
library(clustEff)

## ----waves--------------------------------------------------------------------
set.seed(1234)
tt <- seq(0, 1, length.out = 200)
sdt <- 0.3 * sqrt(pmax(4 * pmax(tt - 0.2, 0) - 8 * pmax(tt - 0.5, 0) + 4 * pmax(tt - 0.8, 0), 0))
W <- cbind(sapply(1:6, function(i) sin(3 * pi * tt) + rnorm(200, 0, sdt)),
           sapply(1:6, function(i) cos(3 * pi * tt) + rnorm(200, 0, sdt)),
           sapply(1:4, function(i) sin(3 * pi * tt) * cos(pi * tt) + rnorm(200, 0, sdt)))

## ----waves-fit----------------------------------------------------------------
fit <- clustEff(W, p = tt, k = c(2, 6), smooth = TRUE)
fit
summary(fit)

## ----waves-plot, fig.height = 7-----------------------------------------------
oldpar <- par(mfrow = c(2, 2))
plot(fit, "dendrogram")
plot(fit, "clusters")
plot(fit, "silhouette")
plot(fit, "criterion")
par(oldpar)

## ----qrcm, eval = requireNamespace("qrcm", quietly = TRUE)--------------------
set.seed(1)
n <- 500
x <- runif(n, 0, 5)
slope <- list(function(u) 0.5 + 0.5 * u + u^2,
              function(u) -0.5 + 0.5 * u + u^2,
              function(u) 0.3 - 0.5 * u - u^2 + 2 * u^3)
fits <- lapply(1:9, function(r) {
  u <- runif(n)
  y <- 1 + qnorm(u) + slope[[(r - 1) %/% 3 + 1]](u) * x
  qrcm::iqr(y ~ x, formula.p = ~ I(p) + I(p^2) + I(p^3))
})
names(fits) <- paste0("y", 1:9)
ec <- extract_curves(fits, which = "x")
ec
fit2 <- clustEff(ec, k = 3)
summary(fit2)
oldpar <- par(mfrow = c(1, 2))
plot(fit2, "dendrogram")
plot(fit2, "clusters")
par(oldpar)

## ----sim----------------------------------------------------------------------
p <- seq(0.01, 0.99, 0.01)
th <- list(c(.5, .5, 1, 2), c(-3, .5, 1, .5), c(.3, -.5, -1, 2), c(2.5, .5, 1, 2))
crits <- c("height", "pi.out", "dist.rel", "silhouette")
set.seed(3)
res <- t(sapply(1:12, function(r) {
  G <- 2 + r %% 3
  B <- do.call(cbind, lapply(1:G, function(g) sapply(1:8, function(i) {
    t <- th[[g]] + rnorm(4, 0, 0.3); t[1] + t[2] * p + t[3] * p^2 + t[4] * p^3 })))
  sapply(crits, function(cr) clustEff(B, B - 0.6, B + 0.6, p = p, k = c(2, 7), criterion = cr)$k == G)
}))
colMeans(res)

## ----fpcac--------------------------------------------------------------------
W2 <- cbind(W[, 1:12], sapply(1:2, function(i) rnorm(200, 0, 0.3)))   # two outliers
set.seed(1)
o <- opt.fpcac(W2, t = tt, k.max = 5, nbasis = 12, alpha = 0.15)
o
f <- fpcac(W2, t = tt, K = o$K.opt, nbasis = 12, alpha = 0.15)
f
plot(f)