In this paper, we effectively solve the inverse source problem of the fractional Poisson equation using MC-fPINNs. We construct two neural networks $u_{NN} (x; θ)$ and $f_{NN} (x;ψ)$ to approximate the solution $u^∗(x)$ and the forcing term $f^∗(x)$ of the fractional Poisson equation. To optimize these networks, we use the Monte Carlo sampling method and define a new loss function combining the measurement data and underlying physical model. Meanwhile, we present a comprehensive error analysis for this method, along with a prior rule to select the appropriate parameters of neural networks. Numerical examples demonstrate the great accuracy and robustness of the method in solving high-dimensional problems up to 10D, with various fractional orders and noise levels of the measurement data ranging from 1% to 10%.