  We compute two- and three-point functions at criticality for the
three-dimensional Ising universality class. To this end we simulate the
improved Blume-Capel model at the critical temperature on lattices of a linear
size up to $L=1600$. As check also simulations of the spin-1/2 Ising model are
performed. We find $f_{\sigma \sigma \epsilon} = 1.051(1)$ and $f_{\epsilon
\epsilon \epsilon} =1.533(5)$ for operator product expansion coefficients.
These results are consistent with but less precise than those recently obtained
by using the bootstrap method. An important ingredient in our simulations is a
variance reduced estimator of $N$-point functions. Finite size corrections
vanish with $L^{-\Delta_{\epsilon}}$, where $L$ is the linear size of the
lattice and $\Delta_{\epsilon}$ is the scaling dimension of the leading
$Z_2$-even scalar $\epsilon$.
