R/rlm_stan.R

# Save this file as `R/rlm_stan.R`

#' Bayesian robust linear regression with Stan
#'
#' @export
#' @param x Numeric matrix of input values. This matrix should not contain
#' a column of 1's.
#' @param y Numeric vector of output values.
#' @param nu Character, either "fixed" for a fixed value of nu (I'm going to make it 4,
#' since this is a fake package), "continuous", in which case the prior distribution
#' for nu is gamma(shape = 2, rate = 0.5), for an expected value of 4, or "discrete", in
#' which case the prior distribution for nu is a discrete uniform(1, 10). In a real package,
#' I'd give more input options (fixed parameter value, shape, rate, upper bound, ...),
#' but I don't care here
#' @param ... Arguments passed to `rstan::sampling` (e.g. iter, chains).
#' @return An object of class `stanfit` returned by `rstan::sampling`
#' @examples
#' \dontrun{
#' x <- runif(30, 0, 10)
#' y <- 1 + 2*x + rt(30, 4)
#' rlmFitFixed <- rlm_stan(x, y, "fixed")
#' rlmFitCont <- rlm_stan(x, y, "continuous")
#' rlmFitDisc <- rlm_stan(x, y, "discrete")
#' }
rlm_stan <- function(x, y, nu = c("fixed", "continuous", "discrete"), ...) {
  if(is.vector(x)) x <- matrix(x)
  nu <- match.arg(nu)
  standata <- list(x = x, y = y, K = ncol(x), N = length(y))
  standata$nu <- switch(nu,
                        fixed = 4,
                        continuous = 1,
                        discrete = 1)
  if(nu == "fixed") {
    out <- rstan::sampling(stanmodels$rlmNuFixed, data = standata, ...)
  }else if(nu == "continuous"){
    out <- rstan::sampling(stanmodels$rlmNuEstCont, data = standata, ...)
  }else{
    out <- rstan::sampling(stanmodels$rlmNuEstDisc, data = standata, ...)
  }
  return(out)
}
cbdavis33/rstanPullTest documentation built on June 4, 2019, 11:48 a.m.