This study assessed the parameter recovery accuracy of Marginal Maximum Likelihood (MML) and two Markov Chain Monte Carlo (MCMC) methods, Gibbs and Hamiltonian Monte Carlo (HMC), under the four-parameter unidimensional binary item response function. Data were simulated under the mixed factorial design with sample size (1,000; 2,500; and 5,000 respondents) and latent trait distribution (normal and negatively skewed) as the between-subjects factors, and estimation method (MML, Gibbs, and HMC) as the within-subjects factor. Results indicated that in general, MML was more heavily impacted by latent trait skewness, but MML also improved its performance more strongly than MCMC when sample size increased. Two MCMC methods remained advantageous with lower root mean square errors (RMSE) of item parameter recovery across all conditions under investigation, but sample size increase brought a correspondingly narrower gap between MML and MCMC regardless of theta distributions. Gibbs and HMC provided nearly identical outcomes across all conditions, and no considerable difference between these two MCMC methods was detected. Sample size and latent trait distribution had little observable effect on trait score estimation by MCMC and Expected a Posteriori following MML (MML-EAP), which were essentially unbiased and had similar RMSE across all conditions. Discussions of the findings and model calibration issues are presented together with practical implications and future research recommendations.