diff --git a/NAMESPACE b/NAMESPACE
index 3b10b79a..180b3e22 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,6 +13,8 @@ S3method(dim,kfold)
 S3method(dim,loo)
 S3method(dim,psis_loo)
 S3method(dim,waic)
+S3method(elpd,array)
+S3method(elpd,matrix)
 S3method(importance_sampling,array)
 S3method(importance_sampling,default)
 S3method(importance_sampling,matrix)
@@ -73,6 +75,7 @@ export(.ndraws)
 export(.thin_draws)
 export(E_loo)
 export(compare)
+export(elpd)
 export(example_loglik_array)
 export(example_loglik_matrix)
 export(extract_log_lik)
diff --git a/R/elpd.R b/R/elpd.R
new file mode 100644
index 00000000..d552acda
--- /dev/null
+++ b/R/elpd.R
@@ -0,0 +1,70 @@
+#' Generic (expected) log-predictive density
+#'
+#' The `elpd()` methods for arrays, and matrices can compute the expected log pointwise predictive density for a new dataset or the log pointwise predictive density of the observed data (an overestimate of the elpd).
+#'
+#' @export
+#' @inheritParams loo
+#'
+#' @details The `elpd()` function is an S3 generic and methods are provided for
+#'   3-D pointwise log-likelihood arrays, pointwise log-likelihood matrices.
+#'
+#'
+#' @examples
+#' ### Array methods for calculating the lpd of the observed data (using example objects included with loo package)
+#' LLarr <- example_loglik_array()
+#' elpd(LLarr)
+#'
+#'
+#'
+elpd <- function(x, ...) {
+  UseMethod("elpd")
+}
+
+#' @export
+#' @templateVar fn elpd
+#' @template array
+#'
+elpd.array <- function(x, ...) {
+
+  ll <- llarray_to_matrix(x)
+  elpd.matrix(ll)
+}
+
+#' @export
+#' @templateVar fn elpd
+#' @template matrix
+#'
+elpd.matrix <-
+  function(x, ...) {
+  pointwise <- pointwise_elpd_calcs(x)
+  elpd_object(pointwise, dim(x))
+}
+
+
+pointwise_elpd_calcs <- function(ll){
+  elpd <- colLogSumExps(ll) - log(nrow(ll))
+  ic <- -2 * elpd
+  cbind(elpd, ic)
+}
+
+elpd_object <- function(pointwise, dims) {
+  if (!is.matrix(pointwise)) stop("Internal error ('pointwise' must be a matrix)")
+
+  cols_to_summarize <- colnames(pointwise)
+  estimates <- table_of_estimates(pointwise[, cols_to_summarize, drop=FALSE])
+
+  out <- nlist(estimates, pointwise)
+
+  structure(
+    out,
+    dims = dims,
+    class = c("elpd_generic", "loo")
+  )
+}
+print_dims.elpd_generic <- function(x, ...) {
+    cat(
+    "Computed from",
+    paste(dim(x), collapse = " by "),
+    "log-likelihood matrix using the generic elpd function\n"
+  )
+}
diff --git a/man/elpd.Rd b/man/elpd.Rd
new file mode 100644
index 00000000..be4b59b6
--- /dev/null
+++ b/man/elpd.Rd
@@ -0,0 +1,50 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/elpd.R
+\name{elpd}
+\alias{elpd}
+\alias{elpd.array}
+\alias{elpd.matrix}
+\title{Generic (expected) log-predictive density}
+\usage{
+elpd(x, ...)
+
+\method{elpd}{array}(x, ...)
+
+\method{elpd}{matrix}(x, ...)
+}
+\arguments{
+\item{x}{A log-likelihood array, matrix, or function. The \strong{Methods (by class)}
+section, below, has detailed descriptions of how to specify the inputs for
+each method.}
+
+\item{...}{For the \code{loo.function()} method and the \code{loo_i()}
+function, these are the data, posterior draws, and other arguments to pass
+to the log-likelihood function. See the \strong{Methods (by class)} section
+below for details on how to specify these arguments.}
+}
+\description{
+The \code{elpd()} methods for arrays, and matrices can compute the expected log pointwise predictive density for a new dataset or the log pointwise predictive density of the observed data (an overestimate of the elpd).
+}
+\details{
+The \code{elpd()} function is an S3 generic and methods are provided for
+3-D pointwise log-likelihood arrays, pointwise log-likelihood matrices.
+}
+\section{Methods (by class)}{
+\itemize{
+\item \code{array}: An \eqn{I} by \eqn{C} by \eqn{N} array, where \eqn{I}
+is the number of MCMC iterations per chain, \eqn{C} is the number of
+chains, and \eqn{N} is the number of data points.
+
+\item \code{matrix}: An \eqn{S} by \eqn{N} matrix, where \eqn{S} is the size
+of the posterior sample (with all chains merged) and \eqn{N} is the number
+of data points.
+}}
+
+\examples{
+### Array methods for calculating the lpd of the observed data (using example objects included with loo package)
+LLarr <- example_loglik_array()
+elpd(LLarr)
+
+
+
+}
diff --git a/tests/testthat/reference-results/elpd.rds b/tests/testthat/reference-results/elpd.rds
new file mode 100644
index 00000000..f8bcfb36
Binary files /dev/null and b/tests/testthat/reference-results/elpd.rds differ
diff --git a/tests/testthat/test_loo_and_waic.R b/tests/testthat/test_loo_and_waic.R
index 519f00eb..982b1334 100644
--- a/tests/testthat/test_loo_and_waic.R
+++ b/tests/testthat/test_loo_and_waic.R
@@ -2,7 +2,7 @@ library(loo)
 options(mc.cores = 1)
 set.seed(123)
 
