The following program is a least mean square and Tukey's Biweight function written in octave. The progrem is not tested yet enouth. Please use it by your own risk.
weight for Tukey Biweight estimation function retval = weight(d, W) a = ( 1 - ((d/W) .* (d/W)) ); retval = min( max(sign(W - abs(d)), 0), a .* a); endfunction N = 3; x = [1 2 3]; y = [4 5 6]; W = 1; # least mean square sx = sum(x); sy = sum(y); sxx = sum(x .* x); sxy = sum(x .* y); a = ( (N * sxy) - (sx * sy) ) / ( (N * sxx) - (sx * sx) ); b = ( (sxx * sy) - (sxy* sx) ) / ( (N * sxx) - (sx * sx) ); # Tukey's Biweight estimation d = y - ( (a * x) + b ); w = weight(d, 1); # least mean square with Tukey's Biweight estimation s1 = sum(w); sx = sum(w .* x); sy = sum(w .* y); sxx = sum(w .* x .* x); sxy = sum(w .* x .* y); a = ( (s1 * sxy) - (sx * sy) ) / ( (s1 * sxx) - (sx * sx) ); b = ( (sxx * sy) - (sxy* sx) ) / ( (s1 * sxx) - (sx * sx) );