Estimation of distribution programming is an iterative method to evolve program trees. It estimates the distribution of the most suitable program trees and then produces a new generation of program trees by sampling from the distribution. This paper proposes a hierarchical estimation of multimodal distribution programming (HEMMDP). First, the population is divided into K subpopulations by a clustering algorithm where the distribution of each subpopulation is modified according to an objective function. Then, at each generation, a new subpopulation is generated from the modified distribution. The objective function aims to gradually improve the fitness of the program trees in each subpopulation. Finally, the appropriate program trees are added as new terminal nodes to the terminal set, resulting in a new hierarchy. The best-fitting program trees from each subpopulation with high synergistic value are chosen as basis functions. The proposed approach uses a linear function of the basis functions to solve the regression problem. The proposed method is evaluated on several real-world benchmark datasets. The datasets are divided into four classes: small-difficult, small-easy, large-difficult, and large-easy. The proposed method improves the results of the best methods for the regression problem by 232% and 62% for small difficult data sets and large difficult data sets, respectively.