Robust M-estimator of location and dispersion by hand in R
<p>I try to estimate a location and dispersion model with R, as described <a href="http://www.wiley.com/WileyCDA/WileyTitle/productCd-0470010924.html" rel="nofollow">Maronna et al (2006, pp. 56)</a>. However, my estimate dispersion does not converge to the desired value. Do I have an error in the code?</p>

<pre><code>    #set.seed(12345)
    n &lt;- 100
    error&lt;-rnorm(n,0,1)
    x &lt;- 10+5*error

    # funcion bisquare
    psi_bisq &lt;- function(y){
      k &lt;- 2.2
      ifelse(abs(y)&lt;=k,ind &lt;- 1,ind &lt;- 0)
      return(y*ind*(1-(y/k)^2)^2)
    }

    # derivada de la funcion bisquare
    psi_bisq_der &lt;- function(y){
      k &lt;- 2.2
      ifelse(abs(y)&lt;=k,ind &lt;- 1,ind &lt;- 0)
      return(ind*(1-(y/k)^2)*(1-5*(y/k)^2))
    }

    # funcion rho_bisquare
    rho_bisq &lt;- function(y){
      k &lt;- 2.2
      ifelse(abs(y)&lt;=k, sal &lt;- 1-(1-(y/k)^2)^3, sal &lt;- 1)
      return(sal)
    }

    # segunda derivada de funcion rho_bisquare 
    rho_bisq_der2 &lt;- function(y){
      k &lt;- 2.2
      ifelse(abs(y)&lt;=k, sal &lt;- 6*(1-(y/k)^2)*(1-5*(y/k)^2)/k^2, sal &lt;- 0)
      return(sal)
    }

    # funcion de peso
    W &lt;- function(dato){
      ifelse(dato==0, ww &lt;- psi_bisq_der(0), ww &lt;- psi_bisq(dato)/dato)
      return(ww)
    }

    Wdis&lt;-function(dato){
      ifelse(dato==0,ww&lt;-rho_bisq_der2(0),ww&lt;-rho_bisq(dato)/dato^2)
    }

    #======== Estimador simultáneo de localización y dispersión 

    # estimador robusto de dispersi?n
    sigma_est_a &lt;- median(abs(x-median(x)))/0.6745

    # estimador robusto de la media, inicial
    mu_est_a &lt;- median(x)

    epsilon &lt;- .00001
    err1 &lt;- epsilon*sigma_est_a
    err2&lt;-epsilon
    error &lt;- 1000
    error2&lt;- 1000
    delta &lt;- .9
    iter&lt;-0

    while((error&gt;=err1) &amp;&amp; (error2&gt;=err2)){
      pesos_loc &lt;- rep(0,n)
      pesos_dis &lt;- rep(0,n)
      for(i in 1:n){
        r &lt;- (x[i]-mu_est_a)/sigma_est_a
        pesos_loc[i] &lt;- W(r)
        pesos_dis[i]&lt;-Wdis(r)
      }
      iter&lt;-iter+1
      mu_est_n &lt;-sum(pesos_loc*x)/sum(pesos_loc)
      sigma_est_n&lt;-sqrt(((sigma_est_a^2)/(n*delta))*sum(pesos_dis*r^2))
      error &lt;- abs(mu_est_n - mu_est_a)
      error2 &lt;- abs(sigma_est_n/sigma_est_a -1)
      mu_est_a&lt;-mu_est_n
      sigma_est_a&lt;-sigma_est_n
    }
    sigma_est_a
    mu_est_a
</code></pre>