Fundamental to our understanding of planetary bulk compositions is the relationship between their masses and radii, two properties that are often not simultaneously known for most exoplanets. However, while many previous studies have modeled the two-dimensional relationship between planetary mass and radii, this approach largely ignores the dependencies on other properties that may have influenced the formation and evolution of the planets. In this work, we extend the existing nonparametric and probabilistic framework of MRExo to jointly model distributions beyond two dimensions. Our updated framework can now simultaneously model up to four observables, while also incorporating asymmetric measurement uncertainties and upper limits in the data. We showcase the potential of this multidimensional approach to three science cases: (i) a four-dimensional joint fit to planetary mass, radius, insolation, and stellar mass, hinting of changes in planetary bulk density across insolation and stellar mass; (ii) a three-dimensional fit to the California Kepler Survey sample showing how the planet radius valley evolves across different stellar masses; and (iii) a two-dimensional fit to a sample of Class-II protoplanetary disks in Lupus while incorporating the upper limits in dust mass measurements. In addition, we employ bootstrap and Monte Carlo sampling to quantify the impact of the finite sample size as well as measurement uncertainties on the predicted quantities. We update our existing open-source user-friendly MRExo Python package with these changes, which allows users to apply this highly flexible framework to a variety of data sets beyond what we have shown here.