Produce violin-plot using base graphics.

violplot.R

# The syntax is very much as in `boxplot` (a large portion of the code is taken 
# from graphics::boxplot.default and graphics::bxp).
# The argument `density.args` should be a named list and is used to pass arguments to `density`.

violplot <-
function (x, ...) 
UseMethod("violplot")

registerS3method("violplot", "formula", 
function (formula, data = NULL, ..., subset, na.action = NULL, 
    xlab = mklab(horizontal), ylab = mklab(!horizontal), 
    add = FALSE, ann = !add, horizontal = FALSE,
    drop = FALSE, sep = ".", lex.order = FALSE
    ) {
    if (missing(formula) || (length(formula) != 3L)) 
        stop("'formula' missing or incorrect")
    if (missing(xlab) || missing(ylab)) 
        mklab <- function(y_var) if (y_var) names(mf)[response] else
            paste(names(mf)[-response], collapse = " : ")
    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval(m$data, parent.frame()))) 
        m$data <- as.data.frame(data)
    m$... <- m$drop <- m$sep <- m$lex.order <- NULL
    m$xlab <- m$ylab <- m$add <- m$ann <- m$horizontal <- NULL
    m$na.action <- na.action
    m[[1L]] <- quote(stats::model.frame.default)
    mf <- eval(m, parent.frame())
    response <- attr(attr(mf, "terms"), "response")
    violplot(split(mf[[response]], mf[-response], drop = drop, 
        sep = sep, lex.order = lex.order), xlab = xlab, ylab = ylab, 
        add = add, ann = ann, horizontal = horizontal, ...)
})

