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

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

Общий и линейный методы

2. Общий метод наименьших квадратов

Нередко при анализе экспериментальных данных требуется их сопоставить с
результатами физической модели, описывающей изучаемое явление. Чаще всего
речь идет о сравнении экспериментально полученной функциональной
зависимости с соответствующей аналитической кривой, полученной из теории.
При этом возникают, по крайней мере, две задачи: имеется ли соответствие между
кривой и экспериментально полученной зависимостью и как определить
неизвестные параметры, входящие в выражение для теоретической кривой.
Рассмотрим решение этих задач.
Проблему можно сформулировать следующим образом. Пусть в результате
измерения некоторой физической величины y, зависящей от x, для ряда значений
xi, известных точно, получены независимо друг от друга значения yi (i = 1, 2, ...., n)
с погрешностью σi (среднеквадратичные отклонения). Предположим, что эта
экспериментальная зависимость может быть описана (аппроксимирована)
выражением y = f(x, 1, 2, ..., k), где 1,..., k - неизвестные параметры, а
аналитический вид функции f (x, 1, 2, ..., k), выбирается на основе физической
модели какого либо явления или процесса. Как правило n ≥ k.

3. Общий метод наименьших квадратов

Аппроксимация экспериментальных данных.

4. Общий метод наименьших квадратов

Это означает, что существуют такие значения параметров 01,..., 0k, что
y0i = f(xi, 01, ..., 0k) являются математическими ожиданиями случайных величин
yi (i = 1, ..., n), то есть E(yi) = y0i = f(xi, 01, ..., 0k), где E - знак операции
вычисления математического ожидания случайной величины. Требуется найти
наилучшие оценки параметров 1, ..., k, используя данные yi ± σi. Важным
шагом в решении этой задачи является выбор меры близости между
теоретической и экспериментальной кривыми. В регрессионном анализе
постулируется выбор среднеквадратичной меры близости, согласно которой
расстояние между набором теоретически вычисленных значений f(xi, 1, …, k) и
набором экспериментальных значений определяется формулой:
где Wi - статистический вес i-го измерения. Обычно веса выбираются аналогично
случаю неравноточных измерений (для получения эффективной оценки), то есть
Wi =1/ i2, где i - среднеквадратичная погрешность измерений yi, которая
считается известной. Для отыскания значений параметров j выбираются те
значения, при которых сумма взвешенных квадратов отклонений теоретических
значений от экспериментальных значений в точках xi минимальна.

5. Общий метод наименьших квадратов

Таким образом, оценки по методу наименьших квадратов (МНК - оценки)
параметров j могут быть найдены как решение системы уравнений
Возможным способом решения этой задачи является следующий.
Выбираем число экспериментальных значений n, равным числу неизвестных
параметров k. Подставляем пары чисел (xi,yi) в систему уравнений
y1 = f(x1, 1, 2, ..., k)
..................
yk = f(xn, 1, 2, ..., k)
Получим k уравнений с k неизвестными параметрами j. Подразумевается, что
решение такой системы существует. Однако существенным недостатком такого
способа является то обстоятельство, что решение такой системы не всегда
существует, а также никак не учитывается возможное различие в погрешностях
отдельных измерений.

6. Общий метод наименьших квадратов

Поэтому выбирается другой способ решения задачи. Пусть число
экспериментальных точек n>k. В этом случае получаем избыточную систему
уравнений, не имеющую совместного решения. Следовательно, существует
единственная возможность – найти такие значения параметров ( 01,..., 0k), при
которых значения теоретической кривой (4.1) в точках xi с учетом погрешности
измерений yi максимально близко подходили бы к экспериментальным значениям.
Таким образом, необходимо найти такие значения параметров ( 01,..., 0k), при
которых функция
достигает своего минимального значения, а значения параметров 0j могут быть
найдены как решение системы уравнений

7. Проверка качества аппроксимации

О качестве аппроксимации экспериментальных данных с помощью выбранной
зависимости f (x, 1, ...., k) можно судить по остаточной сумме квадратов
Sост = S( 01, ...., 0k).
.
Известно, что если yi распределены по нормальному закону, а аппроксимирующая
функция такова, что E(yi) = f(xi, 01, ...., 0k), причем дисперсии i2 значений yi
известны точно, то в этом случае величина
или
подчиняется распределению 2 с n - k степенями свободы ( n переменных и k
уравнений связи для нахождения параметров 0k ).
Таким образом, с помощью критерия согласия 2 выполняется проверка гипотезы
о корректности выбора аппроксимирующей функции y = f (x, 1, ...., k).

8. Линейный метод наименьших квадратов

Пусть функция f (x, 1,...., k) зависит от параметров j линейно, то есть
Этот случай является одним из самых распространенных на практике. Часто
требуется представить экспериментальные данные в виде разложения в ряд по
каким-либо функциям (степенным функциям, тригонометрическим функциям и
т.д.). Тогда систему уравнений
y1 = f (x1, 1, 2, ..., k)
..................
yn = f (xn, 1, 2, ..., k)
можно записать в виде:
y1 = 1φ1(x1) + 2φ2(x1) +…+ kφk(x1)
...............................
yn = 1φ1(xn) + 2φ2(xn) +…+ kφk(xn)
Для решения этой задачи вводится матрица A, с элементами
Aij = j(xi) (i = 1, 2, ..., n; j = 1, 2, ..., k). Она называется конструкционной матрицей.

