Базовые сведения

Постановка задачи

Суть задачи калибровки состоит в следующем. Пусть имеется некоторая переменная (свойство) $y$, величину которой необходимо установить. Однако, по ряду обстоятельств (недоступность, дороговизна, длительность), прямое измерение величины y невозможно. В то же время можно легко (быстро, дешево) измерить другие величины: $\mathbf{x}=(x_1, x_2, x_3,…)$, которые связаны с искомой величиной $y$. Задача калибровки состоит в установлении количественной связи между переменными $\mathbf{x}$ и откликом $y$ –

$$y=f(x_1, x_2, x_3, \dots | a_1, a_2, a_3, \dots) + \epsilon$$

На практике это означает:

1. подбор вида зависимости $f$, и ...
2. оценку неизвестных параметров $a_1, a_2, a_3, \dots$ в этой калибровочной зависимости.

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

Простейший пример калибровки дает общеизвестный прибор, называемый безменом, т.е. пружинные весы. Искомая величина $y$ – это вес образца, а $x$ – это удлинение пружины весов. Закон Гука —

$$y = ax + ε$$

связывает $y$ и $x$ через коэффициент жесткости пружины $a$, который априори неизвестен. Процедура калибровки очень проста – взвешиваем стандартный образец весом в 1 кг и отмечаем на шкале удлинение пружины, затем используем образец в 2 кг, и т.д. В результате этой калибровки (точнее, градуировки) получается шкала, по которой можно определить вес нового, нестандартного образца.

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

Разумеется, на практике все обстоит не так просто, как в этом элементарном примере. Например, в калибровке может участвовать не один показатель $y$ (отклик), а несколько откликов $y_1, y_2,\dots y_K$, которые нас интересуют. Все возможные особенности, различные трудности, сопутствующие процедуре калибровке будут рассмотрены далее. Сейчас же мы подведем первые итоги и сформулируем задачу калибровки в общем виде.

Пусть имеется матрица $\mathbf{Y}$, размерностью ($I×K$), где $I$ – это число стандартных образцов (сравнения), использованных в калибровке, а $K$ – это число одновременно калибруемых откликов. Матрица $\mathbf{Y}$ содержит значения откликов $y$, известные из независимых экспериментов (референтные или стандартные значения). Пусть, с другой стороны, имеется соответствующая матрица переменных $\mathbf{X}$ размерностью ($I×J$), где $I$ – это, естественно, тоже число образцов, а $J$ – это число независимых переменных (каналов), используемых в калибровке. Матрица $\mathbf{X}$ состоит из альтернативных, как правило, многоканальных ($J>>1$) измерений. Используя калибровочные данные $(\mathbf{X}, \mathbf{Y})$, требуется построить функциональную связь между $\mathbf{Y}$ и $\mathbf{X}$.

Итак, задача калибровки состоит в построении математической модели, связывающей блоки $\mathbf{X}$ и $\mathbf{Y}$, с помощью которой можно в дальнейшем предсказывать значения показателей $y_1, y_2, \dots y_K$, по новой строке значений аналитического сигнала $\mathbf{x}$.

Линейная и нелинейная калибровки

Теоретически функциональная связь между переменными $x$ и $y$ может быть сложной, например,

$$y = b_0 \mathrm{exp}(b_1x+b_2x^2) + ε$$

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

$$y = b_0 + b_1x_1 + b_2x_2+ \dots + b_Jx_J + ε$$

Заметим, что термин "линейность" употребляется в контексте этого пособия вместо более правильного термина "билинейность", означающего линейность как по отношению к переменным $x$, так и в отношении коэффициентов $b$. Поэтому, калибровка

$$y = b_0 + b_1x + b_2x^2 + ε$$

является нелинейной, несмотря на то, что она легко "линеаризуется" введением новой переменной $x_2=x^2$.

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

$$\mathbf{Y} = \mathbf{XB} + \mathbf{E}$$

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

