# Robust Fitting Iteratively Reweighted Least Squares

I'm interested in Mathematica's capabilities with regard to robust fitting vs. outliers. In the Statistical Model Analysis Tutorial under the Generalized Linear Models section heading, there is this sentence in a paragraph talking about Options for GeneralizedLinearModelFit:

Parameter estimates are obtained via iteratively reweighted least squares with weights obtained from the variance function of the assumed distribution.

This seems to say that GeneralizedLinearModelFit does iteratively reweighted least squares fitting. However, I can't see any options for iteratively reweighting data points, nor any indication that the weight applied to the points is anything other than 1.

For example:

$Version (*11.2.0 for Linux x86 (64-bit) (September 2, 2017)*) data = RandomVariate[NormalDistribution[100, 30], 200]; modData=ReplacePart[data,1->10^6]; Through@{Mean,StandardDeviation}@data (*{99.529,30.395}*) Through@{Mean,StandardDeviation}@modData (*{5098.99,70703.6}*) Through@{Median,SnDispersion}@modData (*{99.2765,31.6749}*) Through@{LinearModelFit,GeneralizedLinearModelFit}[data,{},x]//Normal (*{99.529,99.529}*) Through@{LinearModelFit,GeneralizedLinearModelFit}[modData,{},x]//Normal (*{5098.99,5098.99}*) GeneralizedLinearModelFit[modData,{},x,DispersionEstimatorFunction->Function[{y,yhat,w},SnDispersion[y-yhat]]]//Normal (*5098.99*) GeneralizedLinearModelFit[modData,{},x,DispersionEstimatorFunction->Function[{y,yhat,w},0]]//Normal (*5098.99*)  Clearly, the general method isn't reweighting according to the variance. The variance is gigantic in the modified data set and no reweighting occurred. Secondly, the dispersion estimator function appears to do nothing. Lastly, I thought about modifying the Weights function, but it only takes one original data point as its argument in functional form: Weights->func associates weight func[xi1,xi2,…,yi] with the i^(th) data element. The weight function, as documented, does not have access to the current iteration's prediction$\hat{y}_i\$ or at least the current best fit parameters/model function, which is something that would be needed for reweighting (see section 2.2 Table 1). Certainly this information could be provided via a side effect in the DownValues of a named function or a call to a cleverly constructed side function that happens to have the current state of the model, but given the construction of the options for Weight, the weight function is unlikely to be called each iteration of GeneralizedLinearModelFit anyway.

So, what am I misunderstanding? How do we get GeneralizedLinearModelFit to do iteratively reweighted least squares regression to ignore outliers?

Should we just repeatedly call the function and manually update the weights?

I figured out this can be done using the NormFunction option of FindFit. Using the Huber weighting function of the last reference in the question, the solution would be
modData=ReplacePart[RandomVariate[NormalDistribution[100, 30], 200],1->10^6];