It is an interesting topic to predict how dynamics change in disease progression. Our goal is to predict dynamics of the vital cells in the human immune system using the hierarchical HIV models. We assumed that the parameters are random variables with intra-level noise and inter-level noise. Individual-level parameters of HIV model are estimated by generalized least squares method based on the field data. At the population-level, it is not desirable to estimate joint distribution of all the parameters because of computational cost and overfitting. Thus, we build a model by estimating the distribution of only a few selected by the parameter reduction using the sensitivity matrices. Then, we estimate the joint probability distribution of remaining parameters employing algorithms for population parameter estimation. Finally, we solve ordinary differential equations with random coefficients to predict the dynamics of target cells.