Non-local viscoelasticity is a subject of great interest in the context of non-local theories. In a recent study, the authors have proposed a non-local fractional beam model where non-local effects are represented as viscoelastic long-range volume forces and moments, exchanged by non-adjacent beam segments depending on their relative motion, while local effects are modelled by elastic classical stress resultants. Long-range interactions have been given a fractional constitutive law, involving the Caputo's fractional derivative. This paper introduces a comprehensive numerical approach to calculate the stochastic response of the non-local fractional beam model under Gaussian white noise. The approach combines a finite-element discretization with a fractional-order state-variable expansion and a complex modal transformation to decouple the discretized equations of motion. While closed-form expressions are derived for the finite-element matrices associated with elastic and fractional terms, fractional calculus is used to solve the decoupled equations of motion, in both time and frequency domain. Remarkably, closed-form expressions are obtained for the power spectral density, cross power spectral density, variance and covariance of the beam response along the whole axis. Time-domain solutions are obtained by time-step numerical integration methods involving analytical expressions of impulse response functions. Numerical examples show versatility of the non-local fractional beam model as well as computational advantages of the proposed solution procedure.