The importance of modeling geological units and mineralization factor, as a key step in evaluating mineral deposits, is undeniable. Multivariate geostatistical methods are very useful tools in the case of geochemical modeling. In this study, the minimum/maximum autocorrelation factor (MAF), a multivariate geostatistical method, and the sequential indicator simulation (SIS) were used in order to model the Sury Gunay epithermal gold deposit, NW of Iran. The analyses were done using core samples obtained from drill holes. The MAF is a geospatially modified version of principal component analysis, which decorrelates factors for all lag distances. By applying the MAF on alr-transformed data, four MAF were obtained and the fourth factor represents the mineralization factor. Joint simulation of the mineralization elements was carried out by applying the sequential Gaussian simulation to the fourth factor scores. The main rock types and alterations of the deposit were also modeled using SIS and the probabilistic models of the rock types, and alterations were obtained using E-type maps resulting from 20 realizations. The results indicate that there are better relations with higher values of the mineralization factor and the existence of the volcanogenic breccia and dacite porphyry rock types. Furthermore, the simulated alterations demonstrate that the higher probability of silicification existence may have better correlation with higher concentrations of the mineralization elements.