Предпочтение (би)линейных, формальных моделей, не отягощенных дополнительным физико-химическим смыслом, является магистральным направлением развития хемометрики. Такой подход позволяет исключить влияние субъективного фактора, проявляющегося при выборе калибровочной зависимости. Однако за это приходится расплачиваться – все линейные модели имеют ограниченную область применения.

Суть дела иллюстрирует Рис. 2, где показаны четыре участка, на которых сложная нелинейная зависимость приближается простыми линейными моделями. Каждая из этих моделей работает только на своей области переменной $x$ и выход за ее границы приводит к грубым ошибкам. Принципиальным моментом здесь является то, какую область можно считать допустимой – иначе говоря, насколько широко можно применять калибровочную модель. Ответ на этот вопрос дают методы проверки (валидации).

Калибровка и проверка

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

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

Некоторую фиксированную долю (например, первые 10% образцов) исключают из исходного набора данных. Затем строят модель, используя только оставшиеся 90% данных, и применяют ее к исключенному набору. На следующем цикле исключенные данные возвращаются, и удаляется уже другая порция данных (следующие 10%), и опять строится модель, которая применяется к исключенным данным. Эта процедура повторяется до тех пор, пока все данные не побывают в числе исключенных (в нашем случае – 10 циклов). Наиболее (но неоправданно) популярен вариант перекрестной проверки, в котором данные исключаются по одному (Leave-One-Out – LOO).

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

"Качество" калибровки

Результатом калибровки являются величины $\mathbf{\hat{Y}}_\mathrm{c}$ — оценки стандартных откликов $\mathbf{Y}_\mathrm{c}$, найденные по модели, построенной на обучающем наборе. Результатом проверки служат величины $\mathbf{\hat{Y}}_\mathrm{t}$ – оценки проверочных откликов $\mathbf{Y}_\mathrm{t}$, вычисленные по той же модели. Отклонение оценки от стандарта (проверочного значения) вычисляют как матрицу остатков: в обучении $\mathbf{E}_\mathrm{c} = \mathbf{Y}_\mathrm{c} - \mathbf{\hat{Y}}_\mathrm{c}$, и в проверке $\mathbf{E}_\mathrm{t} = \mathbf{Y}_\mathrm{t} - \mathbf{\hat{Y}}_\mathrm{t}$.

Следующие величины характеризуют "качество" калибровки в среднем.

A. Полная дисперсия остатков в обучении ($TRVC$) и в проверке ($TRVP$) –

$$TRVC = \frac{1}{KI_\mathrm{c}} \sum_{k=1}^{K} \sum_{i=1}^{I_\mathrm{c}} e^2_{ki}$$ $$TRVP = \frac{1}{KI_\mathrm{t}} \sum_{k=1}^{K} \sum_{i=1}^{I_\mathrm{t}} e^2_{ki}$$

Эти величины вычисляются по формуле указанное в пособии Статистика, для $m = 0$.

B. Объясненная дисперсия остатков при обучении ($ERVC$) и при проверке ($ERVP$) –

$$ ERVC = 1 - \sum_{k=1}^{K} \sum_{i=1}^{I_\mathrm{c}} e^2_{ki} \bigg/ \sum_{k=1}^{K} \sum_{i=1}^{I_\mathrm{c}} y^2_{ki} $$ $$ ERVP = 1 - \sum_{k=1}^{K} \sum_{i=1}^{I_\mathrm{t}} e^2_{ki} \bigg/ \sum_{k=1}^{K} \sum_{i=1}^{I_\mathrm{t}} y^2_{ki} $$

C. Среднеквадратичные остатки калибровки ($RMSEC$) и проверки ($RMSEP$) –

$$RMSEC(k) = \sqrt{\frac{1}{I_\mathrm{c}} \sum_{i=1}^{I_\mathrm{c}} e^2_{ki}}$$ $$RMSEP(k) = \sqrt{\frac{1}{I_\mathrm{t}} \sum_{i=1}^{I_\mathrm{t}} e^2_{ki}}$$

