Pointers to understand "rlm" in R better? Could you please point me to tutorial/notes that can help me understand "rlm" better? 
Here is an example: summary(rlm(stack.loss ~ ., stackloss, weights=myweights))
My questions are:
a. How does the iterative procedure work?
   b. What's the relation between the weights I supplied above "myweights" and the final used weights?
My understanding is that the final weights used are actually the Huber loss weights? So are myweights not used at all? Ultimately, I want to use myweights... I have a feeling that actually I should put "myweights" into the "w" argument not the "weights" argument...
c. Could you please give an example showing whether the final used weights are actually the Huber weights or my weights?
I am also confused about the manual:
The input arguments:

wt.method are the weights case weights (giving the relative importance of case, so a  weight of 2 means there are two of these) or the inverse of the variances, so a weight of two means this error is half as variable? 
w (optional) initial down-weighting for each case. 
init (optional) initial values for the coefficients OR a method to find initial values OR the result of a fit with a coef component. Known methods are "ls" (the default) for an initial least-squares fit using weights w*weights, and "lts" for an unweighted least-trimmed squares fit with 200 samples. 
The returned values:
w the weights used in the IWLS process 
wresid a working residual, weighted for "inv.var" weights only. 

Anybody please shed some light?
Thank you! 
 A: *

*a. How does the iterative procedure work? 


It's best explained here (first two equations). 
The key elements are the starting $\beta^0\in\mathbb{R}^{p+1}$'s.
 In this case (rlm), the algorithm uses a numbers ($M$) of starting 
$\beta^0$'s $\{\beta^{0,1},...,\beta^{0,M}\}$ and carries 
the iterative procedure on each of them until convergence. 
Normally, each of these $\beta^{0,m}$ should be the slope
 of the regression line passing by $p+1$ observations 
drawn at random from the $n$ you have (denoting $n\times p$
 the dimensions of your design matrix). But in the MASS implementation this is 
only true if init is set to lts. Using the default init=ls uses as starting point a single ($M=1$) weighted lm fit (with weights w or if none is supplied a single, un-weighted, lm fit). From a robustness perspective the choice of this default behavior is hard to comprehend but consistent with the generally low quality of the implementation of the robust functions in MASS. 
--if you are looking for good implementation of robust procedures, i would advise having a look at the dedicated task view--
When $M>1$, the final $\beta^{F,*}$ shown is chosen among all the $\beta^{F,m}$ 
as the one which has smallest value of the sum of the $n/2$ smallest normalized residuals.


*

*b. What's the relation
between the weights I supplied above "myweights" and the final used
weights? 


At each iteration (and therefore also at the final one) the weights used in the "myweights"-weighted IRLS are equal to "myweights" times the normal IRLS weights. 


*

*c. Could you please give an example showing whether the
final used weights are actually the Huber weights or my weights?


Neither: a product of both. 
