Метод наименьших квадратов

6.4. Метод наименьших квадратов

Аппроксимация по методу наименьших квадратов (МНК) доставляет минимум функционалу

где (x_i,y_i) - исходные данные,
f(x) - аппроксимирующая функция.

Если число определяемых параметров функции равно числу отсчетов n, аппроксимация вырождается в интерполяцию.

Линейная регрессия:

Расчетные формулы:

a, b - параметры аппроксимирующей прямой;
s - среднеквадратическая ошибка.

Программа 103. Линейная регрессия

10 VAC
20 PRINT C+1;:INPUT "x",A,"y",B
30 C=C+1:D=D+A:E=E+B:F=F+A*A:G=G+B*B:H=H+A*B:GOTO 20
40A=(D*E-C*H)/(D*D-C*F):B=(E-A*D)/C
50 PRINT "a";A,"b";B,"s";SQR ((G-B*E-A*H)/C):GOTO 20

Размер: 149, H

Пример:

Завершение ввода данных: [AC] RUN 40[EXE].

Программа 104. Регрессия, сводимая к линейной. Трансцендентная двухпараметрическая зависимость y=f(a,b,x) путем замены переменных сводится к линейной

Рассмотрим частный случай y=px/(q+x).

Линеаризация: 1/y=(q/p)(1/x)+1/p;

откуда p=1/b, q=a/b, y=1/y^*, x=1/x^*.

10 VAC
20 PRINT C+1;:INPUT "x",A,"y",B:A=1/A:B=1/B
30 C=C+1:D=D+A:E=E+B:F=F+A*A:G=G+B*B:H=H+A*B:GOTO 20
40A=(D*E-C*H)/(D*D-C*F):B=(E-A*D)/C
50 PRINT "p";1/B,"q";A/B,"s";SQR ((G-B*E-A*H)/C):GOTO 20

Размер: 165, H

Пример:

Завершение ввода данных: [AC] RUN 40[EXE].

Вычисленная оценка погрешности s относится к линеаризованной зависимости.

Программа 105. Нелинейная регрессия. Представленный здесь упрощенный алгоритм МНК принимает все отсчеты равноточными, кроме того, не оценивает погрешность.

В качестве примера рассмотрим зависимость:

Эквивалентное уравнение, линейное по параметрам:

в общем виде: y=a₁y+a₂x²+a₃x+a₄;
откуда: b=a₁, d=a₄, c=-a₂/b, a=a₃-c-db.

Уравнения для частных производных функционала МНК по неизвестным параметрам

образуют систему:

Формулы вычисления элементов матриц:

Для решения системы использован компактный алгоритм метода Гаусса для симметричных матриц.

Размещение в памяти:

В программной области P0:

10 VAC
20 INPUT "x",C,"y",D:A=C*C:B=D*D:E=E+C*B:F=F+A*D:G=G+C*D:H=H+D
30 I=I+A*B:M=M+A*C*D:N=N+A*A:R=R+A*C:S=S+A:W=W+C:X=X+1:GOTO 20

В программной области P1:

10 Q=F:U=G:V=S:A=4
20 FOR D=2 TO A:D(A+D)=E(A*D)/E(A)
30 FOR B=2 TO A:IF D>B;D(A*B+D)=D(A*D+B)/D(A*B+B):GOTO 50
40 FOR C=1 TO D-1:D(A*B+D)=D(A*B+D)-D(A*B+C)*D(A*C+D):NEXT C
50 NEXT B:NEXT D:E=E/E(A):FOR B=2 TO A:FOR C=1 TO B-1
60 D(B)=D(B)-D(A*B+C)*D(C):NEXT C:D(B)=D(B)/D(A*B+B):NEXT B
70 FOR B=A-1 TO 1 STEP -1:FOR C=B+1 TO A
80 D(B)=D(B)-D(A*B+C)*D(C):NEXT C:NEXT B
90 J=E:L=H:K=-F/J:I=G-K+J*L
100 FOR B=1 TO A:PRINT CHR (64+B);H(B):NEXT B

Пример:

Ввод данных: [P0] .5[EXE] 14[EXE] .8[EXE] 17.2[EXE] 1[EXE] 20[EXE] 1.5[EXE] 34[EXE] [AC] [P1]

6.5. Гармоническая аппроксимация

Программа 106. Вычисление коэффициентов Фурье

Отсчеты y_i аппроксимируемой функции делят ее период на m равных интервалов. Порядок аппроксимации n=n_a+n_b.

10 VAC:MODE 5:INPUT "n",A,"m",B
20 FOR D=0 TO B-1:PRINT D;:INPUT E:FOR C=1 TO A
30 E(C)=E(C)+E*SIN (2*p*D/B*INT (C/2)+p*FRAC (C/2)):NEXT C
40 NEXT D:FOR C=1 TO A
50 PRINT CHR (66-2*FRAC (C/2));INT (C/2);2/B*E(C):NEXT C

Размер: 138, E(n)

Пример: n=2; m=3; y₀=0; y₁=0.8660254; y₂=-0.8660254.

Программа 107. Гармоническая интерполяция. В отличие от предыдущей программы предполагается n=m. Интерполяционная формула:

где a_k, b_k - коэффициенты Фурье (см. описание программы 106);