Величины $RMSE$ зависят от $k$ – номера отклика.

D. Смещение в калибровке ($BIASC$) и в проверке ($BIASP$) –

$$BIASC(k) = \frac{1}{I_\mathrm{c}} \sum_{i=1}^{I_\mathrm{c}} e_{ki}$$ $$BIASP(k) = \frac{1}{I_\mathrm{t}} \sum_{i=1}^{I_\mathrm{t}} e_{ki}$$

Величины $BIAS$ зависят от $k$ – номера отклика.

E. Стандартные ошибки в калибровке ($SEC$) и в проверке ($SEP$) –

$$ SEC(k) = \sqrt{\frac{1}{I_\mathrm{c}} \sum_{i=1}^{I_\mathrm{c}} \big( e_{ki} - BIASC(k) \big)^2} $$ $$ SEP(k) = \sqrt{\frac{1}{I_\mathrm{t}} \sum_{i=1}^{I_\mathrm{t}} \big( e_{ki} - BIASP(k) \big)^2} $$

Величины $SE$ зависят от $k$ – номера отклика.

F. Коэффициенты корреляции $R^2$ между стандартными $y_{ki}$ и оцененными $\hat{y}_{ki}$ откликами. Они также вычисляется отдельно для обучающего $R_\mathrm{c}^2(k)$ и проверочного $R_\mathrm{y}^2(k)$ наборов.

$$ R(k) = \frac{ I \sum y_{ki} \hat{y}_{ki} - \sum y_{ki}\sum \hat{y}_{ki} }{ \sqrt{I \sum y^2_{ki} - (\sum y_{ki})^2}
\sqrt{I \sum \hat{y}^2_{ki} - (\sum \hat{y}_{ki})^2}
} $$

Величины $R^2$ зависят от $k$ – номера отклика.

Во всех этих формулах величины $e_{ki}$ – это элементы матриц $\mathbf{E}_\mathrm{c}$ или $\mathbf{E}_\mathrm{t}$. Для характеристик, наименование которых оканчивается на $C$ (например, $RMSEC$), используется матрица $\mathbf{E}_\mathrm{c}$ (обучение), а для тех, которые оканчиваются на $P$ (например, $RMSEP$), берется матрица $\mathbf{E}_\mathrm{t}$ (проверка).

Неопределенность, точность и воспроизводимость

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

Воспроизводимость (precision) характеризует то, насколько близко находятся друг от друга независимые повторные измерения.

Точность (accuracy) определяет степень близости оценок к истинному (стандартному) значению y.

Рис. 6 объясняет суть дела. Графики a) и b) представляют оценки с хорошей воспроизводимостью. Вариант a) , кроме того, обладает и высокой точностью, что, разумеется, нельзя сказать о графике c). В последнем случае и воспроизводимость, и точность оставляют желать лучшего.

Показатели $SEC$ и $SEP$ характеризуют воспроизводимость калибровки, тогда как $RMSEC$ и $RMSEP$ показывают ее точность. Величины $BIASC$ и $BIASP$ определяют смещение калибровки относительно истинного значения (Рис. 6b). Можно показать, что –

$$RMSE^2 = SE^2 + BIAS^2$$

Таким образом, при построении калибровки предпочтение следует отдавать показателям $RMSE$, а не $SE$.

Показатели $TRV$ и $ERV$ характеризуют ситуацию "в целом", без различения калибруемых откликов, т.е.

$$TRV = \frac{1}{K} \sum_{k=1}^{K} RMSE^2(k)$$

Недооценка и переоценка

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

Из курса школьной физики известно, что расстояние $L$, на которое летит снаряд, выпущенный со скоростью $v$ из орудия, наклоненного под углом α к горизонту, равно –

$$L = v^2 \mathrm{sin}(2α)/g$$

Предположим, что у нас имеются экспериментальные данные $L(α)$.

