In this paper we present an algorithm of the EM-type that performs fitting of PH-distributions with a given number of states to empirical distribution. In contrast to known approaches, the algorithm first generates a discretized representation of the observed data and then performs the fitting. In this way the algorithm is more efficient than other known EM-type methods. By an appropriate discretization, the main part as well as the tail of the distribution can be approximated sufficiently accurate.