|||
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));
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2025-1-10 03:33
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社