Bulk density (BD) is a key soil indicator for the assessment of the sustainability of production from agricultural land. However, it is difficult to obtain data on BD from field studies because sampling is an extremely labor-intensive, time-consuming and expensive task at the large scale needed to generate meaningful amounts of data. Therefore, the purpose of this study was to evaluate different pedotransfer functions (PTFs) for their usefulness in the estimation of the BD values of an intensely cultivated area through the use of different geostatistical methods. PTFs are estimation functions of certain soil properties derived from raw data for other soil properties. Different approaches have been used for the development of PTFs, including regression methods and artificial neural networks. In this study, the multivariate analysis methods, general and stepwise regression, were used for BD estimation. In addition, different learning algorithms were evaluated with artificial neural networks (ANN) for their efficacy in estimating the BD. Seven basic soil properties, namely, sand, silt, clay, organic matter, pH, EC and CaCO₃, were used in the development of models. The estimation power of the general regression model (normalized root mean square error (NRMSE): 7.10%) was higher than that of stepwise regression (NRMSE: 19.93%). Additionally, the lowest NRMSE (6.74%) and the highest R² for BD estimation determined with different learning algorithms through artificial neural networks (ANN) were obtained with the Levenberg-Marguardt algorithm. Moreover, for BD estimation, ANN performed better than the multivariate regression equations. Spatial distribution maps of soil BD were generated with the commonly used ordinary kriging method by utilizing the real BD values in combination with BD values obtained from estimation models. Maps of BD values produced by stepwise regression estimation deviated significantly from maps generated with real values whereas ANN-II (Bayesian regularization algorithm) values were closest to the real values and that was reflected in the increased accuracy of mapping.