Density function theory (DFT) calculations are integrated into a continuum mechanics formulation for simulating macroscale ferroelectric constitutive behavior. Theoretical relations associated with the Hellman-Feynman stress theory are used to guide the energy and electromechanical constitutive relations in the homogenized, continuum approximation of the stored energy function. Challenges associated with addressing uncertainty in length scales between the atomic lattice scale and the continuum mesoscale are treated using stochastic homogenization and Bayesian statistics. This is achieved by quantifying DFT energy, stress, and polarization over a range of constrained lattice dimensions and internal atomic configurations. The DFT energy calculations are fit to a continuum scale stored energy function using Markov chain Monte Carlo algorithms and Bayesian statistics. Probabilities distributions of the material parameters are used to guide the development of stochastic based homogenization of the stored energy. Ferroelectric hysteresis curves are numerically simulated on lead titanate and compared to data in the literature.