使用最小二乘法擬合比較簡單：
x_r=[abscissa ones(size(abscissa))]\ordinates;
求出來即為題中的x和γ。如果不限方法，也可以使用多項式擬合：
p = polyfit(abscissa, ordinates,1);
得到的結果是壹致的（但二者分別是列向量和行向量）。
使用絕對值最小的擬合方法稍微復雜壹些：
e = ones(size(abscissa));
f = [0; 0; e];
A1 = [-abscissa -e -eye(length(abscissa))];
b1 = -ordinates;
A2 = [abscissa e -eye(length(abscissa))];
b2 = ordinates;
c = linprog(f,[A1;A2],[b1;b2]);
求出來的即為所需的系數。
完整的代碼如下（含繪圖）：
% Linear fit with least square and least absolute error
slope=3; intercept=-2;
abscissa = (-5:5)'; m = length(abscissa);
WhiteNoise = 5*randn(m,1);
ordinates = slope*abscissa + intercept + WhiteNoise;
GrossError=80;
ordinates(6)=ordinates(6)+GrossError;
ordinates(10)=ordinates(10)-GrossError;
x_r=[abscissa ones(size(abscissa))]\ordinates;
y2=[abscissa ones(size(abscissa))]*x_r;
e = ones(size(abscissa));
f = [0; 0; e];
A1 = [-abscissa -e -eye(length(abscissa))];
b1 = -ordinates;
A2 = [abscissa e -eye(length(abscissa))];
b2 = ordinates;
c = linprog(f,[A1;A2],[b1;b2]);
y1=[abscissa ones(size(abscissa))]*c(1:2);
plot(abscissa, ordinates, 'ro',abscissa, y2,'b-',abscissa, y1,'g--')
legend('Original data', 'Least L^2 error', 'Least L^1 error')
某次運行的結果如下（因數據為隨機生成，每次運行結果不同）：
對於L1最小擬合，也可以使用優化工具箱的函數fminunc來做：
objfun = @(x)sum( abs(abscissa*x(1)+x(2)-ordinates) );
c2=fminunc(objfun,[2 -2])
y3=[abscissa ones(size(abscissa))]*c2.';
plot(abscissa, ordinates, 'ro',abscissa, y2,'b-',abscissa, y1,'g--',abscissa, y3,'m:')
legend('Original data', 'Least L^2 error', 'Least L^1 error', 'fminunc')
擬合的結果可能與前述線性規劃的結果不同（但目標函數基本都能達到最小），這也說明這個問題本質上是多極值的，存在多個局部最優點。