Nuclear fusion reactor design crucially depends on numerical simulation. The plasma can usually be modeled using fluid equations (for mass, momentum and energy). However, the reactor also contains neutral (non-charged) particles (which are important in its operation), of which both the position and velocity distribution is important. This leads to a Boltzmann-type transport equation that needs to be discretised with a Monte Carlo method. One then obtains a coupled finite-volume/Monte-Carlo simulation, of which the results possess both a bias and a variance. In this talk, I introduce the problems associated with the simulation of the plasma edge region in a fusion reactor. I discuss how to couple a finite volume discretisation of the plasma equations with a Monte Carlo simulation of the neutral particles, and show how the Monte Carlo errors affect convergence of steady state computations and reliability of gradient computations (necessary during optimization).