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

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

Краткая инструкция по методу наименьших квадратов для чайников: определяем вид предполагаемой зависимости (чаще всего берется линейная регрессия вида $y(x)=ax+b$), выписываем систему уравнений для нахождения параметров $a, b$. По экспериментальным данным проводим вычисления и подставляем значения в систему, решаем систему любым удобным методом (для размерности 2-3 можно и вручную). Получается искомое уравнение.

Иногда дополнительно к нахождению уравнения регрессии требуется: найти остаточную дисперсию, сделать прогноз значений, найти значение коэффициента корреляции, проверить качество аппроксимации и значимость модели. Примеры решений вы найдете ниже. Удачи в изучении!

Видео:Как работает метод наименьших квадратов? Душкин объяснитСкачать

Примеры решений МНК

Пример 1. Методом наименьших квадратов для данных, представленных в таблице, найти линейную зависимость

Пример 2. Прибыль фирмы за некоторый период деятельности по годам приведена ниже:
Год 1 2 3 4 5
Прибыль 3,9 4,9 3,4 1,4 1,9
1) Составьте линейную зависимость прибыли по годам деятельности фирмы.
2) Определите ожидаемую прибыль для 6-го года деятельности. Сделайте чертеж.

Пример 3. Экспериментальные данные о значениях переменных х и y приведены в таблице:
1 2 4 6 8
3 2 1 0,5 0
В результате их выравнивания получена функция Используя метод наименьших квадратов, аппроксимировать эти данные линейной зависимостью (найти параметры а и b). Выяснить, какая из двух линий лучше (в смысле метода наименьших квадратов) выравнивает экспериментальные данные. Сделать чертеж.

Пример 4. Данные наблюдений над случайной двумерной величиной (Х, Y) представлены в корреляционной таблице. Методом наименьших квадратов найти выборочное уравнение прямой регрессии Y на X.

Пример 5. Считая, что зависимость между переменными x и y имеет вид $y=ax^2+bx+c$, найти оценки параметров a, b и c методом наименьших квадратов по выборке:
x 7 31 61 99 129 178 209
y 13 10 9 10 12 20 26

Пример 6. Проводится анализ взаимосвязи количества населения (X) и количества практикующих врачей (Y) в регионе.
Годы 81 82 83 84 85 86 87 88 89 90
X, млн. чел. 10 10,3 10,4 10,55 10,6 10,7 10,75 10,9 10,9 11
Y, тыс. чел. 12,1 12,6 13 13,8 14,9 16 18 20 21 22
Оцените по МНК коэффициенты линейного уравнения регрессии $y=b_0+b_1x$.
Существенно ли отличаются от нуля найденные коэффициенты?
Проверьте значимость полученного уравнения при $alpha = 0,01$.
Если количество населения в 1995 году составит 11,5 млн. чел., каково ожидаемое количество врачей? Рассчитайте 99%-й доверительный интервал для данного прогноза.
Рассчитайте коэффициент детерминации

Видео:Метод наименьших квадратов. Линейная аппроксимацияСкачать

Решение систем линейных уравнений методом наименьших квадратов

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

в матричной форме имеющая вид

Если система (8.70) описывает реальный физический процесс, но оказывается несовместной, то естественно ввести коррективы в математическую модель и, в частности, изменить трактовку понятия «решение». Например, часто возникает задача найти такие значения неизвестных ад, Х2, . х_п, при которых функция

где х = (ад, ад. х_п) т , принимает наименьшее значение. Этот подход к понятию решения линейной системы составляет суть метода наименьших квадратов. При этом найденный набор значений х = (ад, ад, . х_п) Т называют псевдорешением или обобщенным решением системы. Очевидно, что если система Ах = b совместна, то ее псевдорешение представляет собой обычное решение, а любое решение системы можно искать как ее псевдорешение.

Описанный подход к понятию решения линейной системы объясняется следующими соображениями. Пусть некоторый процесс описывается некоторым набором величин од, 2 имела наименьшее значение, длина вектора z должна быть минимальной. Это возможно, лишь когда в равенстве Ах + z = Ъ вектор Ах = А х + А^ Х2 +. + +А_пх_п является ортогональной проекцией (см. разд. 8.5) вектора b на пространство Ь(А) столбцов А, А2, . А_п матрицы А, а вектор z — ортогональной составляющей вектора Ъ. Поэтому вектор 2 должен быть ортогональным к каждому вектору-столбцу Aj матрицы А. Следовательно, (2, А,) = 0, j = 1, 2. п. В матричной форме это означает, что A* z = 0. Отсюда, умножив равенство А х + z = b слева на матрицу А*, получим систему

