A probabilistic method is developed to predict the uncertainty bounds on Frequency Response Functions (FRFs) developed from Finite Element models. A non-intrusive Polynomial Chaos Expansion (PCE) method is used to predict uncertainty regression models of the various parameters that make up a curvefit of the FRF: natural frequencies, damping ratios, complex amplitudes, mass and stiffness residuals, by making use of an efficient Latin Hypercube technique. These uncertainty models are then combined to efficiently determine PDFs of the parameters and also the uncertainty bounds of the FRFs. The approach is demonstrated using two examples; a simple beam containing uncertainty in Young's Modulus, and a full-scale aircraft composite wing model containing uncertainties in both Young's modulus and the shear modulus. The results were compared with Monte Carlo Simulation (MCS) and it was found that the parameter PDFs and FRF error bounds obtained using a 2nd-order PCE model agreed very well whilst requiring significantly less computation.