diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index bdc816a341..ec340d9d51 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -6,8 +6,9 @@
 
 - Track the model log-likelihood as a sampler stat for NUTS and HMC samplers
   (accessible as `trace.get_sampler_stats('model_logp')`) (#3134)
-- Add Incomplete Beta function `incomplete_beta(a, b, value)` 
+- Add Incomplete Beta function `incomplete_beta(a, b, value)`
 - Add log CDF functions to continuous distributions: `Beta`, `Cauchy`, `ExGaussian`, `Exponential`, `Flat`, `Gumbel`, `HalfCauchy`, `HalfFlat`, `HalfNormal`, `Laplace`, `Logistic`, `Lognormal`, `Normal`, `Pareto`, `StudentT`, `Triangular`, `Uniform`, `Wald`, `Weibull`.
+- Behavior of `sample_posterior_predictive` is now to produce posterior predictive samples, in order, from all values of the `trace`. Previously, by default it would produce 1 chain worth of samples, using a random selection from the `trace` (#3212)
 
 ### Maintenance
 
diff --git a/pymc3/sampling.py b/pymc3/sampling.py
index 5806a42ee5..ac54def005 100644
--- a/pymc3/sampling.py
+++ b/pymc3/sampling.py
@@ -1098,7 +1098,7 @@ def sample_posterior_predictive(trace, samples=None, model=None, vars=None, size
         nchain = 1
 
     if samples is None:
-        samples = len(trace)
+        samples = sum(len(v) for v in trace._straces.values())
 
     model = modelcontext(model)
 
@@ -1108,7 +1108,7 @@ def sample_posterior_predictive(trace, samples=None, model=None, vars=None, size
     if random_seed is not None:
         np.random.seed(random_seed)
 
-    indices = np.random.randint(0, nchain * len_trace, samples)
+    indices = np.arange(samples)
 
     if progressbar:
         indices = tqdm(indices, total=samples)
@@ -1126,9 +1126,9 @@ def sample_posterior_predictive(trace, samples=None, model=None, vars=None, size
         for slc, idx in enumerate(indices):
             if nchain > 1:
                 chain_idx, point_idx = np.divmod(idx, len_trace)
-                param = trace._straces[chain_idx].point(point_idx)
+                param = trace._straces[chain_idx % nchain].point(point_idx)
             else:
-                param = trace[idx]
+                param = trace[idx % len_trace]
 
             values = draw_values(vars, point=param, size=size)
             for k, v in zip(vars, values):