This paper presents the method of polynomial chaos expansion (PCE) for the forward kinematic analysis of nondeterministic multibody systems. Kinematic analysis of both open-loop and closed-loop systems are presented. 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 basis functions. This paper presents the detailed formulation of the kinematics of constrained multibody systems 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 and a SCARA robot. Also, the convergence of the PCE and Monte Carlo methods is analyzed in this paper. 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 computation 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 computational complexity.