Nuclear Parton Distribution Functions (nPDFs) are critical for understanding nuclear structure and making heavy-ion collision predictions. nPDFs have been determined via ‘global QCD analyses’, which is a statistical approach based on fitting nPDF-dependent theoretical predictions to the relevant experimental data. One of the crucial aspects of nPDF determination is uncertainty estimation. Typically, the Hessian method is used to propagate experimental uncertainties into predictions for collisions of nuclei. However, due to the nature of nPDF fits (such as limited data constraints, non-gaussianity, and possible multiple minima), this method does not always provide reliable results. In this work, we introduce the application of Markov Chain Monte Carlo (MCMC) methods as a statistically sophisticated alternative for estimating nPDF uncertainties by sampling directly from the probability distribution of the nPDF parameters.