We present a numerical method for utilizing stochastic models with differing fidelities to approximate parameterized functions. A representative case is where a high-fidelity and a low-fidelity model are available. The low-fidelity model represents a coarse and rather crude approximation to the underlying physical system. However, it is easy to compute and consumes little simulation time. On the other hand, the high-fidelity model is a much more accurate representation of the physics but can be highly time consuming to simulate. Our approach is nonintrusive and is therefore applicable to stochastic collocation settings where the parameters are random variables. We provide sufficient conditions for convergence of the method, and present several examples that are of practical interest, including multifidelity approximations and dimensionality reduction.