# 零截距线性模型的T-型稳健估计的matlab代码

function [beta sigma fvalue] = t_estimator40lm(X,Y,v)
% 普通线性模型T-型估计的EM算法
%%%初始化%%%%%%
[n,p] = size(X);
%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);

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

http://blog.sciencenet.cn/blog-292361-335612.html

## 全部精选博文导读

GMT+8, 2020-10-29 18:49