Probabilistic Boolean networks (PBNs) is a widely used computational framework for modelling biological systems. The steady-state dynamics of PBNs is of special interest in the analysis of biological systems. However, obtaining the steady-state distributions for such systems poses a significant challenge due to the state space explosion problem which often arises in the case of large PBNs. The only viable way is to use statistical methods. We have considered the two-state Markov chain approach and the Skart method for the analysis of large PBNs in our previous work. However, the sample size required in both methods is often huge in the case of large PBNs and generating them is expensive in terms of computation time. Parallelising the sample generation is an ideal way to solve this issue. In this paper, we consider combining the German & Rubin method with either the two-state Markov chain approach or the Skart method for parallelisation. The first method can be used to run multiple independent Markov chains in parallel and to control their convergence to the steady-state while the other two methods can be used to determine the sample size required for computing the steady-state probability of states of interest. Experimental results show that our proposed combinations can reduce time cost of computing stead-state probabilities of large PBNs significantly.