Observations of the neutral Hydrogen (HI ) 21-cm signal hold the potential of allowing us to map out the cosmological large scale structures (LSS) across the entire post-reionization era ($z leq 6$). Several experiments are planned to map the LSS over a large range of redshifts and angular scales, many of these targeting the Baryon Acoustic Oscillations. It is important to model the HI distribution in order to correctly predict the expected signal, and more so to correctly interpret the results after the signal is detected. In this paper we have carried out semi-numerical simulations to model the HI distribution and study the HI power spectrum $P_{HI}(k,z)$ across the redshift range $1 le z le 6$. We have modelled the HI bias as a complex quantity $tilde{b}(k,z)$ whose modulus squared $b^2(k,z)$ relates $P_{HI}(k,z)$ to the matter power spectrum $P(k,z)$, and whose real part $b_r(k,z)$ quantifies the cross-correlation between the HI and the matter distribution. We study the $z$ and $k$ dependence of the bias, and present polynomial fits which can be used to predict the bias across $0 le z le6$ and $0.01 le k le 10 , {rm Mpc}^{-1}$. We also present results for the stochasticity $r=b_r/b$ which is important for cross-correlation studies.