  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 \texttt{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 multi-dimensional approach to three science cases: (i) a
4-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 3-dimensional fit to the California Kepler Survey sample
showing how the planet radius valley evolves across different stellar masses;
and (iii) a 2-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 \texttt{MRExo}
\texttt{Python} package with these changes, which allows users to apply this
highly flexible framework to a variety of datasets beyond what we have shown
here.
