Consider a function \(\widetilde f\) defined on the real line whose values are fuzzy numbers. The fuzzy Riemann integral of \(\widetilde f\) is a fuzzy set whose membership function is \(\xi(r)= \sup_{0\leq\alpha\leq 1}\alpha 1_{A_\alpha}(r)\), where \[ A_\alpha= \Biggl[\int^b_a\widetilde f^L_\alpha(x) dx,\;\int^b_a\widetilde f^R_\alpha(x) dx\Biggr] \] and \([\widetilde f^L_\alpha(x),\widetilde f^R_\alpha(x)]\) is the \(\alpha\)-cut of \(\widetilde f(x)\). After surveying several properties of these integrals, the author shows how to compute them numerically using the Monte Carlo method. The convergence properties are derived using the strong law of large numbers for fuzzy random variables. The membership function of the integral is then evaluated by transforming it into a standard optimization problem.