In this paper we present a method for estimation of functionals depending on one or several phase-type distributions. This could for example be the ruin probability in a risk reserve process where claims are assumed to be of phase-type. The proposed method uses a Markov chain Monte Carlo simulation to reconstruct the Markov jump processes underlying the phase-type variables in combination with Gibbs sampling to obtain a stationary sequence of phase-type probability measures from the posterior distribution of these given the observations. This enables us to find quantiles of posterior distributions of functionals of interest, thereby representing estimation uncertainty in a flexible way. We compare our estimates to those obtained by the method of maximum likelihood and find a good agreement. We illustrate the statistical potential of the method by estimating ruin probabilities in simulated examples.