This paper deals with the construction of numerical solutions of random hyperbolic models with a finite degree of randomness that make manageable the computation of its expectation and variance. The approach is based on the combination of the random Fourier transforms, the random Gaussian quadratures and the Monte Carlo method. The recovery of the solution of the original random partial differential problem throughout the inverse integral transform allows its numerical approximation using Gaussian quadratures involving the evaluation of the solution of the random ordinary differential problem at certain concrete values, which are approximated using Monte Carlo method. Numerical experiments illustrating the numerical convergence of the method are included.