9. Линейный метод наименьших квадратов

Введем вектор-столбец с компонентами
и вектор-столбец Y с
компонентами Y =
(по аналогии с E(yi) = y0i = f (xi, 01, ..., 0k),
обозначим Y0 =
).
Чтобы решить эту систему уравнений, то есть найти вектор , необходимо
решить матричное уравнение
. Однако для корректного решения уравнения
необходимо учесть точность измерений величин yi , то есть вклад каждого
уравнения в решение необходимо учесть со своим статистическим весом. Для
этого обозначим через W диагональную матрицу с отличными от нуля элементами
Wii :
Умножим слева обе части уравнения
на матрицу W :
. Чтобы найти решение для вектора параметров
нужно умножить слева обе части этого уравнения на обратную матрицу
.
Однако матрица
не является квадратной и для нее обратной матрицы не
существует. Поэтому сначала умножим это уравнение слева на матрицу
:
. Теперь матрица
стала квадратной и симметричной
. Тогда решение для вектора параметров примет вид:
.
Для дисперсии отдельных параметров получается выражение
D( j)= [(AT W A)-1]jj .

10. Линейный метод наименьших квадратов

Поскольку значения случайных величин yi получены независимо друг от друга с
дисперсиями i2, ковариационная матрица D (или матрица ошибок) вектора Y
может быть представлена как диагональная матрица:
Вспоминая, что Wi =1/ i2, имеем
, где
матрица, обратная к W.

11. Свойства МНК – оценок

Оценки параметров 1, 2, ..., k, получаемые с помощью МНК, - случайные
величины, поскольку они основаны на результатах измерений случайных величин
yi и являются линейными комбинациями последних. По определению
ковариационная матрица случайного вектора :
D( ) = E[ ( - E( ))·( - E( ))]
Используя правило транспонирования матриц (ABC)T = CT BT AT, получаем:
D( ) = E [(B-1 AT W ( Y – Y0))·(B-1 AT W ( Y – Y0))T ]=
= B-1 AT W E [( Y – Y0 )T·( Y – Y0 )]( B-1 AT W )T , = B-1 AT W D( Y ) ( B-1 AT W )T.
Подставляя вместо D(Y) выражение
, а также учитывая, что
и
, получаем:
D( ) = B-1ATWW-1WTA( B-1)T = B-1AT W A B-1 = B-1BB-1 = B-1.
Таким образом
D( ) = B-1 = (AT W A)-1.
Для дисперсии отдельных параметров получается выражение
D( j)= [(AT W A)-1]jj .

12. Свойства МНК – оценок

Заметим, что дисперсия оценок параметров j определяется непосредственно в
процессе отыскания самих значений параметров, при обращении матрицы B.
Следует отметить, что большие значения погрешностей оценок параметров j
относительно самих значений j, отнюдь не свидетельствуют о плохом качестве
метода оценивания, а говорят лишь о том, что в достаточно широких пределах
изменения параметра остаточная сумма Sост изменяется незначительно. Иными
словами, наблюдается слабая зависимость Sост от данного параметра. Это может
говорить о том, что в теоретической модели этот параметр избыточен.
Иногда погрешности измерений 1, ..., n заранее не известны и подлежат оценке
по результатам измерений y1, ..., yn. Тем не менее, как правило, удается
определить отношения дисперсий ошибок i2/ i+12. В этом случае целесообразно
ввести статистические веса следующим образом i = 02/ i2, где величину 02,
называемую дисперсией наблюдения с единичным весом, следует определить.
Перейдя к нормированным статистическим весам, нетрудно увидеть, что такое
преобразование никоим образом не изменит вид МНК-оценки параметров ,
который не зависит от 02. Однако ковариационная матрица оценок параметров
D( ) будет содержать неизвестный множитель 02, так как теперь ковариационная
матрица результатов измерений Y будет иметь вид D(Y) = 02 W-1 и может быть
представлена как:
D( ) = 02 B-1 = 02 (AT W A)-1 .

13. Свойства МНК – оценок

Чтобы получить матрицу ошибок D( ), необходимо либо знать точно 02, либо
иметь ее оценку.
Можно показать, что несмещенная оценка 02 в МНК дается выражением
02 = Sост /(n-k) = VT W V /(n-k),
где V = y – A α - вектор остатков, Sост - остаточная сумма квадратов
невязок; n - число экспериментальных точек; k - число неизвестных параметров
апроксиматора.
Окончательно, несмещенная ковариационная матрица несмещенных оценок
параметров имеет вид:
МНК-оценки параметров линии регрессии обладают меньшей дисперсией, чем
любые другие оценки линейного вида, то есть являются эффективными. Это
утверждение известно под названием теоремы Гаусса-Маркова. Кроме того,
МНК-оценки являются состоятельными, а для линейного случая они еще и не
смещены.

