## silence output and warnings

SW <- function(expr) suppressMessages(suppressWarnings(expr))

## dataset

set.seed(1)

N <- 50

P <- 10

U <- 3

x <- matrix(rnorm(N * P), nrow=N, ncol=P)

b <- runif(P) - 0.5

y.gauss <- rnorm(N, mean=x %*% b, sd=runif(1, 1, 2))

y.binom <- factor(rbinom(N, 1, binomial()$linkinv(x %*% b)))

df <- data.frame(x, y.gauss=y.gauss, y.binom=y.binom)

df[, 1] <- factor(letters[rbinom(N, 2, 0.5) + 1])

df$X1b_X3 <- df$X3 * (df$X1 == "b")

df$X1c_X3 <- df$X3 * (df$X1 == "c")

df$X3_X2 <- df$X3 * df$X2

## model options

unp <- paste0("X", 1:U)

pen <- setdiff(paste0("X", 1:P), unp)

mod.gauss <- reformulate(unp, "y.gauss")

mod.binom <- reformulate(unp, "y.binom")

folds <- c(rep(1, N / 2), rep(2, N / 2))

iters <- 500

chains <- 2

## numerical tolerance

tol <- 0.000001

## wrapper to set commonly used options

hs <- function(model, family, ...)

    hsstan(df, model, pen, iter=iters, chains=chains, family=family,

           refresh=0, ...)

message("Running hsstan models...")

SW({

    hs.base <- hsstan(df, mod.gauss, iter=iters, chains=chains, refresh=0)

    hs.gauss <- hs(mod.gauss, gaussian)

    hs.binom <- hs(mod.binom, binomial)

    hs.inter <- hs(y.gauss ~ X1 * X3 + X2 * X3, gaussian)

})

message("Running cross-validated hsstan models...")

SW({

    cv.gauss <- kfold(hs.gauss, folds=folds, chains=2)

    cv.binom <- kfold(hs.binom, folds=folds)

    cv.nofit <- cv.gauss

    cv.nofit$fits <- cv.nofit$data <- NULL

})