A new approach for estimation of population phylogeny together with a demographic model of population sizes and migration rates is presented. By including population phylogeny in a general Isolation-with-Migration (IM) model we address a long-standing problem in population genomics. The method utilizes a Markov chain Monte Carlo simulation over both genealogies and population phylogenies and provides a maximum posterior probability estimate of the population phylogeny by integrating over demographic history. By including an unsampled ghost population as an outgroup in the phylogeny, the method accounts for many kinds of interactions between sampled and unsampled populations. We applied the new method to genomic data from Human hunter gatherer populations from Africa. .