Nuclear parton distribution functions (nPDFs) are essential in understanding nuclear structure and predicting outcomes in heavy-ion collisions. Global QCD analyses, a statistical method that relies on fitting the nPDF-dependent theory predictions to the relevant experimental data, have been the primary method for extracting nPDFs. One of the challenges in these analyses is the estimation of uncertainties in the results. Typically, the Hessian method has been applied in order to propagate experimental uncertainties into the predictions for nuclei collisions. However, Hessian does not always provide reliable results due to the nature of nPDF fits—limited data constraints, non-gaussianity, and the potential for multiple minima in the parameter space. In this study, we introduce the application of Markov Chain Monte Carlo (MCMC) methods as a statistically sophisticated alternative for the estimation of nPDF uncertainties by direct sampling from the probability distribution of the nPDF parameters. This method aims to approximate the posterior distribution of the nPDF parameters based on the prior distribution and the likelihood of the observed data. The process iteratively generates a series of samples that converge to this posterior distribution, which integrates all information about the parameters after considering the observed data. This allows for comprehensive statistical inference regarding the parameters' values, reflecting their uncertainties accurately.