Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/~rastor/Software/linsvd.m Дата изменения: Sun Mar 3 23:38:04 2013 Дата индексирования: Fri Feb 28 00:30:15 2014 Кодировка: Windows-1251 Поисковые слова: molecular cloud

function [solution,sigma,residual]=linsvd(matrix,y,plotarg);

%% Программа вычисляет коэффициенты линейной модели МНК
%% (метода наименьших квадратов) вида
%
% y(i) = a(1)*x(1) + ... +a(M)*x(M), i=1,2,...,N > M,
%
% которая может быть записана в матричном виде A * X = Y,
% где вектор-столбец X - неизвестные [x(1) x(2) ... x(M)]',
% а Y - вектор-столбец "левых частей" [y(1) y(2) ... y(N)]',
% матрица A (M столбцов, N > M строк) - матрица известных
% коэффициентов a(i,j) при неизвестных x(j), т.е.
% |a(1,1) a(1,2) ... a(1,M) |
% |a(2,1) a(2,2) ... a(2,M) |
% A = |..............................|
% |a(N-1,1) a(N-1,2) ... a(N-1,M)|
% |a(N,1) a(N,2) ... a(N,M) |
%
% K - число неизвестных, N > M - число условных уравнений
%
% Примечание: по умолчанию данные считаются равноточными,
% но для неравноточных наблюдений/измерений следует
% перед вызовом программы поделить все коэффициенты a(i,j)
% i-й строки и левые части y(i) на ошибки левых частей sigma_y(i)
% В такой модификации программа реализует chi^2-минимизацию
%
% Обращение к программе:
%
% [solution,sigma[,residual]]=linsvd(matrix,y[,plotarg]);
%
% (в скобках указаны необязательные входные (выходные) параметры
%% Вход:
% matrix - матрица A, столбцы которой являются коэффициентами
% при неизвестных для системы условных уравнений
% y - вектор-столбец "левой части" y
% [plotarg] - если задан, выводится график зависимости от него
% остаточных уклонений (residual)(необязательный
% параметр)
%% Выход:
% solution(:,1) - SVD МНК-решение линейной системы
% solution(:,2) - ошибки линейных коэффициентов
% sigma - среднеквадратичная невязка (в смысле chi^2,
% т.е. средневзвешенная)
% [residual] - массив уклонений от решения (необязательный параметр)
%
%% Автор: А.С.Расторгуев, физфак МГУ, ГАИШ МГУ, 2000-2013
%%

[n,m]=size(matrix);

%% SVD-разложение исходной матрицы:
[U,S,V]=svd(matrix,0);
for i=1:m
% Сингулярные числа:
w(i)=S(i,i);
S(i,i)=1/S(i,i);
end
solution=zeros(m,2);
% Решение SVD:
solution(:,1)=V*S*U'*y;
zsol=matrix*solution(:,1);
%% Вычислим сумму квадратов невязок chi^2:
chi2=0;
residual(:,1)=y(:)-zsol(:);
for i=1:n
chi2=chi2+residual(i)^2;
end
chi2=sqrt(chi2/(n-m));
sigma=chi2;
%% Ошибки неизвестных:
for j=1:m
solution(j,2)=0;
for i=1:m
solution(j,2)=solution(j,2)+(V(j,i)*S(i,i))^2;
end
solution(j,2)=sqrt(solution(j,2))*chi2;
end
%
%% График остаточных уклонений
if nargin>2,
% subplot(211)
% plot(plotarg,zsol,'g.',plotarg,y,'r.');grid;xlabel('Plot argument');ylabel('Dependent var');
% subplot(212)
plot(plotarg,residual,'b.');grid;xlabel('Plot argument');ylabel('Residuals');
end;