以下存为一个MATLAB文件t_estimator40lm(X,Y,v)，自变量矩阵X,因变量列向量Y，参数v(一般取4左右)必须先在主程序中输入，再调用该子函数。

function [beta sigma fvalue] = t_estimator40lm(X,Y,v)
% 普通线性模型T-型估计的EM算法
%%%初始化%%%%%%
[n,p] = size(X);
beta0 = ladreg0lm(X,Y) ;% 用0截距线性模型L1估计作初值
sigma0 = mad(Y- X*beta0)/0.6745;
%Y = MAD(X) returns the mean absolute deviation of the values in X.  For vector input, Y is MEAN(ABS(X-MEAN(X)).
w0 = ones(n,1);
w1 = (v+1)./(v+w0.^2.*(Y- X*beta0).^2/sigma0^2);
k = 1;
for ii = 1:1000
    T1 = zeros(p);
    for i = 1:n
        T1=T1+w1(i)*w0(i)^2*X(i,:)'*X(i,:);
    end
    T2 = zeros(p,1);
    for i = 1:n
        T2=T2+w1(i)*w0(i)^2*X(i,:)'*Y(i,:);
    end
    beta1 = pinv(T1)*T2;
    T3 = 0;
    for i = 1:n
        T3 = T3 + w1(i)*w0(i)^2*(Y(i,:) -X(i,:)*beta1)^2;
    end
    sigma1 =  sqrt(T3/n);
    w2 = (v+1)./(v+w0.^2.*(Y- X*beta1).^2/sigma1^2);
    if max(abs(beta1-beta0))>1e-6
        beta0 = beta1;    sigma0 = sigma1;   w1=w2; 
    else
        fprintf('Iteration number k = %d ',k);
        break
    end
    k = k+1;
end
beta = beta0; sigma = sigma0;
canshu = [beta0;sigma0];
fvalue = goalls(X,Y,w0,v,canshu);

%%%子程序
function f = goalls(X,Y,w,v,canshu)%无常数项
[n p]= size(X);
beta = canshu(1:p,:); sigma = canshu(p+1,:);
f = 0;
for i = 1:n
    f = f +rhols(w(i)*(Y(i)-X(i,:)*beta)/sigma,v) +log(sigma);
end

function f = rhols(x,v)
f = (v+1)/2*log(1+x^2/v);

function [beta L1norm]= ladreg0lm(X,Y)  %用线性规划求零截距线性模型L1估计
[n p] = size(X);
f = [zeros(2*p,1); ones(2*n,1)];
%共有n个等式约束条件
b = Y ;%因变量
A = [X,-X,eye(n),-eye(n)];
lb = zeros(2*(n+p),1);
options = optimset('Display','off','MaxIter',100000,'LargeScale', 'off', 'Simplex', 'on');
[xs,fval] =  linprog(f,[],[],A,b,lb,[],[],options);
beta = xs(1:p)-xs(p+1:2*p);
L1norm = sum(abs(Y-X*beta));