防灾数学分享 http://blog.sciencenet.cn/u/fzmath 防灾科技学院数学教研室

博文

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

已有 4675 次阅读 2010-6-15 11:33 |个人分类:生活点滴|系统分类:科研笔记| 线性模型, T-型稳健估计

以下存为一个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));




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

上一篇:今天编出了线性模型和线性EV模型T-型估计的MATLAB程序
下一篇:高斯消去法的matlab程序

0

发表评论 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-10-2 06:38

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部