In a fork-join network each incoming job is split into K tasks and the K tasks are simultaneously assigned to K parallel service stations for processing. For the distributions of response times and queue lengths of fork-join networks, no explicit formulae are available. Existing methods provide only analytic approximations for the response time and the queue length distributions. The accuracy of such approximations may be difficult to justify for some complicated fork-join networks. In this paper we propose a perfect simulation method based on coupling from the past to generate exact realisations from the equilibrium of fork-join networks. Using the simulated realisations, Monte Carlo estimates for the distributions of response times and queue lengths of fork-join networks are obtained. Comparisons of Monte Carlo estimates and theoretical approximations are also provided. The efficiency of the sampling algorithm is shown theoretically and via simulation.
