In this paper the mathematical framework of an advanced algorithm is presented to efficiently form and solve the equations of motion of a multibody system involving uncertainty in the system parameters and/or the excitations. The uncertainty is introduced to the system through the application of the polynomial chaos expansion. In this scheme, states of the system, nondeterministic parameters, and the constraint loads are expanded using modal values as well as orthogonal basis functions. Computational complexity of the application of traditional methods to solve the stochastic equations of motion of a multibody system drastically grows as a cubic function of the number of the states of the system, uncertain parameters and the maximum degree of the polynomial chosen for the basis function. The presented method forms the equation of motion of the system without forming the entire mass and Jacobian matrices. In this strategy, the stochastic governing equations of motion of each individual body as well as the one associated with the kinematic constraint at the connecting joint are developed in terms of the basis functions and modal coordinates. Then sweeping the system in two passes assembly and disassembly, one can form and solve the stochastic equations of motion. In the assembly pass the non-deterministic equations of motion of the assemblies are obtained. In the disassembly process, these equations are then recursively solved for the modal values of the spatial accelerations and the constraints loads. In the serial and parallel implementations, computational complexity of the method increases as a linear and logarithmic functions of the number of the states of the system, uncertain variables, and the maximum degree of the basis functions used in the expansion.