Метки

, , , , ,

В эксперименте были намеряны значения \{y_i\} = \{y_1, y_2, \ldots, y_n\}, которым соответствуют \{x_i\} = \{x_1, x_2, \ldots, x_n \}. После получения апроксимирующего полинома (в нашем случае прямой Y = ax+b), нужно узнать погрешности коэффициентов и самой аппроксимации.

Остаточная дисперсия (она же дисперсия адекватности) имеет вид
S^2_{add}(y) = \frac{\sum (Y_i - y_i)^2}{n - k},
где k — количество коэффициентов полинома. Считается, что модель выбрана хорошая, если остаточная дисперсия сравнима по своему значению с дисперсией воспроизводимости, равной
S^2_{disp}(y) = \frac{\sum (\overline y_i - y_i)^2}{n - k - 1}.

В общем случае для f(x_1, x_2, \ldots, x_n), S^2(f) = \sum\limits_i \left( S^2(x_i) \left( \frac{\partial f}{\partial x_i} \right)^2\right).

Далее для линейной зависимости
S^2(a) = \frac{n}{n\sum\limits_i x_i^2 - \left(\sum\limits_i x_i \right)^2} S^2(y);
S^2(b) = \frac{\sum\limits_i x_i^2}{n\sum\limits_i x_i^2 - \left(\sum\limits_i x_i \right)^2} S^2(y).

Соответствующий матлабовский скрипт:
function PP = LSM(x, y, k)

% — Усреднение по методу наименьших квадратов —-
%
% Входные параметры:
% x, y — векторы абсциссы и ординаты, где y = f(x);
% k — количество параметров МНК, фактически равно степени полинома,
% которым приближается зависимость, плюс 1 (для линейной апроксимации
% k = 2).
%
% Выход функции:
% P — массив коефициентов МНК. Для линейной зависимости
% y = a * x + b = P(1) * x + P(2). Остальные параметры несут смысл лишь
% при линейной зависимости (k = 2).
% S_a, S_b — дисперсии a и b.
% S_X и S_Y — дисперсии X и Y (обратная и прямая регресионная задача).

n = length(x);
if (n — length(y) ~= 0)
disp(‘Векторы должны иметь одинаковую размерность!’);
return;
end

if (n == k)
disp(‘Размерность апроксимирующего полинома должна быть меньше’);
disp(‘количества длины входного вектора данных!’);
return;
end

P = polyfit(x, y, k — 1);
Y = polyval(P, x);

S_y2 = 0;
temp2 = 0;
temp1 = 0;
for i=1:1:n
S_y2 = S_y2 + (Y(i) — y(i)).^2;
temp1 = temp1 + x(i);
temp2 = temp2 + x(i).^2;
end

S_y2 = S_y2 / (n — k);
if (k == 2)
temp3 = n * temp2 — temp1^2;
S_a = sqrt(n * S_y2 / temp3);
S_b = sqrt(temp2 / temp3 * S_y2);
S_Y = sqrt(S_y2);
S_X = S_Y * sqrt(1/n + 1 + n * (n — 1) * std(y)^2 / temp3 / P(1)^2) / P(1);
else
S_a = -1;
S_b = -1;
S_Y = -1;
S_X = -1;
end

disp(‘Коэффициенты аппроксимирующего полинома:’);
disp(P);
if (k == 2)
disp(‘Линейная аппроксимация Y = ax + b:’);
disp(‘Дисперсия a S_a = ‘);
disp(S_a);
disp(‘Дисперсия b S_b = ‘);
disp(S_b);
disp(‘Обратная регрессионная задача S_X = ‘);
disp(S_X);
disp(‘Прямая регрессионная задача S_Y = ‘);
disp(S_Y);
end
PP = [P S_a S_b S_X S_Y];
return;

Реклама