Метод наименьших квадратов для полиномиальных моделей

порядка

Обозначим алгебраический полином й степени как:

(2.29)

 

Вектор неизвестных параметров имеет вид:

 

(2.30)

 

(2.31)

 

(2.32)

 

(2.33)

 

(2.34)

 

Изменив порядок суммирования, получим уравнение вида:

 

. (2.35)

 

В конечном виде уравнение можно представить как:

 

(2.36)

 

Обозначим как матрицу базисных функций

 

(2.37)

 

Транспонированная матрица базисных функций:

 

(2.38)

 

 

В матричной форме система линейных алгебраических уравнений для метода наименьших квадратов будет иметь вид:

 

(2.39)

 

отсюда следует, что вектор может быть вычислен по формуле:

 

(2.40)

 

Рассмотрим пример линейной регрессионной зависимости:

 

матрица базисных функций. (2.41)

транспонированная матрица

базисных функций. (2.42)

 

Тогда систему линейных алгебраических уравнений для этого примера можно представить в виде:

 

(2.43)

 

,

 

 

(2.44)

 

Решая систему уравнений получим значения коэффи-

циентов и линейного регрессионного уравнения

 

2.4. Нормы вектора и матрицы

 

Нормой вектора называется число, обладающее следующими свойствами:

· Норма только в случае

· Для любого вектора и числа выполняется условие .

· неравенство треугольника.

 

Норма является количественной характеристикой вектора .

Наиболее распространенными нормами являются:

 

(2.45)

 

Вторая норма, называемая Евклидовой, используется наиболее часто.

Нормой матрицы называется величина

(2.46)

Норма матрицы обладает следующими свойствами:

· Норма только в случае

· Для любой матрицы и числа выполняется условие .

· для любых матриц и .

· для любых матриц и .

· для любой матрицы и любого вектора .

Существуют следующие нормы матриц:

(2.47)

Здесь собственные числа матрицы Собственные числа и их вычисление будут рассмотрены ниже.

 

2.5. Обусловленность системы линейных алгебраических уравнений

 

Рассмотрим систему из двух алгебраических уравнений

 

(2.48)

Дадим геометрическую интерпретацию решения этой системы уравнений, рисунок 2.2.

 

Рисунок 2.2. Графическая интерпретация решений систем линейных алгебраических

уравнений

а – единственное решение, ; б – нет решения, ;

в – матрица А плохо обусловлена,

 

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

На практике применения метода наименьших квадратов на ЭВМ может оказаться, что определитель матрицы при сравнительно малых степенях аппроксимирующего полинома близок к нулю [23]. Особенно неблагоприятная ситуация возникает в случае, когда точки равномерно расположены на отрезке . Без потери общности можно предположить, что Тогда элементы матрицы системы нормальных уравнений имеют вид:

 

(2.49)

 

Определитель матрицы системы уравнений равен:

 

(2.50)

где определитель матрицы Гильберта порядка.

 

матрица Гильберта. (2.51)

 

При увеличении определитель Гильберта очень быстро стремится к нулю. Определитель Гильберта вычисляется по формуле:

 

(2.52)

 

Матрица Гильберта является плохо обусловленной.

Обусловленность матрицы системы алгебраических уравнений определяется числом обусловленности

 

(2.53)

 

где нормы вектора и матрицы соответственно.

Число обусловленности Чем больше число обусловленности тем матрица ближе к вырожденной ().

Рассмотрим, как влияет изменение правой части системы уравнений на решение . Для этого представим систему уравнений в виде

(2.54)

 

где – изменение вектора решения при изменении вектора правых

частей ;

– изменение (ошибка) вектора правых частей.

Из уравнения (3.4) следует, что

 

(2.55)

 

Из определения числа обусловленности (3.3) следует, что

(2.56)

 

 

Из этих неравенств можно получить неравенство:

 

(2.57)

 

 

где – относительное изменение вектора решения;

-относительное изменение вектора правых частей.

 

Неравенство (2.57) показывает, что относительное изменение вектора правых частей увеличивает относительную ошибку решения

на величину , то есть число обусловленности играет

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

То же самое верно и в отношении изменений в коэффициентах матрицы , то есть

 

(2.58)

В частности число обусловленности единичной матрицы .

Число обусловленности ортогональной матрицы () также равно единице.

Если диагональная матрица, то число обусловленности

 

(2.59)

 

Число обусловленности матрицы Гильберта является величиной порядка где размерность матрицы. Определитель матрицы , обратной к матрице , оценивается величиной Поэтому обратная матрица также плохо обусловлена.

Пример. Дана система линейных алгебраических уравнений

 

(2.60)

 

Точным решением этой системы уравнений будет

 

