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 \to2 $. The eigenvalues and eigenfunctions of these four operators are different, and the $ k $-th (for $ k \in {\mathbb N} $) eigenvalue of the spectral fractional Laplacian is always larger than those of the fractional Laplacian and regional fractional Laplacian. For any $ \alpha \in (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 $ {\mathcal O}(\delta ^{-\alpha }) $. In contrast, although the regional fractional Laplacian can be used to approximate the fractional Laplacian as $ \alpha \to2 $, it generally provides inconsistent result from that of the fractional Laplacian if $ \alpha \ll 2 $. Moreover, some conjectures are made from our numerical results, which could contribute to the mathematics analysis on these operators.