РОБАСТНОЕ ОЦЕНИВАНИЕ
ПРОВЕРКА АДЕКВАТНОСТИ МОДЕЛИ
Принцип расчета параметров. Критериальная функция. Робастные оценки
Расчет неизвестных параметров методами Ньютона и Гаусса-Ньютона:
Вычислительная схема алгоритмов
Проверка адекватности модели:
Остаточная дисперсия и критерий c 2
Асимметрия, эксцесс, средний остаток и среднее модулей остатков
Оценка погрешностей и коэффициенты корреляции параметров
Выявление избыточности модели при анализе матрицы Якоби
1) логарифмы р неизвестных констант равновесия lg
bi;2) факторы интенсивности первых
W реагентов в стехиометрической матрице для каждой из L аналитических позиций. В методе спектрофотометрии факторы интенсивности – молярные коэффициенты поглощения, аналитические позиции – длины волн.Общее число искомых параметров z = p +
L ґ W .В программе CLINP логарифмы уточняемых констант устойчивости комплексов (химических форм, species) L
i рассматриваются как функции подлежащих расчету параметров LKu (их количество р равно числу неизвестных констант bi):,
где lg
bi0 – неменяемый вклад в lg bi, LKu – искомые параметры, tiu – элементы матрицы T размером pґ p преобразования параметров LKu в lg bi.См. раздел
ДанныеКРИТЕРИАЛЬНАЯ ФУНКЦИЯ. РОБАСТНЫЕ ОЦЕНКИ
В качестве оценки искомых параметров принимают такой вектор
q (набор подлежащих определению bi, ali), который обращает в минимум выбранную критериальную функцию U:|q > = arg min U(q).
Если измеренные величины A
lk независимы и известен закон распределения их погрешностей измерений e , то вид минимизируемого критерия выбирают, исходя из метода максимального правдоподобия.Если распределение
e подчиняется закону Гаусса, то оценки q находят нелинейным методом наименьших квадратов (МНК), минимизируя функционал U(q), задаваемый формулой,
где w
lk – статистический вес, связанный с оценкой дисперсии s2(Dlk), Dlk – невязка между вычисленной и измеренной величинами некоторого свойства равновесной системы:D
lk = Alkвычислено - Alkизмерено.При отклонении распределения погрешностей
e от нормального оценки метода наименьших квадратов (МНК-оценки) являются смещенными, несостоятельными и неэффективными. Способы оценивания параметров, мало чувствительные к нарушениям стандартных предположений, лежащих в основе МНК, называют устойчивыми (робастными).Конструируя робастный критерий для параметрической идентификации, предполагаем распределение
e с хвостами, более тяжелыми, чем у нормального распределения. Наличие таких хвостов делает более вероятным появление погрешностей e , значительно отличающихся от нуля. Количественной мерой длины хвостов выступает эксцесс g2. Для нормального распределения g2 = 0, для распределений с тяжелыми хвостами g2 > 0.Пусть погрешности измерений
e – независимые случайные величины с плотностью распределенияр(
e) = [(100 - d)Ч j (e) + d Ч h(e)] / 100,где
j (e) – плотность нормального распределения с нулевым средним и единичной дисперсией (функция плотности вероятности основных наблюдений), h(e) – плотность вероятности грубых промахов, d – интенсивность промахов, %. Такая гипотеза об экспериментальных погрешностях носит название “модели грубых промахов”. В этом случае, как показал П.Хьюбер [Хьюбер П. Робастность в статистике. М.: Мир, 1984], оценка максимального правдоподобия соответствует минимуму критерия,
где
xk – взвешенная невязка wk1/2Ч Dk, r(xk, q) – некоторая “функция потерь” (формула записана для случая единственной аналитической позиции, L = 1). Оценки, соответствующие минимуму указанного функционала, называют М-оценками. Они обладают свойствами несмещенности, состоятельности и высокой эффективности, если функция r(x; q) растет медленнее, чем x2. При этом М-оценки менее чувствительны к грубым выбросам, чем МНК-оценки. Хьюбер ввел функцию потерь,
где c – константа, зависящая от интенсивности грубых промахов
d и определяемая как решение уравнения.
М-оценки Хьюбера сохраняют свои оптимальные свойства при малом числе опытов и устойчивы к несимметричным засорениям нормального распределения погрешностей.
Критерий M представляет собой гибрид критериальных функций, соответствующих задачам МНК и метода наименьших модулей (МНМ): при с =
Ґ (интенсивность промахов d = 0) функция М переходит в,
при
d ® 100 % (параметр c ® 0) функция М эквивалентна критерию.
М-оценки не инвариантны относительно преобразования масштаба. Чтобы обеспечить такую инвариантность и обрабатывать данные с произвольной дисперсией функции распределения основных наблюдений, решают минимизационную задачу:
,
одновременно с параметрами
q оценивая и неизвестный масштабный параметр s,,
где z – количество искомых параметров,
j(x) – функция плотности вероятности стандартного нормального распределения.Величина
s2 – дисперсия функции распределения основных наблюдений (оценка остаточной дисперсии в отсутствие загрязнений):.
Задачу параметрической идентификации легко обобщить на многомерный случай, когда в каждой экспериментальной точке одновременно регистрируют
L свойств равновесной системы. Тогда минимизируемый критерий имеет вид,
.
Анализ данных с помощью М-оценок обладает высокой гибкостью: по сравнению с методом наименьших квадратов появляется новая переменная – процент промахов
d . Варьирование d (априорная информация о его величине отсутствует) можно рассматривать как адаптацию процедур параметрической идентификации к обрабатываемым данным.Расчет неизвестных параметров методами Ньютона и Гаусса-Ньютона
Вычислительная схема алгоритмов
Точечные оценки неизвестных параметров |lg bi>, <ali> находим как величины, обращающие в минимум критериальную функцию. Наряду с ними вычисляется и масштабный параметр s .
Для минимизации
Итеративную вычислительную схему уточнения текущих значений искомых параметров lg
bc (с – от current) в методах Ньютона и Гаусса-Ньютона (при фиксированной величине s) передает формула:lg b+ = lg bc + l Ч s,
где вектор
s = -H-1 Ч g
задает направление
движения от lg bc к lg b+ – значениям параметров на следующей итерации, Н – матрица Гессе, g – вектор градиента минимизируемой функции, 0<lЈ1 – длина шага. В методе Ньютона матрицу Гессе рассчитывают точно, а в методе Гаусса-Ньютона пользуются ее приближением. На каждой итерации по lg b М-оценки неизвестных факторов интенсивности ali рассчитываем методом линейной регрессии [Бугаевский А.А., Холин Ю.В., Коняев Д.С. Расчет констант комплексообразования по данным спектрофотометрии на персональных ЭВМ // Журн. неорг. химии. 1993. Т.38, № 2. С. 350-356; Мерный С.А., Коняев Д.С., Холин Ю.В. Робастное оценивание параметров в задачах количественного физико-химического анализа // Вестник Харьковского университета. Химические науки. № 2. 1998. С. 112-120].Глобальную сходимость расчетного алгоритма и высокую скорость его локальной сходимости обеспечивают 1) коррекция знаконеопределенной матрицы Гессе, основанная на ее спектральном разложении; 2) движение в направлении отрицательной кривизны при обходе седловых точек; 3) применение стратегии линейного поиска длины шага
l при движении вне окрестностей седловых точек (правило Армио-Гольдстайна [Armijo L., Pasific J. Math. 1966 V. 16, P.1-3] и Гольдстeйном [Goldstein A.A. Constructive Real Analysis, Harper&Row, New York, 1967]) .Упрощенный расчет матрицы Гессе в методе Гаусса-Ньютона приводит к потере локальной сходимости в случаях а) больших невязок при неудачном выборе начального приближения искомых параметров или неудачном выборе набора химических реакций, включенных в модель или б) при существенно нелинейном характере минимизируемого критерия.
Используются два основных критерии останова итераций по lg b .
Первый связан с исследованием градиента минимизируемой функции
g(lg b). При итерациях вычисляется вектор относительного градиента gr(lg b+) с компонентами gr(lgb+i), вычисляемыми по формуле.
где
Итерации могут быть остановлены и по другому критерию: при малой величине всех приращений
l Ч si. Итерации прекращаются, если значениегде typ
b – “типичное значение” искомых параметров (рекомендуется 1 – 10), становится меньше минимально допустимой относительнаой длины шага при итерациях (рекомендуется 10-6 – 10-4).Кроме того, пользователь может ограничить число выполняемых итераций уточнения lg
b, задав максимально допустимое число итераций Itermax. При Itermax = 0 программа рассчитает равновесный состав системы, взвешенные невязки и характеристики адекватности модели для заданных пользователем значений параметров.Остаточная дисперсия и критерий
c2Если взвешенные невязки
xlk = wlk1/2Ч Dlk – независимые случайные нормально распределенные величины с нулевым средним и единичной дисперсией, то остаточная дисперсия,
где f = N
ЧL- z, z – общее число рассчитанных параметров, представляет собой случайную величину, распределенную как c2/f с f степенями свободы и математическим ожиданием 1. Модель следует признать адекватной, если выполняется неравенство,
где
cf2(a) – 100a-процентная точка распределения c2 для f степеней свободы. Это условие является точным для моделей, линейных по параметрам. Его применяют и для проверки адекватности нелинейных по параметрам моделей комплексообразования. При этом считают уровень значимости a заданным лишь приближенно. При параметрической идентицикации модели с использованием М-оценок Хьюбера найденные значения параметров ali и lg bi могут не соответствовать минимуму остаточной дисперсии. Следовательно, необходимо смягчить стандартную процедуру проверки адекватности. Необходимая коррекция критерия достигается при уменьшении f.При отличном от нормального распределении взвешенных невязок
xlk с эксцессом g2 скорректированное число степеней свободы находят по формулеf = (NЧ L - z) Ч {1 + 0.5Ч g2 Ч (NЧ L - z)/NЧ L }-1.
В программе CLINP 2.1, применяя вместо
g2 его выборочную оценку, находят f и отбрасывают дробную часть.Асимметрия, эксцесс, средний остаток и среднее модулей остатков
По величинам взвешенных невязок
xlk можно найти выборочные оценки параметров их распределения. В программе CLINP 2.1 вычисляются асимметрия,
где
m2 и m3 – соответственно второй (дисперсия) и третий центральные моменты распределения остатков xlk, эксцессg
2 = m4 / (m2)2 - 3,где
m4 – четвертый центральный момент, среднее значение взвешенных остатков (mean residual, average residual)и среднее значение модулей взвешенных остатков (residual mean)
.
Если
x lk подчиняются нормальному распределению с нулевым средним и единичной дисперсией, математические ожидания асимметрии, эксцесса и средней невязки равны нулю, а математическое ожидание среднего значения модулей невязокВыборочные значения
Таблица 1. Процентные точки распределения статистики
Количество наблюдений |
Уровень значимости a |
||||
0.01 |
0.05 |
0.10 |
|||
11 |
0.94 |
0.91 |
0.89 |
||
16 |
0.92 |
0.89 |
0.87 |
||
21 |
0.90 |
0.88 |
0.86 |
||
26 |
0.89 |
0.87 |
0.86 |
||
31 |
0.88 |
0.86 |
0.85 |
||
36 |
0.88 |
0.86 |
0.85 |
||
41 |
0.87 |
0.85 |
0.84 |
||
46 |
0.87 |
0.85 |
0.84 |
||
51 |
0.86 |
0.85 |
0.84 |
||
71 |
0.85 |
0.84 |
0.83 |
||
101 |
0.85 |
0.83 |
0.83 |
||
201 |
0.83 |
0.82 |
0.82 |
||
µ |
0.81 |
0.81 |
0.81 |
Таблица 2. Процентные точки распределения статистики
Количество наблюдений |
Уровень значимости a |
||
0.01 |
0.05 |
||
25 |
1.06 |
0.71 |
|
30 |
0.98 |
0.66 |
|
35 |
0.92 |
0.62 |
|
40 |
0.87 |
0.59 |
|
50 |
0.79 |
0.53 |
|
70 |
0.67 |
0.49 |
|
100 |
0.57 |
0.39 |
|
150 |
0.46 |
0.32 |
|
200 |
0.40 |
0.28 |
|
400 |
0.285 |
0.20 |
|
1000 |
0.18 |
0.13 |
Таблица 3. Процентные точки распределения статистики
g2
Количество наблюдений |
Уровень значимости a |
||
0.01 |
0.05 |
||
50 |
1.92 |
1.01 |
|
100 |
1.40 |
0.77 |
|
200 |
0.98 |
0.57 |
|
300 |
0.79 |
0.47 |
|
500 |
0.60 |
0.37 |
|
1000 |
0.41 |
0.26 |
Оценка погрешностей и коэффициенты корреляции параметров
Программа по окончании итераций рассчитывает симметричную ковариационную матрицу рассчитанных параметров D(lg b). На диагонали ковариационной матрицы стоят оценки дисперсий компонентов вектора lg b – s2(lg bj), а над и под диагональю – их ковариации cov(lg bi, lg bj):
.
Эта матрица используется для определения доверительных интервалов параметров и расчета коэффициентов их корреляции.
Частные коэффициенты корреляции
rij измеряют взаимное влияние параметров lg bi и lg bj при условии, что остальные параметры фиксированы:.
Общие коэффициенты корреляции
sij служат мерой коррелированности параметров lg bi и lg bj, когда другие параметры рассматриваются как уточняемые:.
Множественные коэффициенты корреляции
Ri измеряют степень коррелированности параметра lg bi со всеми остальными параметрами:.
Коэффициенты корреляции принимают значения в интервале [-1;1]. Нулевые значения соответствуют полной независимости параметров, значения
± 1 – линейной зависимости между параметрами (модель переопределена). Считают, что пересмотр химической модели (выведение из нее одного из комплексов) необходим в случае, когда коэффициент корреляции по абсолютной величине превышает пороговое значение, выбранное в интервале 0.95-0.98.Выявление избыточности модели при анализе матрицы Якоби
Программа сообщает о результатах сингулярного разложения матрицы Якоби
J = ¶ A / ¶ lg b = U S VT,
где
U и V – ортогональные матрицы, S – прямоугольная матрица с диагональным и нулевым блоками, на диагонали которой находятся сингулярные числа ki. Малые сингулярные числа (ki / k макс < 10-4 – 10-6) соответствуют линейным комбинациям искомых параметров lg bj, к изменению которых критериальная функция не чувствительна. Коэффициенты этих комбинаций содержатся в столбцах матрицы V, номера которых соответствуют номерам малых сингулярных чисел. Зачастую уже одних этих коэффициентов достаточно, чтобы выявить избыточный комплекс. В сложных случаях устранить избыточность помогает исследование равновесные концентрации реагентов: неразрешимая линейная комбинация логарифмов констант появляется, если в выражения искомых констант входят концентрации реагентов, не представленных в системе.