The method of weighted least squares is shown to be an appropriate way of fitting variogram models. The weighting scheme automatically gives most weight to early lags and downweights those lags with a small number of pairs. Although weights are derived assuming the data are Gaussian, they are still appropriate when data are a smooth transform of the Gaussian case. The method of iterated generalized least squares, which accounts for correlation between variogram estimators at different lags, offers more statistical efficiency but at the cost of increased complexity. Weighted least squares for the robust estimator, based on square root differences, is less of a compromise. This article is written from a statistician's perspective, highlighting their role in geostatistical studies. The statistician typically leads the team through stages such as graphing and summarizing data, detecting nonstationarity, estimating spatial relationships, estimating in situ resources via kriging, and assessing recoverable reserves. Cressie (1984) presents a resistant approach to stages 1 and 2, using exploratory data analysis techniques adaptable to spatial data. Robust estimation of the variogram in the presence of contaminated data has been discussed by several authors. Here, the focus is on fitting a model to various variogram estimators, both classical and robust. Previous fitting procedures have been subjective or based on least squares, but statistical criteria can improve these methods by weighting the influence of different parts of the estimator. The intrinsic hypothesis assumes that the grade of an ore body is a realization of a random function. This hypothesis leads to the variogram, a crucial parameter in geostatistics. When data are nonstationary in the drift, naive variogram estimation can lead to bias. However, resistant techniques can reduce this bias. Nonresistant fitting of low-order polynomials in disjoint regions is another technique, but it still suffers from residual bias. The classical variogram estimator is based on data {Z_{t_i}, i=1,...,n}.The method of weighted least squares is shown to be an appropriate way of fitting variogram models. The weighting scheme automatically gives most weight to early lags and downweights those lags with a small number of pairs. Although weights are derived assuming the data are Gaussian, they are still appropriate when data are a smooth transform of the Gaussian case. The method of iterated generalized least squares, which accounts for correlation between variogram estimators at different lags, offers more statistical efficiency but at the cost of increased complexity. Weighted least squares for the robust estimator, based on square root differences, is less of a compromise. This article is written from a statistician's perspective, highlighting their role in geostatistical studies. The statistician typically leads the team through stages such as graphing and summarizing data, detecting nonstationarity, estimating spatial relationships, estimating in situ resources via kriging, and assessing recoverable reserves. Cressie (1984) presents a resistant approach to stages 1 and 2, using exploratory data analysis techniques adaptable to spatial data. Robust estimation of the variogram in the presence of contaminated data has been discussed by several authors. Here, the focus is on fitting a model to various variogram estimators, both classical and robust. Previous fitting procedures have been subjective or based on least squares, but statistical criteria can improve these methods by weighting the influence of different parts of the estimator. The intrinsic hypothesis assumes that the grade of an ore body is a realization of a random function. This hypothesis leads to the variogram, a crucial parameter in geostatistics. When data are nonstationary in the drift, naive variogram estimation can lead to bias. However, resistant techniques can reduce this bias. Nonresistant fitting of low-order polynomials in disjoint regions is another technique, but it still suffers from residual bias. The classical variogram estimator is based on data {Z_{t_i}, i=1,...,n}.