金子邦彦研究室プログラミングOctave の活用最小二乗法と Tukey の Biweight 推定法

最小二乗法と 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) );