The demand for efficient energy harvesting from vibratory environments has notably increased with the proliferation of piezoelectric devices. However, mismatches between model predictions and experimental observations underpin a critical challenge: accurately modeling piezoelectric energy harvesters (PEHs) while accounting for the inherent uncertainties in their electromechanical properties. A comprehensive Bayesian inference framework is posited for identifying electromechanical properties and facilitating model class selection of nonlinear PEHs. This approach addresses uncertainties by integrating prior knowledge and experimental data to update model parameters, offering the potential to balance model complexity and predictive precision systematically. The influence of different excitation amplitudes and the number of observations employed in the framework is studied. The results demonstrated the ability to capture both linear and nonlinear behaviors, while the model class selection effectively determined the simplest constitutive models sufficient for accurate predictions, validated against varying excitation intensities. By elucidating the trade-off between model simplicity and accuracy, this method not only refines the understanding of PEH behavior under diverse operating conditions but also equips designers with a robust predictive tool. Consequently, this methodology serves to optimize the deployment of PEHs in real-world applications, such as structural health monitoring, ensuring reliable energy supply with quantifiable confidence levels.