-context("loo and waic")
+context("loo, waic and elpd")
 
 LLarr <- example_loglik_array()
 LLmat <- example_loglik_matrix()
@@ -13,6 +13,7 @@ r_eff_mat <- relative_eff(exp(LLmat), chain_id = chain_id)
 
 loo1 <- suppressWarnings(loo(LLarr, r_eff = r_eff_arr))
 waic1 <- suppressWarnings(waic(LLarr))
+elpd1 <- suppressWarnings(elpd(LLarr))
 
 test_that("using loo.cores is deprecated", {
   options(mc.cores = NULL)
@@ -22,9 +23,10 @@ test_that("using loo.cores is deprecated", {
   options(mc.cores = 1)
 })
 
-test_that("loo and waic results haven't changed", {
+test_that("loo, waic and elpd results haven't changed", {
   expect_equal_to_reference(loo1, "reference-results/loo.rds")
   expect_equal_to_reference(waic1, "reference-results/waic.rds")
+  expect_equal_to_reference(elpd1, "reference-results/elpd.rds")
 })
 
 test_that("loo with cores=1 and cores=2 gives same results", {
@@ -86,6 +88,24 @@ test_that("loo returns object with correct structure", {
   expect_equal(dim(loo1), dim(LLmat))
 })
 
+
+test_that("elpd returns object with correct structure", {
+  expect_true(is.loo(elpd1))
+  expect_named(
+    elpd1,
+    c(
+      "estimates",
+      "pointwise"
+    )
+  )
+  est_names <- dimnames(elpd1$estimates)
+  expect_equal(est_names[[1]], c("elpd", "ic"))
+  expect_equal(est_names[[2]], c("Estimate", "SE"))
+  expect_equal(colnames(elpd1$pointwise), est_names[[1]])
+  expect_equal(dim(elpd1), dim(LLmat))
+})
+
+
 test_that("two pareto k values are equal", {
   expect_identical(loo1$pointwise[,"influence_pareto_k"], loo1$diagnostics$pareto_k)
 })
@@ -111,9 +131,15 @@ test_that("waic.array and waic.matrix give same result", {
   expect_identical(waic1, waic2)
 })
 
-test_that("loo and waic error with vector input", {
+test_that("elpd.array and elpd.matrix give same result", {
+  elpd2 <- suppressWarnings(elpd(LLmat))
+  expect_identical(elpd1, elpd2)
+})
+
+test_that("loo, waic, and elpd error with vector input", {
   expect_error(loo(LLvec), regexp = "no applicable method")
   expect_error(waic(LLvec), regexp = "no applicable method")
+  expect_error(elpd(LLvec), regexp = "no applicable method")
 })
 
 
diff --git a/vignettes/elpd.Rmd b/vignettes/elpd.Rmd
new file mode 100644
index 00000000..d4f943c3
--- /dev/null
+++ b/vignettes/elpd.Rmd
@@ -0,0 +1,215 @@
+---
+title: "Holdout validation and K-fold cross-validation of Stan programs with the loo package" 
+author: "Bruno Nicenboim"
+date: "`r Sys.Date()`"
+output:
+  html_vignette:
+    toc: yes
+params:
+  EVAL: !r identical(Sys.getenv("NOT_CRAN"), "true")
+---
+<!--
+%\VignetteEngine{knitr::rmarkdown}
+%\VignetteIndexEntry{Holdout validation and K-fold cross-validation of Stan programs with the loo package}
+-->
+```{r, child="children/SETTINGS-knitr.txt"}
+```
+
+```{r, child="children/SEE-ONLINE.txt", eval = if (isTRUE(exists("params"))) !params$EVAL else TRUE}
+```
+
+# Introduction
+
+This vignette demonstrates how to do holdout validation and K-fold cross-validation with __loo__ for a Stan program.
+
+
+# Example: Eradication of Roaches using holdout validation approach
+
+This vignette uses the same example as in the vignettes
+[_Using the loo package (version >= 2.0.0)_](http://mc-stan.org/loo/articles/loo2-example.html) and [_Avoiding model refits in leave-one-out cross-validation with moment matching_](https://mc-stan.org/loo/articles/loo2-moment-matching.html).
+
+
+
+## Coding the Stan model
+
+Here is the Stan code for fitting a Poisson regression model:
+
+```{r stancode}
+stancode <- "
+data {
+  int<lower=1> K;
+  int<lower=1> N;
+  matrix[N,K] x;
+  int y[N];
+  vector[N] offset;
+
+  real beta_prior_scale;
+  real alpha_prior_scale;
+}
+parameters {
+  vector[K] beta;
+  real intercept;
+}
+model {
+  y ~ poisson(exp(x * beta + intercept + offset));
+  beta ~ normal(0,beta_prior_scale);
+  intercept ~ normal(0,alpha_prior_scale);
+}
+generated quantities {
+  vector[N] log_lik;
+  for (n in 1:N)
+    log_lik[n] = poisson_lpmf(y[n] | exp(x[n] * beta + intercept + offset[n]));
+}
+"
+```
+
+Following the usual approach recommended in 
+[_Writing Stan programs for use with the loo package_](http://mc-stan.org/loo/articles/loo2-with-rstan.html), 
+we compute the log-likelihood for each observation in the 
+`generated quantities` block of the Stan program.
+
+
+## Setup
+
+In addition to __loo__, we load the __rstan__ package for fitting the model.
+We will also need the __rstanarm__ package for the data.
+
+```{r setup, message=FALSE}
+library("rstan")
+library("loo")
+seed <- 9547
+set.seed(seed)
+```
+
+# Holdout validation
+
+For this approach, the model is first fit to the "train" data and then is evaluated on the held-out "test" data.
+
+## Splitting the data between train and test 
+
+The data is divided between train (80% of the data) and test (20%):
+
+```{r modelfit-holdout, message=FALSE}
+# Prepare data
+data(roaches, package = "rstanarm")
+roaches$roach1 <- sqrt(roaches$roach1)
+roaches$offset <- log(roaches[,"exposure2"])
+# 20% of the data goes to the test set:
+roaches$test <- 0
+roaches$test[sample(.2 * seq_len(nrow(roaches)))] <- 1
+# data to "train" the model
+data_train <- list(y = roaches$y[roaches$test == 0],
+                   x = as.matrix(roaches[roaches$test == 0,
+                                         c("roach1", "treatment", "senior")]),
+                   N = nrow(roaches[roaches$test == 0,]),
+                   K = 3,
+                   offset = roaches$offset[roaches$test == 0],
+                   beta_prior_scale = 2.5,
+                   alpha_prior_scale = 5.0
+                   )
+# data to "test" the model
+data_test <- list(y = roaches$y[roaches$test == 1],
+                   x = as.matrix(roaches[roaches$test == 1,
+                                         c("roach1", "treatment", "senior")]),
+                   N = nrow(roaches[roaches$test == 1,]),
+                   K = 3,
+                   offset = roaches$offset[roaches$test == 1],
+                   beta_prior_scale = 2.5,
+                   alpha_prior_scale = 5.0
+                   )
+```
+
+
+## Fitting the model with RStan
+
+Next we fit the model to the "test" data in Stan using the __rstan__ package:
+
+```{r fit-train}
+# Compile
+stanmodel <- stan_model(model_code = stancode)
+# Fit model
+fit <- sampling(stanmodel, data = data_train, seed = seed, refresh = 0)
+```
+
+We recompute the generated quantities using the posterior draws conditional on the training data, but we now pass in the held-out data to get the log predictive densities for the test data. Because we are using independent data, the log predictive density coincides with the log likelihood of the test data.
+
+```{r gen-test}
+gen_test <- gqs(stanmodel, draws = as.matrix(fit), data= data_test)
+log_pd <- extract_log_lik(gen_test)
+```
+
+## Computing holdout elpd:
+
+Now we evaluate the predictive performance of the model on the test data using `elpd()`.
+
+```{r elpd-holdout}
+(elpd_holdout <- elpd(log_pd))
+```
+
+When one wants to compare different models, the function `loo_compare()` can be used to assess the difference in performance.
+
+# K-fold cross validation
+
+For this approach the data is divided into folds, and each time one fold is tested while the rest of the data is used to fit the model (see Vehtari et al., 2017).
+
+## Splitting the data in folds
+
+We use the data that is already pre-processed and we divide it in 10 random folds using `kfold_split_random`
+
+```{r prepare-folds, message=FALSE}
+# Prepare data
+roaches$fold <- kfold_split_random(K = 10, N = nrow(roaches))
+```
+
+
+## Fitting and extracting the log pointwise predictive densities for each fold
+
+We now loop over the 10 folds. In each fold we do the following. First, we fit the model to all the observations except the ones belonging to the left-out fold. Second, we compute the log pointwise predictive densities for the left-out fold. Last, we store the predictive density for the observations of the left-out fold in a matrix. The output of this loop is a matrix of the log pointwise predictive densities of all the observations.
+
+```{r}
+# Prepare a matrix with the number of post-warmup iterations by number of observations:
+log_pd_kfold <- matrix(nrow = 4000, ncol = nrow(roaches))
+# Loop over the folds
+for(k in 1:10){
+  data_train <- list(y = roaches$y[roaches$fold != k],
+                   x = as.matrix(roaches[roaches$fold != k,
+                                         c("roach1", "treatment", "senior")]),
+                   N = nrow(roaches[roaches$fold != k,]),
+                   K = 3,
+                   offset = roaches$offset[roaches$fold != k],
+                   beta_prior_scale = 2.5,
+                   alpha_prior_scale = 5.0
+                   )
+  data_test <- list(y = roaches$y[roaches$fold == k],
+                   x = as.matrix(roaches[roaches$fold == k,
+                                         c("roach1", "treatment", "senior")]),
+                   N = nrow(roaches[roaches$fold == k,]),
+                   K = 3,
+                   offset = roaches$offset[roaches$fold == k],
+                   beta_prior_scale = 2.5,
+                   alpha_prior_scale = 5.0
+                   )
+  fit <- sampling(stanmodel, data = data_train, seed = seed, refresh = 0)
+  gen_test <- gqs(stanmodel, draws = as.matrix(fit), data= data_test)
+  log_pd[, roaches$fold == k] <- extract_log_lik(gen_test)
+}
+```
+
+## Computing K-fold elpd:
+
+Now we evaluate the predictive performance of the model on the 10 folds using `elpd()`.
+
+```{r elpd-kfold}
+(elpd_kfold <- elpd(log_pd_kfold))
+```
+
+If one wants to compare several models (with `loo_compare`), one should use the same folds for all the different models. 
+
+# References
+
+Gelman, A., and Hill, J. (2007).  *Data Analysis Using Regression and Multilevel Hierarchical Models.*  Cambridge University Press.
+
+Stan Development Team (2020) _RStan: the R interface to Stan, Version 2.21.1_   https://mc-stan.org
+
+Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. _Statistics and Computing_. 27(5), 1413--1432. \doi:10.1007/s11222-016-9696-4. Links: [published](http://link.springer.com/article/10.1007\%2Fs11222-016-9696-4) | [arXiv preprint](http://arxiv.org/abs/1507.04544).
+