В действительном (вещественном) случае к системе (8.73) можно прийти также используя методы математического анализа, а именно, приравняв к нулю дифференциал

функции

или, что то же самое, приравняв к нулю первые частные производные ПО Х, Х2, . х_п функции F<x).

Примечание. При вычислении дифференциала dF(x) использовано то обстоятельство, что вещественное выражение (Ax — b)* A dx, являясь матрицей первого порядка, не меняется при транспонировании и комплексном сопряжении и потому совпадает с выражением dx* А* (А х — Ь).

Систему А* Ах = А* b называют системой нормальных уравнений. Ее матрица А* А квадратная п-го порядка. Это матрица Грама системы столбцов матрицы А.

Теорема 8.40. Система А* Ах = А* b всегда совместна, причем ее решения, и только они, являются псевдорешениями системы Ах = Ь.

О Сначала убедимся в совместности системы нормальных уравнений. Для этого докажем, что вектор .то = А + Ъ, где А + — псев- дообратная к А матрица, является решением этой системы. С этой целью возьмем какое-нибудь скелетное разложение А = В С матрицы А. Тогда А* = С* В*, Л+ = С* (СС*)- 1 (В* Б)» 1 В* и

Следовательно, Хо = А + b — решение системы нормальных уравнений, а эта система совместная.

В совместности системы нормальных уравнений можно убедиться и следующим образом. Представим столбец b в виде Ь = Ь _<]+ /Д, где

А-, — столбцы матрицы А, Ь 1 — — ортогональная составляющая вектора Ъ. Тогда

Таким образом, система нормальных уравнений равносильна системе А* Ах = А* Ьо, а последняя всегда совместная, поскольку при х = (ki,k₂, . k_n) T

Перейдем к доказательству второй части утверждения теоремы. Пусть х — решение системы нормальных уравнений А* Ах = А* Ь. Эта система равносильна системе А* (А х — Ь) = 0, означающей, что каждый столбец матрицы А ортогонален столбцу г = Ах—b. Следовательно, столбец А х является проекцией вектора b на подпространство Ь(А), порожденное столбцами матрицы А, а вектор х есть псевдорешение системы Ах = Ь.

Эти рассуждения можно провести и в обратном порядке (см. вывод системы нормальных уравнений). Тогда придем к заключению, что псевдорешения системы Ах = b являются решениями системы нормальных уравнений. ?

Следствис8.6. Общее псевдорешение системы Ах = b представимо в виде

где .то — какое-либо частное решение системы А* Ах — А* 6, т.е. какое- либо частное псевдорешение системы А х = 6, а т₀б_Щ — общее решение однородной системы А* А х = 0.

Следствие 8.7. Система А* Ах = А* b имеет единственное решение, если столбцы матрицы А линейно независимы, и бесконечное множество решений, если столбцы матрицы А линейно зависимы.

О Определитель системы А* Ах = А* b есть определитель матрицы Грама системы столбцов матрицы А. Этот определитель отличен от нуля тогда и только тогда, когда столбцы матрицы А линейно независимы. В то же время это условие означает, что система А* Ах = А* b имеет единственное решение. ?

Еще и еще раз подчеркнем, что решения х системы А* Ах = А* b нормальных уравнений, и только они, являются псевдорешениями системы Ах = Ь. Причем, при таких х, и только при них, вектор Ах = Ах 1 + А2Х2 + . + А_пх_п является (см. раздел 8.5) ортогональной проекцией вектора b на подпространство Ь(А), порожденное столбцами А, А2, . А_п матрицы А, а вектор z = Ах — b — ортогональной составляющей вектора Ъ.

Итак, чтобы решить систему Ах = b методом наименьших квадратов, нужно найти общее решение системы А* Ах = А* b и, если требуется, из множества решений выделить нормальное решение, т.е. то, которое имеет наименьшую длину.

Пример 8.36. Для системы

найти общее и нормальное псевдорешения.

Решение. Составим систему А* Ах — А* Ъ нормальных уравнений, которая в данном случае имеет вид:

Решив эту систему, найдем общее псевдорешение Теперь составим функцию

Из равенства нулю ее производной по Х найдем ж 4 = 1, при котором из общего псевдорешения х получается нормальное псевдорешение х° = (3, —5, 3,1) т .