(2.61)

Изменив правую часть системы уравнений , получим следующее решение:

(2.62)

 

Это произошло из-за того, что , то есть матрица плохо обусловлена.

Плохая обусловленность матрицы приводит к тому, что появление даже самых малых ошибок в правой части системы нормальных уравнений вызывает появление больших ошибок в решении системы. Уже при в расчетах на ЭВМ, оперирующих восьмизначными цифрами, получаются совершенно неверные результаты, что видно из приведенного выше примера. Например, при матрица обратная к матрице Гильберта содержит элементы порядка Поэтому даже ошибка порядка в правых частях и коэффициентах системы нормальных уравнений может привести к ошибке порядка 3 в решении системы.

Влияние плохой обусловленности значительно уменьшается, если вместо полиномиальной модели используется модель вида

(2.63)

 

представляющую собой линейную комбинацию элементов подсистемы базисных функций

, (2.64)

 

где ортогональные полиномы.

Система полиномов

(2.65)

 

является ортогональной, если выполняется условие

 

(2.66)

 

2.6. Применение ортогональных полиномов Чебышева и Форсайта

 

Ортогональные полиномы Чебышева имеют вид:

где

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

матрица базисных функций. (2.68)

 

 

 

транспонированная матрица

базисных функций. (2.69)

 

 

Используя уравнение для ортогональных полиномов, получим

 

(2.70)

 

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

 

(2.71)

 

 

Оценки коэффициентов не коррелированы между собой и имеют дисперсии

(2.72)

 

 

где дисперсия случайных ошибок эксперимента.

Дисперсия , как правило, неизвестна, поэтому можно использовать ее оценку

(2.73)

 

где значения выходной величины, полученные по уравнению регрессии в точках Если модель выбрана правильно, то оценка дисперсии является несмещенной, то есть

При решении практических задач степень аппроксимирующего полинома обычно неизвестна, поэтому выбор степени полинома осуществляется следующим образом.

Аппроксимирующая функция имеет вид

 

(2.74)

 

где - коэффициенты многочлена.

Аппроксимирующую функцию можно представить в виде алгебраического полинома

(2.75)

где

 

Ввод данных (экспериментальная зависимость)

 

 

 

 

 

 

 

Величина остаточной дисперсии указывает на хорошее совпадение аппроксимирующей и табличной зависимостей.

 

 

Процедура вычисления коэффициентов алгебраического полинома

путем преобразования полиномов Чебышева в системе Mathcad

 

Ввод исходных данных (экспериментальная зависимость)

 

 

 

 

 

 

 

Ортогональные полиномы Форсайта

 

Ортогональные полиномы Форсайта имеют вид:

 

(2.76)

где

(2.77)

 

Преобразование полиномов Форсайта осуществляется по формулам (4) - (6), аналогично полиномам Чебышева.

Здесь

 

 

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

(2.79)

Процедура вычисления коэффициентов алгебраического полинома

путем преобразования ортогональных полиномов Форсайта в системе Mathcad

 

 

 

 

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

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

 

 

2.7. Регрессионный анализ для многомерных линейных моделей [3]

 

Пусть задана таблица плана эксперимента с факторами, таблица 2.1.

 

Таблица 2.1

 

С использованием метода наименьших квадратов можно найти коэффициенты уравнения регрессии по данным таблицы 1:

 

где фиктивная переменная, равная единице.

 

 

Представим данные в таблице 1 в матричной форме:

 

вектор наблюдений. (2.81)

 

 

вектор неизвестных коэффициентов

регрессионной модели. (2.82)

 

 

Транспонированная матрица будет иметь вид:

 

(2.83)

 

 

Система нормальных уравнений для определения коэффициентов имеет вид:

 

(2.84)

 

В матричной форме система уравнений имеет вид:

 

(2.85)

 

Перемножив матрицы , получим:

 

(2.86)

 

Матрица называется матрицей моментов.

Вектор правых частей в развернутой форме:

 

(2.87)

 

Вектор коэффициентов определяется из уравнения:

 

(2.88)

Обратная матрица может быть представлена в виде:

 

ковариационная матрица. (2.89)

 

Элементы обратной матрицы представляют собой следующие выражения:

(2.90)

где алгебраическое дополнение элемента в матрице

Вычисление обратной матрицы осуществляется по следующей формуле:

 

(2.91)

 

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

Чтобы существовала обратная матрица должна быть невырожденной. Для этого необходимо, чтобы переменные были линейно независимы, то есть должны выполняться условия

 

(2.92)

 

Это означает, что ни один элемент не является линейной комбинацией других элементов В этом случае в матрице независимых переменных элементы одного столбца не будут линейной комбинацией соответствующих элементов других столбцов.

