Modeling genetic regulatory networks is an important problem in genomic research. Boolean Networks (BNs) and their extensions Probabilistic Boolean Networks (PBNs) have been proposed for modeling genetic regulatory interactions. In a PBN, its steady-state distribution gives very important information about the long-run behavior of the whole network. However, one is also interested in system synthesis which requires the construction of networks. The inverse problem is ill-posed and challenging, as there may be many networks or no network having the given properties, and the size of the problem is huge. The construction of PBNs from a given transition-probability matrix and a given set of BNs is an inverse problem of huge size. We propose a maximum entropy approach for the above problem. Newton’s method in conjunction with the Conjugate Gradient (CG) method is then applied to solving the inverse problem. We investigate the convergence rate of the proposed method. Numerical examples are also given to demonstrate the effectiveness of our proposed method.