In this paper, we introduce a bivariate distribution on (Formula presented.) arising from a single underlying Markov jump process. The marginal distributions are phase-type and discrete phase-type distributed, respectively, which allow for flexible behavior for modeling purposes. We show that the distribution is dense in the class of distributions on (Formula presented.) and derive some of its main properties, all explicit in terms of matrix calculus. Furthermore, we develop an effective EM algorithm for the statistical estimation of the distribution parameters. In the last part of the paper, we apply our methodology to an insurance dataset, where we model the number of claims and the mean claim sizes of policyholders, which is seen to perform favorably. An additional consequence of the latter analysis is that the total loss size in the entire portfolio is captured substantially better than with independent phase-type models.