14. Пример использования общего метода наименьших квадратов

Аппроксимировать зависимостью y(x) = a·exp(-((x-b)2)/c)+d данные
x = (0 1 2 3 4 5 6 7 8 9)
y = (3 3.6 4.5 6 9 11.1 8.7 6.3 4.2 3.3)
y = (1 1 1 1.1 1.2 1.2 1.2 1 1 1 ). Проверить качество аппроксимации
2
(
x
b
)
с помощью критерия . Считать Pдов = 0.95.
0
3
1
1
3.6
1
2
4.5
1
3
6
1.1
4
X
Y
5
9
11.1
1.2
s
1.2
6
8.7
1.2
7
6.3
1
8
4.2
1
9
3.3
1
2. c
a. e
f( x , a , b , c , d )
d
В аппроксимирующей функции наблюдается нелинейная
зависимость от параметров. Для вычисления параметров
используем общий метод наименьших квадратов.
Произведем выбор начальных приближений параметров
с использованием графики.
12
N
a
8
b
5
c
2
d
3
i
10
10
0 .. N
Y
i
1
f X , a, b, c, d
i
8
6
4
2
0
2
4
6
i
8
10

15. Пример использования общего метода наименьших квадратов

Определим меру близости экспериментальной и теоретической зависимостей.
N
1
S( a , b , c , d )
i= 0
Yi
f Xi , a , b , c , d
si
2
Given
2
S( a , b , c , d ) 0
1
A
7.46
A=
4.988
2.205
1
1
1
- вычисленные параметры
x
0 , 0.1 .. N
4
6
1
12
Проверка качества аппроксимации с
использованием критерия . P 0.95
qchisq( P , N
1
Minerr( a , b , c , d )
3.199
S A0 , A1 , A2 , A3 = 0.54
1
4) = 12.592
10
Y
i
8
f( x, a , b , c , d )
Качество аппроксимации хорошее, так
f x, A 0 , A 1 , A 2 , A 3 6
как вычисленное значение величины
4
меньше критического значения
= 12.592.
Визуальная проверка качества аппроксимации: 2
0
2
i, x, x
8
10

16. Пример использования линейного метода наименьших квадратов

Аппроксимировать зависимостью y (x) = Ax2 + Bx + C данные
x = (0 1 2 3 4 5 6 7 8 9)
y = (3 3.3 4.0 5.1 6.5 8.9 10.7 12.9 15.8 20.2)
y = (0.6 0.6 0.8 1.0 1.1 1.2 1.3 1.4 1.5 1.7). Проверить качество
0 1
аппроксимации с помощью критерия . Считать Pдов = 0.9.
X
0
3
0.6
1
3.3
0.6
2
4.0
0.8
3
5.1
1
4
6.5
1.1
5
Y
8.9
1.2
6
10.7
1.3
7
12.9
1.4
8
15.8
1.5
9
20.2
1.7
W
2
a. x
y( a , b , c , x)
b. x
00
0
1
11
1
1
24
2
1
39
3
1
4 16 4
1
A = 5 25 5
6 36 6
1
7 49 7
1
8 64 8
1
9 81 9
1
c
Построим конструкционную матрицу А
N last( x) K
3
i
j
0.. N 1
0.. K
x0 1 x2 x
i
i
i
1
либо A augment( x2 augment( x x0) )
и весовую матрицу W
Ai , j
w
diag( w)
Xi
1
2
2
j
0
1
0 2.778 0
2
3
4
5
6
7
8
9
0
0
0
0
0
0
0
0
1 0
2.778 0
0
0
0
0
0
0
0
2 0
0
1.562 0
0
0
0
0
0
0
3 0
0
0
1
0
0
0
0
0
0
4 0
0
0
0
0.826 0
0
0
0
0
W=5 0
6 0
0
0
0
0
0.694 0
0
0
0
0
0
0
0
0
0.592 0
0
0
7 0
0
0
0
0
0
0
0.51 0
0
8 0
0
0
0
0
0
0
0
0.444 0
9 0
0
0
0
0
0
0
0
0
0.346
2
2
1

17. Пример использования линейного метода наименьших квадратов

0.191
1
. A T. W . Y
Вычислим вектор искомых параметров α ={a,b,c}:
Проверка качества аппроксимации с использованием критерия
T.
A W. A
S
T
( Y A. ) . W. ( Y A. )
P
S = 0.437
0.9
= 0.136
2.989
qchisq( P, N K) = 12.017
Качество хорошее, так как вычисленное значение величины S = 0.437 меньше
критического значения
= 12.017.
Визуальная проверка качества аппроксимации
z
1
T
D A W A j j
j
25
D
20
0, 0.1.. 10
y , , , z 15
0 1 2
Y
i
0.191
= 0.136
2.989
10
0.047
5
0.508
0
0.371
0
2
4
6
z, i
8
10

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

Задание:
Данные предыдущей задачи аппроксимировать c помощью функции
f (x, a, b, c, d) = ax3+bx2+cx+d. Определить возможный избыточный параметр.