Аппроксимация. Метод наименьших квадратов

Пусть данные некоторого эксперимента представлены в виде таблицы значений независимой переменной x и зависимой переменной y.

Требуется отыскать аналитическую зависимость f(x,a₀,a₁,…,a_m), являющуюся функцией одной независимой переменной x и параметров a₀,a₁,a₂,…,a_m, которая наилучшим образом описывала бы эти экспериментальные данные в смысле минимума квадратичного критерия рассогласования R(a₀,a₁,…,a_m):

Функцию f(x,a₀,a₁,…,a_m) определим как полином степени m вида: Надо найти такие значения параметров, при которых квадратичный критерий рассогласования имел бы минимальное значение  Вывод формулы для определения параметров в матричном виде рассмотрим на примере полинома второй степени (m=2).  Тогда критерий R будет являться функцией трёх переменных a₀, a₁, a₂:  Необходимые условия минимума критерия R имеют вид:

или

Полученную линейную относительно искомых параметров a0,a1,a2, систему уравнений запишем в матричном виде:  где

 Для удобства формирования матрицы коэффициентов  и столбца свободных членов  введем матрицу  элементы которой определяются через значения независимой переменной x_i, i=0,1,2,…,n  тогда

При аппроксимации полиномами высших порядков матрица  будет иметь вид:

В общем случае количество строк в матрице  равно количеству точек, а количество столбцов равно количеству параметров, где строка состоит из значений частных производных от функции f(x,a₀,a₁,…,a_m) по соответствующему параметру.

 

Программа. Аппроксимация.

 

function DATA

global x y m n;

x=[-4:4];

y=[254;75;12;0;1;2;21;92;287];

n=9;

m=4;

end

 

function [ a,dy,r,yr ] = fun_Approximation(x,y,m,n)

for i=1:n

for j=1:m+1

   F(i,j)=x(i)^(j-1);

end %for j

end %for i

Ft=F';

N=Ft*F;

b=Ft*y;

Nobr=N^(-1);

a=Nobr*b;

yr=F*a;

dy=y-yr;

r=dy'*dy;

 end % function

 

function GLAV_Approximation

global x y m n a dy r yr;

DATA;

[ a,dy,r,yr ] = fun_Approximation(x,y,m,n);

REPORT;

end

 

function REPORT

global x y m n a dy r yr;

disp('Approximation');

disp(' ');

disp('     i    x        y      yr    dy'); 

disp(' ________________________________________________________')

i=0;

for i=1:length(x)

xx=x(i);

yy=y(i);

yrr=yr(i);

dyy=dy(i);

disp(sprintf('%10.1f\t%10.5f\t %10.5f\t%10.5f\t %10.5f',i,xx,yy,yrr,dyy));

end %for

disp(' ');

disp('   i    a '); 

disp('  _______________')

i=0;

for i=1:length(a)

aa=a(i);

disp(sprintf('%10.1f\t%10.5f',i,aa));

end %for

plot(x,y,'r*',x,yr,'k-');

grid on;

xlabel('x *, xx --');

ylabel('y *, yr --');

title('Approximation');

end