Microsimulation of urban systems evolution requires synthetic population as a key input. Currently, the focus is on treating synthesis as a fitting problem and thus various techniques have been developed, including Iterative Proportional Fitting (IPF) and Combinatorial Optimization based techniques. The key shortcomings of fitting based procedures include: a) synthesis of only one population, while there can be many populations matching the available data b) due to cloning rather than true synthesis of the population, losing the heterogeneity that may not have been captured in the microdata c) over reliance on the accuracy of the data to determine the cloning weights d) poor scalability with respect to the increase in number of attributes of the synthesized agents. In order to overcome these shortcomings, we propose a Markov Chain Monte Carlo (MCMC) simulation based approach. Partial views of the joint distribution of agents attributes that are available from various data sources can be used to simulate draws from the original distribution. The real population from Swiss census is used to compare the performance of simulation based synthesis with the standard IPF. The standard root mean square error statistics indicated that even the worst case simulation based synthesis (SRMSE = 0.35) outperformed the best case IPF synthesis (SRMSE=0.85). We also used this methodology to generate the synthetic population for Brussels, Belgium where the data availability was highly limited.