Если для выделения нормального псевдорешения применить второй способ, то сначала нужно построить какую-либо ФСР однородной системы Ах = 0. В данном случае ФСР состоят из одного решения. Одна из этих ФСР состоит из решения X = ( — 1, —1, —1,1) т . Далее из соотношения ортогональности

векторов X и х найдем значение Х4 = 1, при котором из общего псевдорешения х выделяется нормальное псевдорешение х° = (3, —5,3,

Пример 8.37. Найти многочлен p(t) = х + X2t—x^t 2 , с наименьшей квадратичной погрешностью приближающий функцию /(?), которая задана таблицей

Видео:Метод наименьших квадратов, урок 1/2. Линейная функцияСкачать

Методы наименьших квадратов без слёз и боли

Итак, очередная статья из цикла «математика на пальцах». Сегодня мы продолжим разговор о методах наименьших квадратов, но на сей раз с точки зрения программиста. Это очередная статья в серии, но она стоит особняком, так как вообще не требует никаких знаний математики. Статья задумывалась как введение в теорию, поэтому из базовых навыков она требует умения включить компьютер и написать пять строк кода. Разумеется, на этой статье я не остановлюсь, и в ближайшее же время опубликую продолжение. Если сумею найти достаточно времени, то напишу книгу из этого материала. Целевая публика — программисты, так что хабр подходящее место для обкатки. Я в целом не люблю писать формулы, и я очень люблю учиться на примерах, мне кажется, что это очень важно — не просто смотреть на закорючки на школьной доске, но всё пробовать на зуб.

Итак, начнём. Давайте представим, что у меня есть триангулированная поверхность со сканом моего лица (на картинке слева). Что мне нужно сделать, чтобы усилить характерные черты, превратив эту поверхность в гротескную маску?

В данном конкретном случае я решаю эллиптическое дифференциальное уравнение, носящее имя Симеона Деми Пуассона. Товарищи программисты, давайте сыграем в игру: прикиньте, сколько строк в C++ коде, его решающем? Сторонние библиотеки вызывать нельзя, у нас в распоряжении только голый компилятор. Ответ под катом.

На самом деле, двадцати строк кода достаточно для солвера. Если считать со всем-всем, включая парсер файла 3Д модели, то в двести строк уложиться — раз плюнуть.

Видео:Решение системы уравнений методом ГауссаСкачать

Пример 1: сглаживание данных

Давайте расскажу, как это работает. Начнём издалека, представьте, что у нас есть обычный массив f, например, из 32 элементов, инициализированный следующим образом:

А затем мы тысячу раз выполним следующую процедуру: для каждой ячейки f[i] мы запишем в неё среднее значение соседних ячеек: f[i] = ( f[i-1] + f[i+1] )/2. Чтобы было понятнее, вот полный код:

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

Если вам неясно, почему происходит сглаживание, прямо сейчас остановитесь, возьмите ручку и попробуйте порисовать примеры, иначе дальше читать не имеет смысла. Триангулированная поверхность ничем принципиально от этого примера не отличается. Представьте, что для каждой вершины мы найдём соседние с ней, посчитаем их центр масс, и передвинем нашу вершину в этот центр масс, и так десять раз. Результат будет вот таким:

Разумеется, если не остановиться на десяти итерациях, то через некоторое время вся поверхность сожмётся в одну точку ровно так же, как и в предыдущем примере весь массив стал заполнен одним и тем же значением.

Видео:Метод наименьших квадратов. Квадратичная аппроксимацияСкачать

Пример 2: усиление/ослабление характеристических черт

Полный код доступен на гитхабе, а здесь я приведу самую важную часть, опустив лишь чтение и запись 3Д моделей. Итак, триангулированная модель у меня представлена двумя массивами: verts и faces. Массив verts — это просто набор трёхмерных точек, они же вершины полигональной сетки. Массив faces — это набор треугольников (количество треугольников равно faces.size()), для каждого треугольника в массиве хранятся индексы из массива вершин. Формат данных и работу с ним я подробно описывал в своём курсе лекций по компьютерной графике. Есть ещё третий массив, который я пересчитываю из первых двух (точнее, только из массива faces) — vvadj. Это массив, который для каждой вершины (первый индекс двумерного массива) хранит индексы соседних с ней вершин (второй индекс).

Первое, что я делаю, это для каждой вершины моей поверхности считаю вектор кривизны. Давайте проиллюстрируем: для текущей вершины v я перебираю всех её соседей n1-n4; затем я считаю их центр масс b = (n1+n2+n3+n4)/4. Ну и финальный вектор кривизны может быть посчитан как c=v-b, это не что иное, как обычные конечные разности для второй производной.

Непосредственно в коде это выглядит следующим образом:

Ну а дальше мы много раз делаем следующую вещь (смотрите предыдущую картинку): мы вершину v двигаем v := b + const * c. Обратите внимание, что если константа равна единице, то наша вершина никуда не сдвинется! Если константа равна нулю, то вершина заменяется на центр масс соседних вершин, что будет сглаживать нашу поверхность. Если константа больше единицы (заглавная картинка сделана при помощи const=2.1), то вершина будет сдвигаться в направлении вектора кривизны поверхности, усиливая детали. Вот так это выглядит в коде:

Кстати, если меньше единицы, то детали будут наоборот ослабляться (const=0.5), но это не будет эквивалентно простому сглаживанию, «контраст картинки» останется:

Обратите внимание, что мой код генерирует файл 3Д модели в формате Wavefront .obj, рендерил я в сторонней программе. Посмотреть получившуюся модель можно, например, в онлайн-вьюере. Если вам интересны именно методы отрисовки, а не генерирование модели, то читайте мой курс лекций по компьютерной графике.

Видео:Пример решения задачи методом наименьших квадратовСкачать

Пример 3: добавляем ограничения

Давайте вернёмся к самому первому примеру, и сделаем ровно то же самое, но только не будем переписывать элементы массива под номерами 0, 18 и 31:

Остальные, «свободные» элементы массива я инициализировал нулями, и по-прежнему итеративно заменяю их на среднее значение соседних элементов. Вот так выглядит эволюция массива на первых ста пятидесяти итерациях:

Вполне очевидно, что на сей раз решение сойдётся не к постоянному элементу, заполняющему массив, а к двум линейным рампам. Кстати, действительно ли всем очевидно? Если нет, то экспериментируйте с этим кодом, я специально привожу примеры с очень коротким кодом, чтобы можно было досконально разобраться с происходящим.

Видео:Метод наименьших квадратов. ТемаСкачать

Лирическое отступление: численное решение систем линейных уравнений.

Пусть нам дана обычная система линейных уравнений:

Её можно переписать, оставив в каждом из уравнений с одной стороны знака равенства x_i:

Пусть нам дан произвольный вектор , приближающий решение системы уравнений (например, нулевой).

Тогда, воткнув его в правую часть нашей системы, мы можем получить обновлённый вектор приближения решения .

Чтобы было понятнее, x1 получается из x0 следующим образом:

Повторив процесс k раз, решение будет приближено вектором

Давайте на всякий случай запишем рекурретную формулу:

При некоторых предположениях о коэффициентах системы (например, вполне очевидно, что диагональные элементы не должны быть нулевыми, т.к. мы на них делим), данная процедура сходится к истинному решению. Эта гимнастика известна в литературе под названием метода Якоби. Конечно же, существуют и другие методы численного решения систем линейных уравнений, причём значительно более мощные, например, метод сопряжённых градиентов, но, пожалуй, метод Якоби является одним из самых простых.

Видео:11 4 Применение МНК к решению систем линейных уравненийСкачать

Пример 3 ещё раз, но уже с другой стороны

А теперь давайте ещё раз внимательно посмотрим на основной цикл из примера 3:

Я стартовал с какого-то начального вектора x, и тысячу раз его обновляю, причём процедура обновления подозрительно похожа на метод Якоби! Давайте выпишем эту систему уравнений в явном виде:

Потратьте немного времени, убедитесь, что каждая итерация в моём питоновском коде — это строго одно обновление метода Якоби для этой системы уравнений. Значения x[0], x[18] и x[31] у меня зафиксированы, соответственно, в набор переменных они не входят, поэтому они перенесены в правую часть.

Итого, все уравнения в нашей системе выглядят как — x[i-1] + 2 x[i] — x[i+1] = 0. Это выражение не что иное, как обычные конечные разности для второй производной. То есть, наша система уравнений нам просто-напросто предписывает, что вторая производная должна быть везде равна нулю (ну, кроме как в точке x[18]). Помните, я говорил, что вполне очевидно, что итерации должны сойтись к линейным рампам? Так именно поэтому, у линейной функции вторая производная нулевая.

А вы знаете, что мы с вами только что решили задачу Дирихле для уравнения Лапласа?

Кстати, внимательный читатель должен был бы заметить, что, строго говоря, у меня в коде системы линейных уравнений решаются не методом Якоби, но методом Гаусса-Зейделя, который является своебразной оптимизацией метода Якоби:

Видео:Построение уравнения линейной регрессии методом наименьших квадратов.Скачать

Пример 4: уравнение Пуассона

А давайте мы самую малость изменим третий пример: каждая ячейка помещается не просто в центр масс соседних ячеек, но в центр масс плюс некая (произвольная) константа:

В прошлом примере мы выяснили, что помещение в центр масс — это дискретизация оператора Лапласа (в нашем случае второй производной). То есть, теперь мы решаем систему уравнений, которая говорит, что наш сигнал должен иметь некую постоянную вторую производную. Вторая производная — это кривизна поверхности; таким образом, решением нашей системы должна стать кусочно-квадратичная функция. Проверим на дискретизации в 32 сэмпла:

При длине массива в 32 элемента наша система сходится к решению за пару сотен итераций. А что будет, если мы попробуем массив в 128 элементов? Тут всё гораздо печальнее, количество итераций нужно уже исчислять тысячами:

Метод Гаусса-Зейделя крайне прост в программировании, но малоприменим для больших систем уравнений. Можно попытаться его ускорить, применяя, например, многосеточные методы. На словах это может звучать громоздко, но идея крайне примитивная: если мы хотим решение с разрешением в тысячу элементов, то можно для начала решить с десятью элементами, получив грубую аппроксимацию, затем удвоить разрешение, решить ещё раз, и так далее, пока не достигнем нужного результата. На практике это выглядит следующим образом:

А можно не выпендриваться, и использовать настоящие солверы систем уравнений. Вот я решаю этот же пример, построив матрицу A и столбец b, затем решив матричное уравнение Ax=b:

А вот так выглядит результат работы этой программы, заметим, что решение получилось мгновенно:

Таким образом, действительно, наша функция кусочно-квадратичная (вторая производная равна константе). В первом примере мы задали нулевую вторую производную, в третьем ненулевую, но везде одинаковую. А что было во втором примере? Мы решили дискретное уравнение Пуассона, задав кривизну поверхности. Напоминаю, что произошло: мы посчитали кривизну входящей поверхности. Если мы решим задачу Пуассона, задав кривизну поверхности на выходе равной кривизне поверхности на входе (const=1), то ничего не изменится. Усиление характеристических черт лица происходит, когда мы просто увеличиваем кривизну (const=2.1). А если const Скрытый текст

Это результат по умолчанию, рыжий Ленин — это начальные данные, голубая кривая — это их эволюция, в бесконечности результат сойдётся в точку:

А вот результат с коэффицентом 2.:

Домашнее задание: почему во втором случае Ленин сначала превращается в Дзержинского, а затем снова сходится к Ленину же, но большего размера?

Видео:Метод Крамера за 3 минуты. Решение системы линейных уравнений - bezbotvyСкачать

Заключение

Очень много задач обработки данных, в частности, геометрии, могут формулироваться как решение системы линейных уравнений. В данной статье я не рассказал, как строить эти системы, моей целью было лишь показать, что это возможно. Темой следующей статьи будет уже не «почему», но «как», и какие солверы потом использовать.

Кстати, а ведь в названии статьи присутствуют наименьшие квадраты. Увидели ли вы их в тексте? Если нет, то это абсолютно не страшно, это именно ответ на вопрос «как?». Oставайтесь на линии, в следующей статье я покажу, где именно они прячутся, и как их модифицировать, чтобы получить доступ к крайне мощному инструменту обработки данных. Например, в десяток строк кода можно получить вот такое:

📺 Видео

12. Решение систем линейных уравнений методом ГауссаСкачать

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

Метод наименьших квадратов (МНК)Скачать

Метод Наименьших Квадратов (МНК)Скачать

Решение системы линейных уравнений с двумя переменными способом подстановки. 6 класс.Скачать

Система с тремя переменнымиСкачать

Матричный метод решения систем уравненийСкачать

ПОСМОТРИ это видео, если хочешь решить систему линейных уравнений! Метод ПодстановкиСкачать

Cистемы уравнений. Разбор задания 6 и 21 из ОГЭ. | МатематикаСкачать

Решение систем линейных уравнений, урок 5/5. Итерационные методыСкачать