Забудем о том, что функциональная связь между $L$ (т.е., $y$) и $α$ (т.е., $x$) нам известна, и попробуем построить формальную полиномиальную калибровку –

$$L = b_0 + b_1α + b_2α^2 + \dots b_n α^n + ε$$

по этим данным. На Рис. 8а показано, как описываются обучающие (■) и проверочные (○) данные моделями при $n = 0, 1, \dots , 5$. Видно, что при недостаточной сложности модели ($n = 0, 1$) обучающие данные лежат далеко от модели. Это – случай недооценки. Когда сложность модели увеличивается, то согласие между обучающими данными и моделью улучшается. Однако, при излишней сложности ($n = 4, 5$), модель все хуже работает при прогнозе на проверочный набор. Это – случай переоценки.

На Рис. 8b показано, как изменяются показатели $RMSEC$ и $RMSEP$ при увеличении сложности модели. Это – типичный график, в котором $RMSEC$ монотонно падает, а $RMSEP$ проходит через минимум. Именно точка минимума $RMSEP$ позволяет определить оптимальную сложность модели. В нашем случае – это 2. Теперь мы можем вспомнить о содержательной модели, которая хорошо аппроксимируется полиномом второй степени по $α$ –

$$L = v^2/g \mathrm{sin}(2α) ≈ Const α(α – π/4)$$

В задачах многомерной калибровки недооценка и переоценка проявляются через выбор числа скрытых, главных компонент. Когда их число мало, модель плохо приближает обучающий набор и при увеличении сложности модели $RMSEC$ монотонно уменьшается. Однако качество прогноза на проверочный набор может при этом ухудшаться (U-образная форма $RMSEP$). Точка минимума $RMSEP$, или начало плато, соответствует оптимальному числу главных компонент. Проблема сбалансированности описания данных рассматривается во многих работах А. Хоскюлдссона, который в 1988 году ввел новую концепцию моделирования – так называемый H-принцип. Согласно этому принципу точность моделирования ($RMSEC$) и точность прогнозирования ($RMSEP$) связаны между собой. Улучшение $RMSEC$ неминуемо влечет ухудшение $RMSEP$, поэтому их нужно рассматривать совместно. Именно по этой причине множественная линейная регрессия, в которой всегда участвует явно избыточное число параметров, неизбежно приводит к неустойчивым моделям, непригодным для практического применения.

Мультиколлинеарность

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

$$\mathbf{Y} = \mathbf{XB}$$

Разумеется, здесь мы имеем дело с обучающим набором данных. Мы не пишем индекс у матриц ($\mathbf{Y}_\mathrm{с}$, $\mathbf{X}_\mathrm{с}$) только для простоты обозначений. Классическое решение этой задачи находится с помощью метода наименьших квадратов. Минимизируя сумму квадратов отклонений $(\mathbf{Y} – \mathbf{XB})^\mathrm{t}(\mathbf{Y} – \mathbf{XB})$, получаем оценки неизвестных коэффициентов $\mathbf{B}$ –

$$ \hat{\mathbf{B}} = \big( \mathbf{X}^{\mathrm{t}}\mathbf{X} \big)^{-1}\mathbf{X}^\mathrm{t}\mathbf{Y} $$

Далее находится –

$$\hat{\mathbf{Y}} = \mathbf{X}\hat{\mathbf{B}}$$

– оценка искомых откликов. Главная проблема при таком классическом подходе – это обращение матрицы $\mathbf{X}^\mathrm{t}\mathbf{X}$ . Очевидно, что если число стандартных образцов меньше, чем число переменных ($I<J$) то обратной матрицы не существует. Более того, даже при достаточно большом числе образцов ($I>J$), обратной матрицы может и не быть. Рассмотрим простейший пример.

