System reliability is often estimated by the use of components' reliability test results when system test data are not available, or are very scarce. A method is proposed for computing the exact posterior probability density function, cumulative distribution function, and credible intervals for system reliability in a Bayesian setting, with the use of components' prior probability distributions and current test results. The method can be applied to series, parallel, and many mixed systems. Although in theory the method involves evaluating infinite series, numerical results show that a small number of terms from the infinite series are sufficient in practice to provide accurate estimates of system reliability. Furthermore, because the coefficients in the series follow some recurrence relations, our results allow us to calculate the reliability distribution of a large system from that of its subsystems. Error bounds associated with the proposed method are also given. Numerical comparisons with other existing approaches show that the proposed method is efficient and accurate.