reserve estimation. Different factors including physical and structural characteristics of the deposit control the distribution of ore grade value. Metallic deposits and related systems have variable features. Various models has been proposed to describe their structural and spatial patterns. Generally, the data used in grade estimation are complicated, non-linear and of high variance and skewness. Traditional methods including geometrical and geostatistical ones often are not able to accurately handle the relationship among data. On the other hand, these methods need an expert with a lot of knowledge and also some assumption about the data which are not practically affordable. To improve the accuracy of grade estimation and also explore variability pattern hidden in data, complicated machine learning algorithms are used as an alternative estimation methods. In this study by comparing three selective machine learning methods the problem of improving ore grade estimation in the metallic ore deposits is considered and then a proper method to locate additional drillholes is proposed. The used machine learning methods are feed-forward neural networks, random forests and Guassian processes. After exploratory analysis of the data, multi set of training and test partitions based on drillholes are generated to appropriately compare the used estimation methods. The idea is a good estimator have a better estimation on drillholes that are far away from the training data. To overcome some difficulties of back propagation method, such as local search, a simultaneous perturbation artificial bee colony algorithm (SPABC) is introduced, for training the neural networks. SPABC, a variant of artificial bee colony algorithm (ABC), benefits from the local search of simultaneous perturbation stochastic approximation (A) algorithm and some new equations to improve the global and local search ability of original ABC algorithm. The performance of the proposed algorithm in training neural networks are compared with some other conventional optimization methods. A skewed activation function by introducing a new adjustment parameter is proposed to increase the flexibility and improve the accuracy of the fitted model. A local linear neural network structure is used to improve the generalization ability of the model. The generalization performance of the proposed network is compared with conventional MLP neural networks and other competent networks including wavelet neural networks. The ten-fold cross-validation method is used instead of traditional methodologies to avoid overfitting, explore appropriate number of hidden nodes and improving the generalization ability. Random forest regression was performed after optimizing number of random trees and leaf nodes. Gaussian process (GP) regression is the third machine learning methods which has some similarity to kriging methods. The main difference is related to the kernel parameters which are learned from the training data. Our study shows that the separate standardization of each location inputs may introduce spurious distortion in the geometery of the deposit. Therefore, a standardization method proposed to simultaneously standardize three spatial dimensions. An anisotropic kernel is proposed into the structure of GP to consider the effect of anisotropy and improve the grade estimation. The results of machine learning methods are validated using ordinary kriging and indicator kriging methods as representatives of linear and non-linear geostatistical methods, respectively. The normalized mean squared error as one of the main criterion of the performance evaluation is employed to compare the generalization ability of the used methods. The available data from drill holes of Esfordi phosphate and Sarcheshmeh copper deposit while having different nature of deposition, geological and structural setting are used to validate the performance of the proposed methods. The obtained results for Esfordi deposit show the closeness between the accuracy of GP with the provided specification and the proposed neural networks structure. However, GP is the best method for the Sarcheshmeh copper deposit. Another observation using three spatial dimensions of the deposit with rock types show significant improvement in the ore grade estimation. The estimated results of GP regression as the best estimation method is used to locate additional drill holes in the domains of mineralization. A new objective function is defined to consider the effect of geological composition of the deposit, alongside the estimated grade value and the estimation variance.