Material data can often contain multiple sources of uncertainty, such as microstructural variability that causes intrinsic uncertainty and measurement error, resulting in extrinsic uncertainty. Many existing methods rely on predefined model forms and assume that all uncertainties are lumped into a single source, leading to an imprecise characterization of material variability. This work presents a method that generates new model forms while separating and quantifying intrinsic and extrinsic uncertainty to generate improved material yield surface models. Our method relies on the machine learning technique of genetic programming-based symbolic regression (GPSR) to generate new yield surface models based on training data derived from finite element (FE) simulations. A multi-objective fitness function, called MOSR [1], is implemented that considers both an implicit fitness metric (ISR) and a statistical metric, the marginal likelihood, to select models that fit the data while also avoiding overfitting to noise in the dataset caused by extrinsic uncertainty. We also employ a multi-tree parent acyclic graph (MPA) approach [2] to help regularize the GPSR algorithm and ensure that models adhere to the known physics of the material. We then use a hierarchical Bayesian approach to first learn the extrinsic uncertainty by quantifying the added noise distributions for a single microstructure, thus eliminating any intrinsic uncertainty. We can then learn the intrinsic parameter distributions for the entire dataset, which consists of many microstructures. We implement this method to learn an improved yield surface model for porous metals, one that improves upon the more traditional Gurson-Tvergaard-Needleman (GTN) approach by incorporating random void distributions and matrix hardening in our FE models. The result is a yield surface model that is interpretable, as it can be easily compared to the GTN model and also contains parameter distributions that quantify the different sources of uncertainty present in the training data.