Summary: Typically, in order to characterize the homogenized effective macroscopic response of new materials possessing random heterogeneous microstructure, a relation between averages \(\left\langle \sigma \right\rangle_\Omega = \mathbb E^* : \left\langle \epsilon \right\rangle_\Omega\) is sought, where \(\left\langle \cdot \right\rangle_\Omega \overset{\text{def}}= \frac{1}{| \Omega |} \int_\Omega \cdot \,d\Omega\), and where \(\sigma\) and \(\epsilon\) are the stress and strain tensor fields within a statistically representative volume element (SRVE) of volume \(| \Omega |\). The quantity \(\mathbb E^{*}\) is known as the effective property, and is the elasticity tensor used in usual macroscale analyses. In order to generate homogenized responses computationally, a series of detailed boundary value representations resolving the heterogeneous microstructure, posed over the SRVE’s domain, must be solved. This requires an enormous numerical effort that can overwhelm most computational facilities. A natural way of generating an approximation to the SRVE’s response is by first computing the response of smaller (subrepresentative) samples, each with a different random realization of the microstructural type under investigation, and then to ensemble average the results afterwards. Compared to a direct simulation of an SRVE, testing many small samples is a computationally inexpensive process since the number of floating point operations is greatly reduced, as well as the fact that the samples’ responses can be computed trivially in parallel. However, there is an inherent error in this process. Clearly the population’s ensemble average is not the SRVE’s response. However, as shown in this work, the moments on the distribution of the population can be used to generate rigorous upper and lower error bounds on the quality of the ensemble-generated response. Two-sided bounds are given on the SRVE response in terms of the ensemble average, its standard deviation and its skewness.