MCMC is a great framework for sampling posterior distributions of complex statistical models. In order to sample from these distribution, the markov chains need to reach stationarity. This condition is often computationally intensive and (as far as I know) there is no way to fasten this stage with the inclusion of independent chains (that can be parallelised).

It could be, however, that once we reach stationarity, we can sample from independent chains simultaneously using threads or multicore processors or even in independent computers. The way to do this is described in the following example:

from pystan import StanModel
from pystan.constants import MAX_UINT
# bernoulli model
model_code = """
    data {
      int<lower=0> N;
       int<lower=0,upper=1> y[N];
        }
     parameters {
          real<lower=0,upper=1> theta;
        }
      model {
        theta ~ beta(0.5, 0.5);  // Jeffreys' prior
        for (n in 1:N)
           y[n] ~ bernoulli(theta);
        }
    """
data = dict(N=10, y=[0, 1, 0, 1, 0, 1, 0, 1, 1, 1])
sm = StanModel(model_code=model_code)
# initial seed can also be chosen by user
# MAX_UINT = 2147483647
seed = np.random.randint(0, MAX_UINT, size=1)
fit = sm.sampling(data=data, seed=seed)


# reuse tuned parameters
stepsize = fit.get_stepsize()
# by default .get_inv_metric returns a list
inv_metric = fit.get_inv_metric(as_dict=True)
init = fit.get_last_position()
# increment seed by 1
seed2 = seed + 1
control = {"stepsize" : stepsize,
         "inv_metric" : inv_metric,
         "adapt_engaged" : False
           }

fit2 = sm.sampling(data=data,
                  warmup=0,
                  iter=1000,
                  control=control,
                  init=init,
                  seed=seed2)

This example was taken from here .