This paper uses the method of Quadratures in conjunction with the Maximum Entropy principle to investigate the effect of parametric uncertainties on the mean power output and root mean square deflection of piezoelectric vibrational energy harvesting systems. Uncertainty in parameters of harvesters could arise from insufficient manufacturing controls or change in material properties over time. We investigate bimorph based harvesters that transduce ambient vibrations to electricity via the piezoelectric effect. Three varieties of energy harvesters-Linear, Nonlinear monostable and Nonlinear bistable are considered in this research. This analysis quantitatively shows the probability density function for the mean power and root mean square deflection as a function of the probability densities of the excitation frequency, excitation amplitude, initial deflection of the bimorph and magnet gap of the energy harvester. The method of Quadratures is used for numerically integrating functions by propagating weighted points from the domain and evaluating the integral as a weighted sum of the function values. In this paper, the method of Quadratures is used for evaluating central moments of the distributions of rms deflection and mean harvested power and , then, in conjunction with the principle of Maximum Entropy (MaxEnt) an optimal density function is obtained which maximizes the entropy and satisfies the moment constraints. The The computed nonlinear density functions are validated against Monte Carlo simulations thereby demonstrating the efficiency of the approach. Further, the Maximum Entropy principle is widely applicable to uncertainty quantification of a wide range of dynamic systems.