Пример линейной зависимости двух столбцов:

 

.

 

Здесь второй столбец матрицы является линейной комбинации первого столбца, то есть

 

 

Для вычисления остаточной дисперсии вычисляют значения элементов вектора

 

(2.93)

 

Здесь вектор значений функций отклика, вычисленных по

уравнению регрессии.

Числитель остаточной дисперсии получается умножением матриц

 

(2.94)

 

 

Обозначим через вектор-столбец коэффициентов истинной регрессии, при этом математическое ожидание вектора оценок коэффициентов регрессии равно Оценки коэффициентов

Дисперсии коэффициентов регрессии определяются по формуле:

 

(2.95)

 

Коэффициенты являются элементами ковариационной матрицы, стоящими на главной диагонали.

Тогда математическое ожидание определяется выражением:

(2.96)

 

где генеральная дисперсия

ковариационный момент между коэффициентами и

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

В результате получим уравнение:

(2.97)

где случайный нормальный вектор с независимыми элементами, имеющими дисперсии

(2.98)

 

Примечание. Ковариация между элементами и векторов определяется следующим соотношением:

 

(2.99)

 

Матрица коэффициентов системы нормальных уравнений является симметричной и поэтому

 

(2.100)

 

Полагая и, учитывая статистическую независимость, получим:

 

(2.101)

 

Таким образом, имеем зависимость:

 

(2.102)

 

Отсюда получаем следующие зависимости:

 

(2.103)

 

Матрица называется матрицей ошибок или ковариационной матрицей. Так как ковариационная матрица в общем случае не является диагональной, коэффициенты статистически взаимосвязаны. Это приводит к тому, что невозможно проверить значимость каждого коэффициента в отдельности. Поэтому отношение

 

(2.104)

 

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

 

 

2.8. Аппроксимация экспериментальных зависимостей кубическими

сплайнами [23]

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

Рисунок 2.3. Сплайн

 

Вторая производная ( сплайн - функция) примерно равна кривизне. Это означает, что кривизна такого сплайна пропорциональна интегралу

. (2.105)

 

 

Сплайн-функция , для которой интеграл (2.105) имеет минимальное значение. Эта функция является полиномом третьей степени. Каждый отрезок между двумя соседними узлами интерполяции описывается своей сплайн – функцией . Соседние полиномы соединяются непрерывно.

Если на границах всего диапазона [ ] , то такой сплайн называется естественным кубическим сплайном. С математической точки зрения естественный кубический сплайн является самой гладкой функцией, интерполирующей заданные точки.

Иногда вместо естественных граничных условий используются некоторые другие условия вблизи и на границах табличной зависимости. На рисунке 2.4 показана графическая иллюстрация интерполяции табличной функции кубическими сплайнами.

 

Рисунок 2.4. Графическая иллюстрация интерполяции кубическими сплайнами

 

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

 

(2.106)

 

Из соотношения следует, что:

 

(2.107)

 

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

(2.108)

 

Функция непрерывно дифференцируемая, поэтому в любой точке интервала левосторонняя и правосторонняя производные равны между собой.

Первая производная от функции сплайна имеет вид:

 

(2.109)

 

Коэффициенты сплайна являются неизвестными, поскольку неизвестна вторая производная функции

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

 

(2.110)

(2.111)

 

должны быть равны, поэтому можно записать следующее уравнение:

 

(2.112)

 

где .

 

Из уравнения можно получить систему линейных алгебраических уравнений вида:

,

где . (2.113)

 

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

(2.114)

 

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

 

 

 

(2.115)

 

 

Значения вычисляются по формулам:

 

(2.116)

(2.117)

 

Система уравнений имеет решение, поскольку в ее матрице диагональные элементы доминируют (преобладают) над элементами, не лежащими на главной диагонали, то есть

Пример. В статье [20] приведена зависимость адсорбции поливинилового спирта от концентрации его раствора , полученная по результатам эксперимента при производстве декоративно-защитных материалов на бумажной основе таблица 2.2.

 

Таблица 2.2

Вычисление кубического сплайна в системе Mathcad

 

Ввод исходных данных

 

Ниже приведены графики экспериментальной и интерполирующей зависимостей с применением встроенных функций Spline и Interp.

 

Здесь точками обозначена экспериментальная зависимость , сплошной линией – интерполирующая зависимость. Как видно из графиков экспериментальная и интерполирующая зависимости полностью совпадают.

Глава 3

Планирование многофакторных экспериментов

Факторы

Фактором называется измеряемая величина, принимающая в некоторый момент времени определенное значение. Факторы соответствуют способам воздействия на объект исследования.