A few things to note (I could be miss-interpreting your problem so please correct me if I am taking this in the wrong direction): Define your prior (before seeing any data) as $\pi_0(\theta)$. Define a first data set as $Y_1$ and a second data set as $Y_2$. Then define the first posterior $$ \pi_1(\theta|Y_1) = \frac{f(Y_1|\theta)\pi_0(\theta)}{m(Y_1)} $$ You can simulate from $\pi_1(\theta|Y_1)$ with STAN, pyMC, etc. (the tools you mention) To restate what you said in your question, you want to use $\pi_1(\theta|Y_1)$ as the prior for $\pi_2(\theta|Y_2)$ i.e. $$ \pi_2(\theta|Y_2) = \frac{f(Y_2|\theta)\pi_1(\theta|Y_1)}{m(Y_2)} $$ but notice $$ \frac{f(Y_2|\theta)\pi_1(\theta|Y_1)}{m(Y_2)}= \frac{f(Y_2|\theta)f(Y_1|\theta)\pi_0(\theta)}{m(Y_1)m(Y_2)} $$ So under the assumption $Y_1$ and $Y_2$ are independent (this generally goes along with iid assumptions, and if it wasn't true you would still use the final result below) $$ \pi_2(\theta|Y_2)=\pi_2(\theta|Y_2,Y_1) = \frac{f(Y_2,Y_1|\theta)\pi_0(\theta)}{m(Y_2,Y_1)} $$ Which can be simulated using STAN, pyMC, etc. **My point: You don't *have* to empirically estimate the prior with the previous posterior, You can just run the model over again with all the data** I understand that the above recommendation may not be possible in some situations, it could be computationally costly or the data may not be in a convenient form it. In these situations you can use a [self-normalized importance sampler][1] (which can be implemented in Python, R, or something similar). So let's back-up and suppose you drew $\theta^{(1)},...,\theta^{(G)} \sim \pi_1(\theta|Y_1)$. Now you want to sample from $\pi_2(\theta|Y_2,Y_1)$. To do this you can re-weight $\theta^{(1)},...,\theta^{(G)}$ with the importance weights $$ w^{(g)} = \frac{f(Y_2|\theta)}{\sum_{g=1}^G f(Y_2|\theta) } $$ so they effectively become a sample from $\pi_2(\theta|Y_2,Y_1)$. For importance sampling to work, $\pi_1(\theta|Y_1)$, and $\pi_2(\theta|Y_2,Y_1)$ must be similar. If not, the variance of the sampler will be very high. You can approximate the "effective sample size" of the importance sampler via; $$ G.\mathrm{eff} \approx \frac{1}{\sum_{g=1}^G [w^{(g)}]^2} $$ There are other related techniques that reduce the variance of the procedure (i.e. sequential importance sampling/particle filtering and rejection sampling). It would be a good idea for you to research these algorithms on your own and decide for yourself how you would like to use them if you are interested in this type pof approach. Also see @Xian answer to my own question http://stats.stackexchange.com/questions/192088/does-this-monte-carlo-technique-have-a-name ...he knows more about this stuff than I do. [1]: http://statweb.stanford.edu/~owen/mc/Ch-var-is.pdf