We study the randomized approximation of weakly singular integral operators. For a suitable class of kernels having a standard type of singularity and being otherwise of finite smoothness, we develop a Monte Carlo multilevel method, give convergence estimates and prove lower bounds which show the optimality of this method and establish the complexity. As an application we obtain optimal methods for and the complexity of randomized solution of the Poisson equation in simple domains, when the solution is sought on subdomains of arbitrary dimension.