3.3.2. Аппроксимация данных с другими нормами
Аппроксимация данных методом наименьших квадратов основана на компромиссе: модельная кривая не может пройти через все точки данных, и тогда 2-норма используется в качестве посредника при выборе параметров. До того как был разработан метод наименьших квадратов, использовались разные способы аппроксимации. Многие их них соответствуют использованию различных норм в качестве посредника.
Например:
- (1 норма)
- (
-норма).
Оптимальные значения параметров для этих норм будут отличаться от значений для 2-ногрмы. Кроме того, вычисления с этими нормами более сложны. Тем не менее они имеют полезные свойства. Например, 1-норма менее чувствительна к присутствию «выбросов» (ложных данных). По этой причине их продолжают использовать.
3.3.3. Аппроксимация данных многочленом заданной степени
Пусть, например, функция задана в виде
набора значений
.
Задача состоит в аппроксимации неизвестной
функциональной зависимости между
и
многочленом заданной степени
:
Для решения этой задачи воспользуемся
методом наименьших квадратов. Согласно
этому методу коэффициенты многочлена
нужно выбрать такими, чтобы сумма
квадратов отклонений найденного
многочлена от заданных значений функции
была минимальной. Другими словами,
коэффициенты
должны минимизировать функцию
В точке минимума функции
ее производные
обращаются
в нуль. Дифференцируя
и приравнивая нулю производные, получим
так называемую нормальную систему
метода наименьших квадратов:
,
где
-
степень аргумента,
-
степень аппроксимирующего многочлена,
-
количество заданных точек из таблицы
значений функции.
Это система алгебраических уравнений
относительно неизвестных
.
Можно показать, что определитель этой
системы отличен от нуля, т. е. решение
системы существует и единственно. Однако
на практике описанную методику применяют
только для нахождения многочленов,
степень которых не выше четырех-пяти.
При более высоких степенях нормальная
система становится плохо обусловленной
(определитель равен нулю) и погрешность
определения коэффициентов велика.
В данном случае заданную табличную
функцию требуется аппроксимировать
многочленом
-й
степени
.
В этом случае нормальная система
уравнений относительно коэффициентов
и
имеет
вид
Эту систему можно представить в матричном виде:
Здесь как
обозначены суммы соответствующих
степеней
-координат заданных точек:
,
а
,
,
используются для обозначения сумм
произведений
-координат заданных точек на
соответствующие степени их
-координат:
.
Полученную систему линейных алгебраических уравнений можно решить одним из известных методов, например методом Гаусса.
Поскольку метод наименьших квадратов
является приближенным методом, то
необходимо оценить погрешность, с
которой выполнена аппроксимация функции.
Итак, если
обозначить полученное в результате
аппроксимации приближенное значение
числа
,
то всегда можно указать (определить)
граничную величину погрешности, некоторое
значение
,
для которого выполняется неравенство
называют пределом (границей) погрешности,
или предельной абсолютной погрешностью,
или сокращенно абсолютной погрешностью,
а
предельной относительной погрешностью
или сокращенно относительной
погрешностью числа
Относительная
погрешность часто указывается в
процентах.
Кроме абсолютной и относительной погрешности необходимо также определить среднее квадратическое отклонение (СКО) по формуле
.
CKО характеризует разброс параметров, или ширину интервала рассеяния, в котором находится большинство точек исходной зависимости.
Пример 1.
{$N+,E+,F+}
PROGRAM MetNKv;
{Аппроксимация функции по методу наимньших квадрвтов}
{Полиномом заданной степени}
Uses CRT;
Const
N=21; {Количество заданных точек кривой намагничивания}
{NN=4;} {Степень полинома}
NameInt = 'mnk.inp'; {Имя файла с данными}
NameOut = 'mnk4-pu.out';} {Имя выходного файла
Type
RealType = Extended; {базовый вещественный тип}
IntType = Integer; {базовый целый тип}
Vec = Array [1..N] of RealType;
TypeS = Array [1..28] of RealType;
TypeA = Array [1..29,1..29] of RealType;
TypeB = Array [1..29] of RealType;
TypeBool = Boolean;
{Глобальные переменные}
Var
S : TypeS;
h,SS,OS,AS : RealType;
Y,X : Vec;
B,P : TypeB;
A : TypeA;
Flag : TypeBool;
F1, F2 : Text;
i, j, k, NN : IntType;
Procedure KoeffA(NN: IntType; Var S: TypeS);
{Вычисление коэффициентов матрицы A}
Var
i, j : IntType;
BEGIN
For j:=1 to 2*NN do
begin
S[j]:=0;
For i:=1 to N do
S[j]:= S[j] + exp(j*ln(X[i]));
end;
END; {KoeffA}
Procedure KoeffB(NN,N: IntType; Y,X: Vec;Var B: TypeB);
{Вычисление правой части системы уравнений B}
Var
i, j : IntType;
BEGIN
For j:=1 to NN+1 do
begin
B[j]:=0;
For i:=1 to N do
B[j]:= B[j] + Y[i]*exp((j-1)*ln(X[i]));
end;
END;{KoeffB}
Procedure GaussKlass(N: IntType; A: TypeA; B: TypeB; Var X: TypeB);
{Решение системы уравнений классическим методом Гаусса}
{Прямой ход исключения переменных}
Var
i, j, k : IntType;
BEGIN
For i:= 1 to N-1 do
begin
For j:= i+1 to N do
begin
A[j,i]:= -A[j,i]/A[i,i];
For k:= i+1 to N do
A[j,k]:= A[j,k] + A[j,i]*A[i,k];
B[j]:= B[j] + A[j,i]*B[i];
end;
end;
{Обратный ход: последовательное нахождение корней}
X[N]:= B[N]/A[N,N];
For i:= N downto 1 do
begin
h:= B[i];
For j:= i+1 to N do h:= h - A[i,j]*X[j];
X[i]:= h/A[i,i];
end;
END; {GaussKlass}
{Головная программа}
BEGIN
ClrScr;
WriteLn;
Write('Введите степень полинома, NN = ');
ReadLn(NN);
Assign(F1,NameInt);
Reset(F1);
For i:=1 to N do Read(F1,X[i]);
For i:=1 to N do Read(F1,Y[i]);
Assign(F2,NameOut);
Rewrite(F2);
WriteLn(F2);
{Вычисление коэффициентов матрицы A}
KoeffA(NN,S);
WriteLn(F2,'МАТРИЦА СИСТЕМЫ УРАВНЕНИЙ A:');
{Заполнение коэффициентов матрицы}
For i:=1 to NN+1 do
begin
For j:=1 to NN+1 do A[i,j]:=S[i+(j-2)];
A[1,1]:= N;
end;
For i:=1 to NN+1 do
begin
For j:=1 to NN+1 do Write(F2,A[i,j]);
WriteLn(F2);
end;
WriteLn(F2);
WriteLn(F2,'ВЕКТОР-СТОЛБЕЦ ПРАВОЙ ЧАСТИ B:');
{Вычисление правой части системы уравнений B}
KoeffB(NN,N,Y,X,B);
For j:=1 to NN+1 do WriteLn(F2, B[j]);
WriteLn(F2);
GaussKlass(NN+1,A,B,P);
WriteLn(F2);
WriteLn(F2,'КОЭФФИЦИЕНТЫ ПОЛИНОМА ',NN:2,' СТЕПЕНИ:');
j:=0;
For i:=1 to NN+1 do
begin
WriteLn(F2,'P',j,'=',P[i],' ');
j:= j+1;
end;
WriteLn(F2);
WriteLn(F2);
WriteLn(F2,'ТАБЛИЦА ЗНАЧЕНИЙ');
WriteLn(F2,' N ',' ЗНАЧ.АРГ.',' ИСТ.ЗН.ФУНК.', ' АПП.ЗН.ФУНК.',' АБС.ПОГР.', ' ОТН. ПОГР.');
WriteLn(F2,' ',' X ',' Y ', ' Y* ',' (Y-Y*) ', ' (Z) ');
SS:=0;
For i:=1 to N do
begin
{Значения аппроксимированной функции}
h:=0;
For j:=1 to NN+1 do
h:= h + P[j]*exp((j-1)*ln(X[i]));
{Абсолютная погрешность}
AS:= Abs(Y[i] - h);
{Относительная погрешность}
OS:= Abs(AS/Y[i])*100;
SS:= SS + AS*AS;
WriteLn(F2,i:3,' ',X[i]:8:2,' ',Y[i]:9:3,' ',h:11:4,' ',AS:9:5,' ',OS:9:5);
end;
{Среднеквадратическое отклонение}
SS:=Sqrt(SS/(N-1));
WriteLn(F2);
WriteLn(F2,'СРЕДНЕКВАДРАТИЧЕСКОЕ ОТКЛОНЕНИЕ:',SS);
Close(F1);
Close(F2);
ClrScr;{Очищает экран, модуль Crt}
END.
Пример 2.
const nmax=25;