<p><b>Abstract</b>—Until recently, the model of the interior structure of a planet was specified based on the solution of the direct problem with input data on the planetary gravitational field (mass, moment of inertia, tidal Love numbers <i>k</i>₂) and the presumed geochemical composition of the planet. To reconcile the different model parameters with the observed quantities, it is important to solve the inverse problem. One of the goals of this study is to design and implement a computational algorithm that allows for easy and fast addition of new input data. At the first step, a computational algorithm is constructed to determine the radial distributions of the parameters of the planet’s interior from a set of observational data. Using the Bayesian statistics approach, we then formulate the inverse problem and solve it using the Markov chain Monte Carlo (MCMC) method. The probabilistic approach to solving the inverse problem greatly simplifies the matching of model parameters that satisfy the observations and the a priori data. The Bayesian statistics approach allows us to take into account the correspondence between the initial information about the model and the observed data. The developed computational algorithm was tested on the classical model example of gravity data inversion. The results of the numerical experiment are presented graphically. The algorithm for solving the problem has the peculiarity that each Markov chain is computed completely independently of the others. The problem is easily distributed evenly over all the cores of a computer or a cluster. This greatly reduces the running time of the computational algorithm, which is important in the future when the number of input parameters increases. At the second step of the work, it is planned to use the presented computational algorithm to find parameter distributions in the interior of planets from the known observational data.</p>