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 weighting scheme, while there can be many solutions 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 and convergence 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 problem of association of different types of agents (person-households) is then treated as a maximum weight problem of a bipartite graph. 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.64).