Reliable equilibrium reconstruction is indispensable for understanding and controlling hot magnetized plasmas to achieve fusion reactors. In axisymmetric systems, current and pressure profiles that satisfy the force balance conditions are given by the Grad-Shafranov (GS) equation. While many novel approaches have been developed to swiftly and robustly find an optimum solution of the GS equation, approaches based on a single solution search may not be adaptable if diagnostics fail to provide sufficient constraints. Here, we investigate the solution space of the GS equation when only basic edge magnetic measurements are available. By combining Gaussian process regression and Markov chain Monte Carlo sampling within the Bayesian framework, we treat each current element as an independent variable and evaluate the probability distribution that describes all possible solutions. We have applied this inference frame to the geometry of the PLATO tokamak and shown that the flux surface locations can be determined relatively well only from 16 pick-up coils, 4 flux loops and a diamagnetic loop. On the other hand, the toroidal current density is inferred with limited success, and the inferences of the safety factor and pressure profiles are difficult. The characterization of possible choices of equilibria realized by this inference frame will help optimize diagnostic setups for equilibrium reconstruction.