10 VAC:MODE 5:INPUT "m",B
20 FOR D=0 TO B-1:PRINT D;:INPUT E:FOR C=1 TO B
30 F(C)=F(C)+E/B*SIN (2*p*D/B*INT (C/2)+p*FRAC (C/2))
40 NEXT C:NEXT D:INPUT "x0",A,"T",E
50 INPUT "x",D:D=(D-A)/E*2*p:F=G:FOR C=2 TO B
60 F=F+2*F(C)*SIN (D*INT (C/2)+p*FRAC (C/2))
70 NEXT C:PRINT F:GOTO 50

Размер: 184, F(m)

Пример: m=3; y₀=0; y₁=0.8660254; y₂=-y₁; x₀=2p; T=2p; x=p/4.

Программа 108. Интерполяция спектра кусочно-линейной функции с разрывами. Значения действительной (a_f) и мнимой (b_f) составляющих спектра на указанной частоте f определяются по формулам:

где x=pfDt, f - частота;

Dt =t_i-t_i-1 - длина линейного участка аппроксимируемой функции;
y_i^-=y(t_i-0) --значение функции слева от разрыва или излома;
y⁺_i=y(t_i +0) --значение функции справа от разрыва или излома;
n -число интервалов разбиения периода функции.

10 INPUT "n",A:FOR B=1 TO A:PRINT B;
20 INPUT "t",J(3*B),"y-",I(3*B),"y+",K(3*B):NEXT B:K=K(3*A)
30 MODE 5:I=0:J=0:INPUT "f",C:C=2*p*C:FOR B=3 TO 3*A STEP 3
40 D=C*G(B):E=C*J(B):F=(I(B)-H(B))/(E-D):G=SIN D:D=COS D
50 H=SIN E:E=COS E:I=I+F*(E-D)+I(B)*H-H(B)*G
60 J=J+F*(H-G)-I(B)*E+H(B)*D:NEXT B:C=C/2*G(B)
70 PRINT "a";I/C,"b";J/C:GOTO 30

Размер: 262, K(3*n)

При вводе данных y⁺₀ запрашивается в конце, как y_n^-.

Пример: n=2; t₁=1 10^-6; y₁^- =1; y⁺₁=0; t₂=4 10^-6; y₂^- =0; y⁺₂=1.

6.6. Преобразование Фурье

Прямое:

Обратное:

Здесь z_i=u_i=x_i+jy_i - комплексные отсчеты функции времени;
z_k=u_k=x_k+jy_k - комплексные отсчеты функции частоты.

Преобразование справедливо для периодических, либо абсолютно интегрируемых функций. В случае периодических функций отсчеты должны охватывать ровно период, в случае апериодических - всю область ненулевых значений.

Программа 109. Программа допускает разбиение периода функции только на n=2^m интервалов, где m - натуральное.

Ограничение размерности: m≤5.

10 INPUT "m",A,"П",B:C=2-A:DEFM 2*C-12:FOR D=1 TO C:PRINT D;
20 INPUT "x",M(D+D),"y",N(D+D):NEXT D:MODE 5
30 FOR D=1 TO A:E=2-(A-D):G=p/E:F=COS G:G=B*SIN G:H=1:I=0
40 FOR J=1 TO E:FOR K=J+J TO C+C STEP 4*E:L=K+E+E:M=M(K)-M(L)
50 N=N(K)-N(L):M(K)=M(K)+M(L):N(K)=N(K)+N(L)
60 M(L)=M*H-N*I:N(L)=N*H+M*I:NEXT K:K=H*F+I*G:I=I*F-H*G
70 H=K:NEXT J:NEXT D:FOR D=1 TO C-1:IF D≥E THEN 100
80 FOR H=0 TO 1:M(H)=M(E+E+H):M(E+E+H)=M(D+D+H)
90 M(D+D+H)=M(H):NEXT H
100 I=C/2
110 IF I<E;E=E-I:I=I/2:GOTO 110
120 E=E+I:NEXT D:IF B>0;D=1
130 FOR E=1 TO C:PRINT E,"x";M(E+E)/D,"y";N(E+E)/D:NEXT E

Размер: 438, N(2*2-m)

Для прямого преобразования необходимо задать "П"=1, для обратного - "П"=-1.

Пример: m=3; "П"=1;

Проверка по свойству обратимости. Применение обратного преобразования к результатам прямого даст исходную функцию с двойной вычислительной погрешностью. Достаточно подать команды B=-B и RUN 30.

Ниже представлена модификация программы 109, рассчитанная на число отсчетов n=64. Ввод и вывод данных ведется в режиме калькулятора.

10 FOR D=1 TO 6:E=64-2-D:G=p/E:F=COS G:G=SIN G:H=1:I=0
20 FOR J=1 TO E:FOR K=J+J TO 128 STEP 4*E:A=K+E+E:B=J(K)-J(A)
30 C=K(K)-K(A):J(K)=J(K)+J(A):K(K)=K(K)+K(A)
40 J(A)=B*H-C*I:K(A)=C*H+B*I:NEXT K:K=H*F+I*G:I=I*F-H*G
50 H=K:NEXT J:NEXT D:FOR D=1 TO 63:IF D≥E THEN 80
60 FOR H=0 TO 1:J(H)=J(E+E+H):J(E+E+H)=J(D+D+H)
70 J(D+D+H)=J(H):NEXT H
80 I=C/2
90 IF I<E;E=E-I:I=I/2:GOTO 90
100 E=E+I:NEXT D

Размер: 1212

Ввод/вывод данных

Прямое

Обратное