A maximum entropy (MaxEnt) based probabilistic approach is developed to model mechanical systems characterized by symmetric positive-definite matrices bounded from below and above. These matrices are typically encountered in the constitutive modeling of heterogeneous materials, where the bounds are deduced by employing the principles of minimum complementary energy and minimum potential energy. Current random matrix based nonparametric approach is only adapted to the Wishart or matrix-variate Gamma probability model supported over the entire space of the symmetric positive-definite matrices, and therefore, unable to exploit additional information available through the lower and upper bounds when appropriate. Specifically, for a given material, the constitutive matrix is construed as a random matrix. A probability measure that reflects the constraints consistent with the energybased bounds, together with an associated sampling scheme, are constructed to synthesize realizations of this random matrix. An additional constraint of the ensemble mean matrix is also considered, and an appropriate probability model and sampling scheme are developed and illustrated numerically for this case.