На Рис. 9 представлены данные, в которых всего три образа ($I=3$) и две переменные ($J=2$). Используя активный элемент на листе, угол между векторами $1$ и $2$ можно изменять от $90^о$ до $0^о$. Чем меньше угол, тем меньше детерминант матрицы $\mathbf{X}^\mathrm{t}\mathbf{X}$ и тем хуже определяется матрица $(\mathbf{X}^\mathrm{t}\mathbf{X})^{–1}$. В предельном случае, когда угол между векторами равен $0^о$, матрица $\mathbf{X}^\mathrm{t}\mathbf{X}$ не может быть обращена. При этом векторы 1 и 2 колинеарны, т.е. они лежат на одной прямой, совпадающей с вектором 3. Анимация на Рис. 10 иллюстрирует суть проблемы – невозможность построения калибровки классическим способом, когда данные коллинеарны.

В многомерной калибровке (т.е. при $J>>1$), мы все время сталкиваемся с мультиколлинеарностью – множественными связями между векторами $X$-переменных, приводящими к той же проблеме – невозможности обращения матрицы $\mathbf{X}^\mathrm{t}\mathbf{X}$. Далее мы увидим, что можно сделать в случае мультиколлинеарности данных.

Подготовка данных

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

Центрирование – это вычитание из исходной матрицы $\mathbf{X}$ матрицы $\mathbf{M}$, т.е.

$$\tilde{\mathbf{X}} = \mathbf{X} - \mathbf{M}$$

Обычно усреднение проводится по столбцам: для каждого вектора $\mathbf{x}_j$ вычисляется среднее значение –

$$m_i = \big(x_{1j} + \dots + x_{Ij} \big) \big/ I$$

Тогда $\mathbf{M}=(m_1\mathbf{1},\dots,m_J\mathbf{1})$, где $\mathbf{1}$ – это вектор из единиц размерности $I$.

Центрирование – это почти обязательная процедура перед применением проекционных методов. Второе простейшее преобразование данных, нормирование, не является обязательным. Нормирование, в отличие от центрирования, не меняет структуру данных, а просто изменяет вес различных частей данных при обработке. Чаще всего применяется нормирование по столбцам – это умножение исходной матрицы $\mathbf{X}$ справа на матрицу $\mathbf{W}$, т.е.

$$\tilde{\mathbf{X}} = \mathbf{X}\mathbf{W}$$

Матрица $\mathbf{W}$ – это диагональная матрица размерности $J×J$. Обычно диагональные элементы $w_{jj}$ равны обратным значениям стандартного отклонения –

$$d_i = \sqrt{\sum_{i=1}^{I} \big( x_{ij} - m_j \big)^2 \big/ I}$$

вычисленным для каждого столбца $\mathbf{x}_j$. Нормирование по строкам (называемое также нормализацией) – это умножение матрицы $\mathbf{X}$ слева на диагональную матрицу $\mathbf{W}$, т.е.

$$\tilde{\mathbf{X}} = \mathbf{XW}$$

При этом размерность $\mathbf{W}$ равна $I×I$, а ее элементы $w_{ii}$ – это обратные значения стандартных отклонений строк $\mathbf{x}_i^\mathrm{t}$.

Комбинация центрирования и нормирования по столбцам

$$\tilde{x}_{ij} = \big( x_{ij} - m_j \big) \big/ d_j$$

называется автошкалированием. Нормирование данных часто применяют для того, чтобы уравнять вклад в модель от различных переменных (например, в гибридном методе ЖХ-МС), учесть неравномерные погрешности, или для того, чтобы обрабатывать совместно разные блоки данных. Нормирование также можно рассматривать как метод, позволяющий стабилизировать вычислительные алгоритмы. В тоже время, к этому преобразованию нужно относиться с большой осторожностью, т.к. оно может сильно исказить результаты качественного анализа. Любое преобразование данных – центрирование, шкалирование, и т.п. – всегда делается сначала на обучающем наборе. По этому набору находятся значения $m_j$ и $d_j$, которые затем применяются и к обучающему, и к проверочному набору