5.21. Однооткликовая задача наименьших квадратов

Пусть у — доля химического вещества А, оставшаяся к моменту времени в результате реакции первого порядка

Переменная у удовлетворяет дифференциальному уравнению

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

Константа скорости зависит от абсолютной температуры следующим образом:

где - предэкспоненциальный множитель, а -энергия активации (выраженная в соответствующих единицах). Уравнение модели принимает вид

Данные пятнадцати измерений х и у приведены в табл. 5.2.

Таблица 5.2 (см. скан) Данные для задачи наименьших квадратов

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

Следуя найдем составляющие вектор-градиента Ф:

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

Пусть в качестве начального приближения взято следующее.

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

Таблица 5.3. (см. скан) Значения функций для задачи наименьших квадратов при [750, 1200]

Чтобы определить направление первого шага необходимо вычислить решить систёму уравнений В нашей двумерной задаче это можно сделать элементарно на настольном вычислителе. Однако для целей иллюстративности мы применим метод Гринштадта. Для этого нам необходимо вычислить обратное шкалированное разложение матрицы N (индекс 1 опускаем). Будем следовать алгоритму, описанному в разделе

Отсюда

2. Матрица С имеет вид Читатель может проверить, что такая матрица имеет собственные значения и собственные векторы соответственно. Следовательно, преобразование для С таково:

3. Преобразование для матрицы, обратной выражается через где

Следовательно,

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

составляло до или поэтому нет необходимости заботиться об увеличении минимальною собственного значения.

4. Чтобы рассчитать проделаем следующие вычисления:

Производная по направлению при равна:

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

где Это говорит о том, что точка приемлема. Попытаемся найти лучшую точку, приближая параболой. Уравнение параболы таково:

В точках искомая парабола должна удовлетворять условиям: 4

Используя получим, что имеет стационарную точку при

Подставляя в найдем, что что уже значительно лучше Примем теперь в качестве исходной точки для следующей итерации.

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

Таблица 5.4 (см. скан) Задача наименьших квадратов, хорошее начальное приближение (процедура Гаусса)

В этой точке вектор-градиент имеет составляющие

а приближенная матрица Гессе имеет вид

Применяя тест из раздела 5.15, найдем, что

Эти значения пренебрежимо малы по сравнению с 0? и следовательно, можно полагать, что наша процедура сошлась к стационарной точке. Остатки, соответствующие этому решению, приведены в табл. 5.5.

Таблица 5.5 (см. скан) Задача наименьших квадратов; остатки, полученные в конце итерациоиой процедуры

В предыдущем примере мы достаточно удачно выбрали начальное приближение; по сравнению с окончательной оценкой Предположим теперь, что мы начали

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

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

при Так как это значение превышает неприемлемо и мы должны прибегнуть к интерполяции. Здесь

Следовательно, из

Поскольку меньше, чем согласно алгоритму рис. 5.2а положим

что

Это значение также неприемлемо. Повторение процедуры интерполяции приведет к выбору

и

Снова необходима интерполяция:

На этот раз так что берем и

Рис. 5.3. (см. скан) Задача наименьших квадратов

Наконец, нами получена приемлемая точка. Положим и перейдем к следующей итерации. Процедура сходится, хотя и довольно медленно, за 25 итераций, что изображено на рис. 5.3. Тут же показано направление отрицательного градиента в точке Так как вектор, направленный из в , лежит между направлением отрицательного градиента и направлением по методу Гаусса от то из

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

Полагая получим размер шага, при котором значение 0 таково, что невозможно вычислить; аналогично и для При однако, найдем:

Это все еще неприемлемо, поэтому увеличиваем до 10 и далее получая:

что является приемлемой стартовой точкой для следующей итерации. Итерационный процесс, который сходится за 10 итераций, также изображен на рис. 5.3.

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

или

где

Теперь наша модель линейна по параметрам которые можно оценить линейным методом наименьших квадратов:

что соответствует

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