registerS3method("violplot", "default", 
function (x, ..., width = NULL, varwidth = FALSE, 
    names,
    density.args = list(),
    border = par("fg"), col = "lightgray",
    pars = list(),
    ann = !add, axes = !add, horizontal = FALSE, add = FALSE, at = NULL) {

    .names <- function(x) attr(x, "names")
    `.names<-` <- function(x, value) `attr<-`(x, "names", value)

    groups <- if (is.list(x)) x else {
        args <- list(x, ...)
        if(!is.null(.names(args))) {
            args[!.names(args) != ""]
        } else args
    }

    pars <- as.list(pars)
    if (...length()) {
        argnames <- .names(args <- list(...))
        if (!is.null(argnames)) {
            args <- args[namedargs <- argnames != ""]
            argnames <- argnames[namedargs]
            if (anyDuplicated(argnames)) {
                dup <- duplicated(argnames)
                warning(sprintf(ngettext(sum(dup), "Duplicated argument %s is disregarded", 
                    "Duplicated arguments %s are disregarded"), sub("^list\\((.*)\\)", 
                    "\\1", deparse(args[dup]))), domain = NA)
                argnames <- .names(args <- args[!dup])
            }
            pars[argnames] <- args
        }
    }

    if (0L == (n <- length(groups))) stop("invalid first argument")
    if (length(class(groups))) groups <- unclass(groups)
    if (!missing(names))
        .names(groups) <- names else {
        if (is.null(.names(groups))) 
            .names(groups) <- 1L:n
        names <- .names(groups)
    }
    
    if (is.null(at)) at <- seq_len(n)  else if (length(at) != n) 
        stop(gettextf("'at' must have a length of 'n', i.e. %d", n), domain = NA)
        
    range <- range(unlist(lapply(groups, range, na.rm = TRUE)))
  
    pcycle <- function(p, d1, d2 = NULL)
        rep_len(if (length(p)) p else if (length(d1)) d1 else d2, 
            length.out = n)
    p <- function(sym) pars[[sym, exact = TRUE]]

    boxlty <- pcycle(pars$boxlty, p("lty"), par("lty"))
    boxlwd <- pcycle(pars$boxlwd, p("lwd"), par("lwd"))
    boxcol <- pcycle(pars$boxcol, border, par("fg"))
    boxfill <- pcycle(pars$boxfill, col, par("bg"))
    boxangle <- pcycle(pars$boxangle, p("angle"), 45)
    boxdensity <- pcycle(pars$boxdensity, p("density"), NA)

    argdflt <- function(a, b) if(length(a)) a else b
   
    density.args <- as.list(density.args)
    density.args <- density.args[.names(density.args) != ""]
    density.args$from <- argdflt(density.args$from, range[1L])
    density.args$to <- argdflt(density.args$to, range[2L])
    density.args$n <- dn <-
        argdflt(density.args$n[1L], eval(formals(stats::density.default)$n))
    density.args$na.rm <- TRUE

    dens <- do.call("lapply", c(list(groups, "density"), density.args))
    px <- dens[[1L]]$x
    py <- sapply(dens, "[[", "y")
    py <- py / max(py) 

    #xlog <- (par("ylog") && horizontal) || (par("xlog") && !horizontal)
    log <- ""
    
    ng <- vapply(groups, length, 0L)

    boxwex <- pcycle(pars$boxwex, 0.45 * {
        if (n <= 1) 1 else min(diff(sort(at)))
    })    
    width <- if (!is.null(width)) {
        if (length(width) != n | anyNA(width) | any(width <= 0)) 
            stop("invalid boxplot widths")
        boxwex * width / max(width)
    } else if (varwidth) {
        boxwex * sqrt(ng / max(ng))
    } else boxwex
        
    dev.hold()
    on.exit(dev.flush())    
    if (!add) {
        ylim <- if (is.null(pars$ylim)) range else pars$ylim
        xlim <- if (is.null(pars$xlim)) 
            base::range(at, finite = TRUE) + c(-0.5, 0.5) else pars$xlim

        plot.new()
        if (horizontal)
            plot.window(ylim = xlim, xlim = ylim, log = log, xaxs = pars$yaxs)
        else
            plot.window(xlim = xlim, ylim = ylim, log = log, yaxs = pars$yaxs)
    }

    xy <- matrix(0, ncol = 2L, nrow = dn * 2L)
    idx <- if(horizontal) c(1L, 2L) else c(2L, 1L)
    for(i in seq.int(n)) {
        xy[, idx[1L]] <- c(px, rev(px))
        xy[, idx[2L]] <- at[i] + (width[i] * c(py[, i], rev(-py[, i])))
            polygon(xy,
            col = boxfill[i], border = boxcol[i],
            lwd = boxlwd[i], lty = boxlty[i],
            angle = boxangle[i], density = boxdensity[i])
    }
    
    if (axes) {
        ax.pars <- pars[.names(pars) %in% c("xaxt", "yaxt", "xaxp", 
            "yaxp", "las", "cex.axis", "col.axis", "format")]
        if (n > 1) 
            do.call("axis", c(list(side = 1L + horizontal, at = at,
                labels = names), ax.pars))
        do.call("Axis", c(list(x = py, side = 2L - horizontal), ax.pars))
        box()
    }
    if (ann) 
        do.call(title, pars[.names(pars) %in% c("main", "cex.main", 
            "col.main", "sub", "cex.sub", "col.sub", "xlab", 
            "ylab", "cex.lab", "col.lab")])

    invisible(at)
})

Example:

z <- data.frame(y = NA, f = gl(5, 20))
z$y <- rnorm(nrow(z), seq.int(5)[z$f], 2)

par(mfrow = c(1, 3))
# formula interface:
violplot(y ~ f, z, col = 1:3, lwd = 1,
         density = c(35, NA), angle = c(45, 90, 180),
         horizontal = TRUE)
violplot(y ~ f, z, col = "#8e9db0", border = "gray10",
         density.args = list(bw = .25),
         horizontal = FALSE)

# pass group vectors directly as unnamed arguments:
violplot(rnorm(100, 1), rpois(100, 2), runif(1000, 0, 10),
         names = c("normal", "Poisson", "uniform"))