最小二乗法と Tukey の Biweight 法

  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) );