最小二乗法と Tukey の Biweight 法

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