We present a method to automatically approximate momentbased invariants of probabilistic programs with non-polynomial updates of continuous state variables to accommodate more complex dynamics. Our approach leverages polynomial chaos expansion to approximate nonlinear functional updates as sums of orthogonal polynomials. We exploit this result to automatically estimate state-variable moments of all orders in Prob-solvable loops with non-polynomial updates. We showcase the accuracy of our estimation approach in several examples, such as the turning vehicle model and the Taylor rule in monetary policy.