The time-varying joint probability distributions of interest for realistic power systems simulation have high dimensional supports, and trajectory-level nonlinearities induce non-Gaussianity, thereby requiring non-parametric computation. For example, transient stability analysis in the presence of stochastic renewables, and uncertainties in the initial conditions and parameters requires scalable yet rigorous predictive algorithms that do not suffer from the "curse-of-dimensionality". We will present a proximal algorithm that implements a variational recursion on the space of joint probability measures to propagate the stochastic uncertainties in power system dynamics over high-dimensional state space. The proposed algorithm takes advantage of the exact nonlinearity structures in the trajectory-level dynamics of the networked power systems, and shows that by lifting the dynamics in the space of measures, it is possible to design a scalable algorithm that obviates the discretization of the high dimensional state space or function approximation. This is enabled by certain generalized gradient descent structure on the infinite-dimensional manifold of probability measures induced by the networked power system dynamics. Somewhat surprisingly, we show that this measure-valued gradient descent structure exists and can be algorithmically leveraged, even though the underlying sample path dynamics remains mixed-conservative dissipative. We will provide the theoretical details, convergence guarantees, and numerical examples on realistic test systems.