This paper presents the method of polynomial chaos expansion (PCE) for the forward kinematic analysis of nondeterministic multibody systems with kinematically closed-loops. The PCE provides an efficient mathematical framework to introduce uncertainty to the system. This is accomplished by compactly projecting each stochastic response output and random input onto the space of appropriate independent orthogonal polynomial base functions. This paper presents the detailed formulation of the kinematics of a constrained multibody system at the position, velocity, and acceleration levels in the PCE scheme. This analysis is performed by projecting the governing kinematic constraint equations of the system onto the space of appropriate polynomial base functions. Furthermore, forward kinematic analysis is conducted at the position, velocity, and acceleration levels for a non-deterministic four-bar mechanism with single and multiple uncertain parameters in the length of linkages of the system. Time efficiency and accuracy of the intrusive PCE approach are compared with the traditionally used Monte Carlo method. The results demonstrate the drastic increase in the computational time of Monte Carlo method when analyzing complex systems with a large number of uncertain parameters while the intrusive PCE provides better accuracy with much less computation complexity.