Closed form solutions for Bayesian models rarely exist, necessitating numerical approaches. In this paper, we propose and assess three different approaches based on polynomial approximations: sparse grid quadrature, quadratic simplex splines, and high degree multivariate polynomials. Each approach allows interpolating derivatives as well as function values, are conjugate allowing factoring and independent solving of larger models, and allow for easy integration and marginalization. Additionally, they allow for a high degree of parallelism, enabling GPU offloading. Accuracy vs time trade-offs are analyzed for a small dimensional model and compared to MCMC and QMC, and we asses the strengths and weaknesses of the approaches.