数据正则化matlab程序,求助 计量经济模型中用Tikhonov正则化方法参数估计程序修改...

我最近正在做毕业论文，是关于不适定反问题的，最后要用多种正则化方法(如Tikhonov正则化方法)用MATLAB编程数值实现得到结果，问题是在一计量经济模型中有三个因自变量x1 x2 x3对结果y有贡献，要通过一些历史数据估计每个变量对结果的贡献系数，即估计y=a1*x1+a2*x2+a3*x3+u中的参数a1 a2 a3，u为随机干扰项。

现在只有一些算法，不会自己修改，求高手不胜感激啊！

function [x_lambda,rho,eta] = tikhonov(U,s,V,b,lambda,x_0)

%TIKHONOV Tikhonov regularization.

%

% [x_lambda,rho,eta] = tikhonov(U,s,V,b,lambda,x_0)

% [x_lambda,rho,eta] = tikhonov(U,sm,X,b,lambda,x_0) , sm = [sigma,mu]

%

% Computes the Tikhonov regularized solution x_lambda. If the SVD

% is used, i.e. if U, s, and V are specified, then standard-form

% regularization is applied:

% min { || A x - b ||^2 + lambda^2 || x - x_0 ||^2 } .

% If, on the other hand, the GSVD is used, i.e. if U, sm, and X are

% specified, then general-form regularization is applied:

% min { || A x - b ||^2 + lambda^2 || L (x - x_0) ||^2 } .

%

% If x_0 is not specified, then x_0 = 0 is used

%

% Note that x_0 cannot be used if A is underdetermined and L ~= I.

%

% If lambda is a vector, then x_lambda is a matrix such that

% x_lambda = [ x_lambda(1), x_lambda(2), ... ] .

%

% The solution norm (standard-form case) or seminorm (general-form

% case) and the residual norm are returned in eta and rho.

% Per Christian Hansen, IMM, April 14, 2003.

% Reference: A. N. Tikhonov & V. Y. Arsenin, "Solutions of

% Ill-Posed Problems", Wiley, 1977.

% Initialization.

if (min(lambda)<0)

error('Illegal regularization parameter lambda')

end

m = size(U,1);

n = size(V,1);

[p,ps] = size(s);

beta = U(:,1:p)'*b;

zeta = s(:,1).*beta;

ll = length(lambda); x_lambda = zeros(n,ll);

rho = zeros(ll,1); eta = zeros(ll,1);

% Treat each lambda separately.

if (ps==1)

% The standard-form case.

if (nargin==6), omega = V'*x_0; end

for i=1:ll

if (nargin==5)

x_lambda(:,i) = V(:,1:p)*(zeta./(s.^2 + lambda(i)^2));

rho(i) = lambda(i)^2*norm(beta./(s.^2 + lambda(i)^2));

else

x_lambda(:,i) = V(:,1:p)*...

((zeta + lambda(i)^2*omega)./(s.^2 + lambda(i)^2));

rho(i) = lambda(i)^2*norm((beta - s.*omega)./(s.^2 + lambda(i)^2));

end

eta(i) = norm(x_lambda(:,i));

end

if (nargout > 1 & size(U,1) > p)

rho = sqrt(rho.^2 + norm(b - U(:,1:n)*[beta;U(:,p+1:n)'*b])^2);

end

elseif (m>=n)

% The overdetermined or square general-form case.

gamma2 = (s(:,1)./s(:,2)).^2;

if (nargin==6), omega = V\x_0; omega = omega(1:p); end

if (p==n)

x0 = zeros(n,1);

else

x0 = V(:,p+1:n)*U(:,p+1:n)'*b;

end

for i=1:ll

if (nargin==5)

xi = zeta./(s(:,1).^2 + lambda(i)^2*s(:,2).^2);

x_lambda(:,i) = V(:,1:p)*xi + x0;

rho(i) = lambda(i)^2*norm(beta./(gamma2 + lambda(i)^2));

else

xi = (zeta + lambda(i)^2*(s(:,2).^2).*omega)./...

(s(:,1).^2 + lambda(i)^2*s(:,2).^2);

x_lambda(:,i) = V(:,1:p)*xi + x0;

rho(i) = lambda(i)^2*norm((beta - s(:,1).*omega)./...

(gamma2 + lambda(i)^2));

end

eta(i) = norm(s(:,2).*xi);

end

if (nargout > 1 & size(U,1) > p)

rho = sqrt(rho.^2 + norm(b - U(:,1:n)*[beta;U(:,p+1:n)'*b])^2);

end

else

% The underdetermined general-form case.

gamma2 = (s(:,1)./s(:,2)).^2;

if (nargin==6), error('x_0 not allowed'), end

if (p==m)

x0 = zeros(n,1);

else

x0 = V(:,p+1:m)*U(:,p+1:m)'*b;

end

for i=1:ll

xi = zeta./(s(:,1).^2 + lambda(i)^2*s(:,2).^2);

x_lambda(:,i) = V(:,1:p)*xi + x0;

rho(i) = lambda(i)^2*norm(beta./(gamma2 + lambda(i)^2));

eta(i) = norm(s(:,2).*xi);

end

end

相关数据见下：

表二 1978-2008我国宏观经济数据增长率

年度 国内生产总值增长率 最终消费增长率 资本形成总额增长率 货物和服务净出口增长率

1978 - - - -

1979 11% 17% 7% 72%

1980 12% 14% 8% -24%

1981 8% 11% -1% -176%

1982 9% 10% 11% 706%

1983 12% 11% 14% -44%

1984 21% 17% 23% -97%

1985 25% 23% 37% -28323%

1986 14% 13% 14% -30%

1987 17% 14% 12% -105%

1988 25% 26% 27% -1414%

1989 13% 13% 11% 23%

1990 10% 8% 6% -375%

1991 17% 16% 17% 21%

1992 23% 21% 28% -55%

1993 30% 27% 56% -347%

1994 35% 33% 28% -193%

1995 25% 26% 24% 57%

1996 18% 19% 13% 46%

1997 9% 9% 6% 96%

1998 6% 6% 4% 7%

1999 3% 7% 4% -26%

2000 9% 10% 6% 6%

2001 7% 8% 22% -3%

2002 9% 7% 15% 33%

2003 12% 7% 23% -3%

2004 17% 29% 24% 37%

2005 33% 11% 17% 151%

2006 16% 14% 17% 63%

2007 21% 17% 18% 40%

2008 17% 16% 20% 3%

2010-5-18 16:38 上传

点击文件名下载附件

1.16 KB, 下载次数: 438