In this paper, we study four nonlocal diffusion operators, including the fractional Laplacian, spectral fractional Laplacian, regional fractional Laplacian, and peridynamic operator. These operators represent the infinitesimal generators of different stochastic processes, and especially their differences on a bounded domain are significant. We provide extensive numerical experiments to understand and compare their differences. We find that these four operators collapse to the classical Laplace operator as alpha -> 2. The eigenvalues and eigenfunctions of these four operators are different, and the k-th (for k is an element of N) eigenvalue of the spectral fractional Laplacian is always larger than those of the fractional Laplacian and regional fractional Laplacian. For any alpha is an element of (0, 2), the peridynamic operator can provide a good approximation to the fractional Laplacian, if the horizon size delta is sufficiently large. We find that the solution of the peridynamic model converges to that of the fractional Laplacian model at a rate of O(delta(-alpha)). In contrast, although the regional fractional Laplacian can be used to approximate the fractional Laplacian as alpha -> 2, it generally provides inconsistent result from that of the fractional Laplacian if alpha << 2. Moreover, some conjectures are made from our numerical results, which could contribute to the mathematics analysis on these operators.