Р 50.1.098-2014. Рекомендации по стандартизации. Статистические методы. Определение и использование линейных функций при калибровке

(утв. и введены в действие Приказом Росстандарта от 24.10.2014 № 1418-ст)

Утверждены и введены в действие Приказом Федерального агентства по техническому регулированию и метрологии от 24 октября 2014 г. № 1418-ст

РЕКОМЕНДАЦИИ ПО СТАНДАРТИЗАЦИИ

СТАТИСТИЧЕСКИЕ МЕТОДЫ

ОПРЕДЕЛЕНИЕ И ИСПОЛЬЗОВАНИЕ ЛИНЕЙНЫХ ФУНКЦИЙ ПРИ КАЛИБРОВКЕ

Statistical methods. Determination and use of straight-line calibration functions

ISO/TS 28037:2010 Determination and use of straight-line calibration functions (IDT)

Р 50.1.098-2014

Группа Т59

ОКС 03.120.30

Дата введения 1 декабря 2015 года

Предисловие

1 ПОДГОТОВЛЕНЫ Открытым акционерным обществом «Научно-исследовательский центр контроля и диагностики технических систем» (АО «НИЦ КД») на основе собственного аутентичного перевода на русский язык международного документа, указанного в пункте 4

2 ВНЕСЕНЫ Техническим комитетом по стандартизации ТК 125 «Применение статистических методов»

3 УТВЕРЖДЕНЫ И ВВЕДЕНЫ В ДЕЙСТВИЕ Приказом Федерального агентства по техническому регулированию и метрологии от 24 октября 2014 г. № 1418-ст

4 Настоящие рекомендации идентичны международному документу ISO/TS 28037:2010 «Определение и использование линейных функций при калибровке» (ISO/TS 28037:2010 «Determination and use of straight-line calibration functions»).

Наименование настоящих рекомендаций изменено относительно наименования указанного международного документа для приведения в соответствие с ГОСТ Р 1.5-2012 (подраздел 3.5).

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

5 ВВЕДЕНЫ ВПЕРВЫЕ

Правила применения настоящих рекомендаций установлены в ГОСТ Р 1.0-2012 (раздел 8). Информация об изменениях к настоящим рекомендациям публикуется в ежегодном (по состоянию на 1 января текущего года) информационном указателе «Национальные стандарты», а официальный текст изменений и поправок — в ежемесячном информационном указателе «Национальные стандарты». В случае пересмотра (замены) или отмены настоящих рекомендаций соответствующее уведомление будет опубликовано в ближайшем выпуске ежемесячного информационного указателя «Национальные стандарты». Соответствующая информация, уведомление и тексты размещаются также в информационной системе общего пользования — на официальном сайте Федерального агентства по техническому регулированию и метрологии в сети Интернет (gost.ru)

Введение

Калибровка во многих случаях является важной частью процедур измерений и часто включает подбор его результатам измерений калибровочной функции, которая наилучшим образом описывает взаимосвязь переменных. В настоящих рекомендациях рассмотрены калибровочные функции, описывающие зависимую переменную Y как линейную функцию независимой переменной X. Параметрами прямой являются параметры A и B. Целью процедуры калибровки является определение оценок a и b параметров A и B для конкретной измерительной системы на основе результатов измерений (xi; yi), i = 1,..., m, выполненных этой измерительной системой. Поскольку результаты измерений обладают неопределенностью, это означает, что оценки a и b также обладают неопределенностью. В настоящих рекомендациях установлен способ определения оценок a и b и соответствующих им неопределенностей по результатам измерений. Использованные в настоящих рекомендациях методы обработки и распространения неопределенности соответствуют Руководству ИСО/МЭК 98-3:2008 «Неопределенность измерения. Часть 3. Руководство по выражению неопределенности в измерении (GUM:1995)».

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

Как только подобрана линейная модель, наилучшим образом соответствующая результатам измерений и требованию состоятельности модели, ее можно использовать для прогноза значения x величины X, соответствующей результату измерения величины Y, полученному с помощью измерительной системы. Калибровочную функцию также можно использовать для оценки неопределенности параметров калибровочной функции и неопределенности прогнозируемого значения x.

Определение и использование линейной калибровочной функции состоят из пяти этапов:

1 Получение информации о неопределенности и ковариации данных результатов измерений. (В рекомендациях приведены соответствующие примеры.)

2 Определение наилучших оценок параметров линейной калибровочной функции.

3 Валидация модели на ее состоятельность и соответствие данным, использование критерия \chi^{2} . (Совместимы ли данные измерений с соответствующими неопределенностями?)

4 Определение стандартной неопределенности и ковариации оценок параметров прямой.

5 Использование калибровочной функции для прогноза, т.е. определение оценки x величины X и ее неопределенности, соответствующих результату y величины Y и ее неопределенности.

Упомянутые этапы показаны в виде схемы на рисунке 1.

Рисунок 1 — Этапы определения и использования линейных калибровочных функций

Приведенные численные методы основаны на [6].

Главной целью настоящих рекомендаций является рассмотрение этапов 2 - 5. Поэтому при использовании настоящих рекомендаций на этапе 1 пользователь должен определить стандартные неопределенности и ковариации, соответствующие результатам измерений величин X и Y. Следует использовать принцип GUM при оценке неопределенности на основе модели измерений, определенной для рассматриваемой области.

В ИСО 11095:1996 (см. [14]) рассмотрены вопросы линейной калибровки с использованием образцов сравнения. Отличия ИСО 11095:1996 от настоящих рекомендаций приведены в таблице 1.

Таблица 1

Отличия ИСО 11095:1996 и настоящих рекомендаций

Характеристика	ИСО 11095:1996	Настоящие рекомендации
Использование специальных образцов сравнения	Да	Более общий случай
Значения X предполагают известными точно	Да	Более общая информация о неопределенности
Все результаты измерений получены независимо	Да	Более общая информация о неопределенности
Соответствие терминологии GUM	Нет	Да
Рассматриваемые типы неопределенности	Два	Пять, включая наиболее общий случай
Только неопределенность, связанная со случайными ошибками	Да	Более общая информация о неопределенности
Проверка сходимости	ANOVA	Критерий \chi^{2}
Неопределенность, соответствующая прогнозируемым значениям	Специальный случай	В соответствии с GUM

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

1. Область применения

В настоящих рекомендациях рассмотрены линейные калибровочные функции, описывающие взаимосвязь переменных X и Y, а именно, функции вида Y = A + BX. Несмотря на то, что многие из положений, установленных в настоящих рекомендациях, применимы и к более общим видам калибровочной функции, в настоящих рекомендациях везде, где это возможно, использована линейная калибровочная функция.

Значения параметров A и B определяют на основе результатов измерений (xi, yi), i = 1,..., m. Рассмотрены различные случаи, касающиеся неопределенности результатов измерений. Не использовано предположение о том, что ошибки yi являются гомоскедастичными (имеют равную дисперсию) и то же для xi, когда ошибки xi не незначительны.

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

Для рассмотренных случаев приведены методы валидации линейной калибровочной функции и оценки неопределенностей и ковариации параметров калибровочной функции.

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

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

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

2. Нормативные ссылки

В настоящих рекомендациях использованы нормативные ссылки на следующие документы:

Руководство ИСО/МЭК 99:2007 Международный словарь по метрологии. Основные и общие понятия и связанные с ними термины (VIM) [ISO/IEC Guide 99:2007 International vocabulary of metrology — Basic and general concepts and associated terms (VIM)]

Руководство ИСО/МЭК 98-3:2008 Неопределенность измерения. Часть 3. Руководство по выражению неопределенности измерений (GUM:1995) [ISO/IEC Guide 98-3:2008, Uncertainty of measurement — Part 3: Guide to the expression of uncertainty in measurement (GUM:1995)]

Руководство ИСО/МЭК 98-3:2008/Дополнение 1:2008 Неопределенность измерения. Часть 3. Руководство по выражению неопределенности измерений (GUM:1995). Дополнение 1. Трансформирование распределений с использованием метода Монте-Карло [ISO/IEC Guide 98-3:2008/Supplement 1:2008, Uncertainty of measurement — Part 3: Guide to the expression of uncertainty in measurement (GUM:1995) — Supplement 1: Propagation of distributions using a Monte Carlo method]

3. Термины и определения

В настоящих рекомендациях применены термины по Руководству ИСО/МЭК 98-3 и Руководству ИСО/МЭК 99, а также следующие термины с соответствующими определениями.

Перечень использованных обозначений приведен в приложении G.

3.1 измеренное значение величины (measured quantity value): Значение, представляющее собой результат измерения величины.

[Руководство ИСО/МЭК 99:2007, 2.10]

3.2 неопределенность измерения (measurement uncertainty): Неотрицательный параметр, характеризующий разброс значений случайной величины, приписываемых ей на основе имеющейся информации об измеряемой величине.

[Руководство ИСО/МЭК 99:2007, 2.26]

3.3 стандартная неопределенность измерения (standard measurement uncertainty): Неопределенность результатов измерений, выраженная в виде стандартного отклонения.

[Руководство ИСО/МЭК 99:2007, 2.30]

3.4 ковариация двух количественных величин (covariance associated with two quantity values): Характеристика взаимозависимости двух количественных величин, которым на основе имеющейся информации, приписывают две измеряемые величины.

3.5 ковариационная матрица, матрица ковариации результатов измерений (measurement covariance matrix, covariance matrix): Матрица размерности № x N, связанная с вектором оценок векторной величины размерности № x 1, содержащая на своей диагонали квадраты стандартной неопределенности соответствующих компонент вектора оценок векторной величины, а в качестве остальных элементов ковариации пар компонентов вектора оценок векторной величины.

Примечание 1 — Ковариационная матрица Ux размерности № x N, соответствующая вектору оценок x векторной величины X, имеет вид:

,

где cov(xi, xi) = u2(xi) — дисперсия (стандартная неопределенность xi);

cov(xi, xj) — ковариация xi и xj, cov(xi, xj) = 0, если элементы Xi и Xj вектора X являются некоррелированными.

Примечание 2 — Ковариацию называют взаимной неопределенностью.

Примечание 3 — Ковариационную матрицу также называют дисперсионно-ковариационной матрицей.

Примечание 4 — Определение соответствует Руководству ИСО/МЭК 98-3:2008/Дополнение 1:2008, определение 3.11 (см. [13]).

3.6 модель измерений (measurement model): Математическая связь всех величин в измерительной задаче.

[Руководство ИСО/МЭК 99:2007, 2.48]

3.7 функциональная модель (functional model): Статистическая модель, включающая ошибки, соответствующие зависимой переменной.

3.8 структурная модель (structural model): Статистическая модель, включающая ошибки, соответствующие независимым и зависимым величинам.

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

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

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

Примечание 3 — Часто под калибровкой понимают только первый этап, указанный в приведенном определении.

[Руководство ИСО/МЭК 99:2007, 2.39]

3.10 распределение вероятностей (probability distribution): Функция (случайной величины), характеризующая вероятность того, что случайная величина принимает данное значение или принадлежит заданному набору значений.

Примечание 1 — Вероятность, соответствующая всему набору значений случайной величины равна 1.

Примечание 2 — Распределение вероятностей называют одномерным, если оно описывает единственную (скалярную) случайную величину, или многомерным, если оно описывает вектор случайных величин. Многомерное распределение вероятностей описывают также как совместное распределение.

Примечание 3 — Распределение вероятностей может иметь форму функции распределения или плотности распределения.

Примечание 4 — Определения и примечание 1 адаптированы по ИСО 3534-1:1993, определение 1.3, и Руководству ИСО/МЭК 98-3:2008, определение C.2.3; примечания 2 и 3 адаптированы по Руководству ИСО/МЭК 98-3:2008/Дополнение 1:2008, определение 3.1 (см. [13]).

3.11 нормальное распределение (normal distribution): Распределение вероятностей непрерывной случайной величины X такое, что соответствующая плотность распределения для - \infty<\xi< + \infty имеет вид:

g_{X}(\xi) = \frac{1}{\sigma \sqrt{2\pi}} exp[- \frac{1}{2} (\frac{\xi - \mu}{\sigma})2] .

Примечание 1 — \mu — математическое ожидание X, \sigma — стандартное отклонение X.

Примечание 2 — Нормальное распределение также называют распределением Гаусса.

Примечание 3 — Определение и примечание 1 адаптированы по ИСО 3534-1:1993, определение 1.37, примечание 2 адаптировано по Руководству ИСО/МЭК 98-3:2008, определение C.2.14.

3.12 t-распределение (t-distribution): Распределение вероятностей непрерывной случайной величины X, плотность распределения которой для - \infty<\xi< + \infty имеет вид:

g_{X}(\xi) = \frac{\text{Г}(\frac{v + 1}{2})}{\sqrt{\pi v} \text{Г}(v/2)} (1 + \frac{\xi^{2}}{v}) - (v + 1)/2 ,

где v — число степеней свободы (положительное целое число);

Г(z) — гамма-функция,

, z > 0.

[Руководство ИСО/МЭК 98-3:2008/Дополнение 1:2008,3.5]

3.13 \chi^{2} - \text{распределение} , распределение хи-квадрат (chi-squared distribution, \chi^{2} distribution): Распределение вероятностей непрерывной случайной величины X, плотность распределения которой для 0 \leq \xi< + \infty имеет вид:

g_{X}(\xi) = \frac{\xi^{(v/2) - 1}}{2^{v/2}\text{Г}(v/2)} exp(- \frac{\xi}{2}) ,

где v — положительное число; Г — гамма-функция.

Примечание - Сумма квадратов v независимых стандартизованных нормальных величин подчиняется \chi^{2} распределению с параметром v; v — число степеней свободы.

3.14 положительно определенная матрица (positive definite matrix): Матрица M размерности n x n, для которой справедливо неравенство zTMz > 0 для всех ненулевых векторов z размерности n x 1.

3.15 положительно полуопределенная матрица (positive semi-definite matrix): Матрица M размерности n x n, для которой справедливо неравенство zTMz >= 0 для всех ненулевых векторов z размерности n x 1.

4. Пояснения к использованным обозначениям

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

4.1 X — независимая величина, Y — зависимая величина, даже если X является неизвестной величиной, а Y — известной, как, например, в разделе 7.

4.2 A и B называют параметрами линейной калибровочной функции Y = A + BX. Их также используют для обозначения (фиктивных) переменных в выражениях, включающих параметры калибровочной функции.

4.3 Величины Xi и Yi используют в качестве (фиктивных) переменных для обозначения координат i-ой точки.

4.4 Константы A* и B* представляют собой (неизвестные) значения A и B, которые определяют линейную калибровочную функцию Y = A* + B*X для рассматриваемой измерительной системы.

4.5 Константы X_{i}^{*} и Y_{i}^{*} представляют собой (неизвестные) координаты i-й точки, полученные измерительной системой и удовлетворяющие уравнению Y_{i}^{*} = A^{*} + B^{*}X_{i}^{*} .

4.6 xi и yi — результаты измерений значений координат i-й точки.

4.7 a и b — оценки параметров калибровочной функции измерительной системы.

4.8 x_{i}^{*} и y_{i}^{*} — оценки координат i-ой точки, удовлетворяющие уравнению y_{i}^{*} = a + bx_{i}^{*} .

4.9 Вектор размерности m x 1

, xT = [x1,..., xm]

и матрица размерности m x n

, .

Для облегчения понимания размерности вектора и матрицы далее всегда такие.

4.10 Т — означает операцию транспонирования.

4.11 Нулевая матрица обозначена 0, а единичный вектор обозначен 1.

4.12 Некоторые символы имеют более одного значения. Необходимые пояснения приведены в тексте.

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

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

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

5. Принципы линейной калибровки

5.1. Общие положения

5.1.1 В данном разделе показано, как соотношение Y = A + BX, описывающее зависимую переменную Y (также называемую «откликом») как функцию независимой переменной X (также называемый «сигналом»), может быть определено по результатам измерений. При калибровке результаты измерений получают с помощью измерительного прибора, которому соответствуют (неизвестные) значения A* и B* параметров калибровочной функции, выполняя измерения на объектах с калиброванными значениями Xi, заданными в стандартных единицах, а результаты измерений Yi фиксируют. Соотношение позволяет определить отклик Y системы для данного объекта с калиброванным значением X. Этот процесс называют предварительной оценкой. Более полезны на практике соотношения, позволяющие преобразовывать измеренное значение y величины Y в оценку x в стандартных единицах X для исследуемого объекта. Этот процесс называют обратной оценкой или прогнозом.

5.1.2 Калибровка измерительной системы должна учитывать неопределенность результатов измерений и соответствующие ковариации. Результатом процедуры калибровки является калибровочная функция, которую используют для прогноза (и при необходимости, предварительной оценки). Результатами калибровки также являются стандартные неопределенности и ковариации, соответствующие оценкам a и b параметров калибровочной функции, которые используют для оценки стандартной неопределенности прогноза (и предварительной оценки).

5.2. Исходные данные для определения калибровочной функции

5.2.1 Данные измерений

Информацией, необходимой для определения уравнения линейной калибровочной функции, являются результаты измерений и соответствующие им стандартные неопределенности и ковариации. В настоящих рекомендациях результаты измерений обозначены (xi, yi), i = 1,..., m, т.е. m пар результатов измерений X и Y. Предполагается, что m >= 2, а значения xi не все равны друг другу.

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

5.2.2 Неопределенности и ковариации

Стандартные неопределенности, соответствующие xi и yi, обозначены u(xi) и u(yi) соответственно. Ковариация xi и xj обозначена cov(xi, xj). Аналогично ковариации yi и yj, xi и yj обозначены cov(yi, yj) и cov(xi, yj), соответственно. В приложении D показано, как могут быть оценены неопределенности и ковариации, соответствующие результатам измерений сигналов и откликов, и приведена интерпретация информации о неопределенности. Полная информация о неопределенности представлена матрицей U размерности 2m x 2m, содержащей дисперсии (квадраты стандартной неопределенности) u2(xi) и u2(yi) и ковариации:

.

Во многих приложениях некоторые или все ковариации принимают равными нулю (см. 5.3).

Примечание - В данных рекомендациях предполагается, что u(xi) и u(yi) различны.

5.3. Определение калибровочной функции

5.3.1 Исходными данными для определения калибровочной функции являются результаты измерений, соответствующие неопределенности и, возможно, ковариации. На основе параметров A и B и исходных данных определяют отклонение i-й точки (xi, yi) от прямой Y = A + BX. Оценки a и b определяют, минимизируя сумму квадратов этих отклонений или более общей меры, если все ковариации отличны от нуля. Как это получить зависит от структуры неопределенности, соответствующей результатам измерений. Структура неопределенности зависит от ответов на следующие вопросы:

i) Неопределенности результатов измерений xi являются несущественными?

ii) Ковариации, соответствующие парам результатов измерений, являются несущественными?

5.3.2 Следующие ситуации, рассмотренные в настоящих рекомендациях, приведены в соответствии с возрастающим порядком сложности в зависимости от ответов на вопросы, приведенные в 5.3.1.

a) неопределенности, соответствующие значениям yi, и все ковариации, соответствующие данным, являются несущественными (раздел 6);

b) неопределенности, соответствующие значениям xi и yi, и все ковариации, соответствующие данным, являются несущественными (раздел 7);

c) имеются неопределенности, соответствующие значениям xi и yi, а ковариации, соответствующие парам (xi, yi), являются несущественными (раздел 8);

d) имеются неопределенности, соответствующие значениям yi, и ковариации, соответствующие yi и yj (i =/ j) (раздел 9);

e) наиболее общий случай, когда имеются неопределенности, соответствующие результатам измерений xi и yi, и ковариации, соответствующие всем парам значений xi, xj, yk и yl (раздел 10).

5.3.3 В каждом случае, перечисленном в 5.3.2, указаны:

a) установленные результаты измерений и структура неопределенности;

b) соответствующая статистическая модель;

c) соответствующая задача метода наименьших квадратов;

d) этапы вычислений;

e) свойства статистической модели;

f) валидация модели (проверка соответствия модели данным);

g) организация выполнения расчетов на компьютере;

h) алгоритм вычислений;

i) один или несколько примеров.

5.4. Числовая обработка

В приложении C приведен подход, использующий ортогональное разложение (факторизацию) матрицы U для наиболее общего случая e) в 5.3.2. Он может быть использован при рассмотрении всех ситуаций. Подход основан на устойчивых методах вычислений. В случаях a) - c) 5.3.2 могут быть использованы элементарные операции, которые могут быть выполнены с помощью электронных таблиц. В случаях d) - e) 5.3.2 необходимо использовать некоторые матричные операции, которые являются прямыми при применении компьютерного языка, допускающего операции с матрицами, но не очень подходят для вычислений с использованием крупноформатных электронных таблиц.

5.5. Неопределенность и ковариация параметров калибровочной функции

5.5.1 Для всех рассмотренных случаев оценки параметров калибровочной функции могут быть представлены (явно или неявно) в виде функции результатов измерений. Принципы GUM [Руководство ИСО/МЭК 98-3:2008] могут быть применены для распространения неопределенности и определения ковариаций, соответствующих результатам измерений с помощью этих функций для получения оценок параметров калибровочной функции. Таким образом, результаты измерений используют для получения оценок a и b параметров калибровочной функции и оценок стандартной неопределенности u(a), u(b) и ковариации cov(a, b), соответствующей этим оценкам. Для случаев a) и d) в 5.3.2 распространение является точным, так как оценки параметров могут быть представлены в виде линейной комбинации входов yi. В других случаях, когда оценки параметров не могут быть так представлены, распространение неопределенности основано на линеаризации оценок параметров. Во многих случаях аппроксимация с помощью линеаризации является достаточно точной.

Примечание - Если распространение неопределенности является приближенным и особенно, если неопределенности являются большими (например, в случаях биологических измерений), может быть использован подход, основанный на распространении распределений. Этот подход [Руководство ИСО/МЭК 98-3:2008/Дополнение 1:2008] использует метод Монте-Карло (не рассматриваемый в настоящих рекомендациях).

5.5.2 Предварительным результатом определения линейной калибровочной функции является вектор оценок параметров размерности 2 x 1 и матрица ковариации Ua размерности 2 x 2

, , (1)

где u(a) и u(b) — стандартные неопределенности оценок a и b соответственно, а cov(a, b) = cov(b, a) — ковариация оценок a и b.

5.6. Валидация модели

5.6.1 При определении оценок a и b параметров линейной калибровочной функции предполагается, что модель Y = A + BX справедлива, а неопределенность, соответствующая результатам измерений, является достоверной мерой отклонения результатов измерений от прямой. После определения a и b, фактическое отклонение точек от наиболее подходящей линейной калибровочной функции может быть определено и сопоставлено с прогнозируемыми отклонениями. При сопоставлении используют совокупную меру отклонений в виде суммы квадратов \chi_{obs}^{2}m взвешенных остатков. В этом случае i-й взвешенный остаток является мерой отклонения i-й точки от прямой. Если ковариация, соответствующая i-й точке (xi, yi), отлична от нуля, может быть использована мера отклонения в более общей форме. Если \chi_{obs}^{2} существенно больше среднего статистических отклонений, есть основание для сомнений в правильности предположения об используемой модели.

5.6.2 Со статистической точки зрения результаты измерений могут быть рассмотрены как реализация случайных величин. Если распределение вероятностей, характеризующее эти случайные величины известно, то можно определить распределение вероятностей для совокупной меры отклонений в 5.6.1. Затем может быть вычислена вероятность того, что \chi_{obs}^{2} (для этого совокупного распределения) превышает заданный квантиль распределения. Однако, поскольку информация об этих величинах часто ограничивается лишь результатами измерений (в виде оценок математического ожидания и дисперсии случайных величин, характеризующихся этими распределениями), этой информации недостаточно для определения распределения вероятностей этой меры. Вместо этого оценку справедливости предположений выполняют предполагая, что распределения этих величин являются нормальными. В этом случае, по крайней мере, для целей валидации, использованным распределением этой меры является распределение \chi_{v}^{2} с v = m — 2 степенями свободы. Соответственно, вероятность того, что \chi_{obs}^{2} превышает заданный квантиль \chi_{v}^{2} , может быть определена (см. 6.3, 7.3, 9.3, 10.3). Обычно используют квантиль уровня 95%.

Примечание 1 — Если \chi_{obs}^{2} превышает квантиль \chi_{v}^{2} уровня 95% это означает, что калибровочная функция не соответствует данным в достаточной мере. В таком случае данные и соответствующая им неопределенность должны быть проверены на наличие ошибок. Может быть использована функция в виде полинома X в степени 2 или более высокой степени или в другой математической форме. Выбор вида калибровочной кривой в настоящих рекомендациях не рассмотрен.

Примечание 2 — Существует возможность, что модель «слишком хороша» в том смысле, что наблюдаемое значение \chi_{obs}^{2} значительно меньше математического ожидания. В этом случае, как правило, неопределенность результатов измерений указана слишком большой. Такой случай в настоящих рекомендациях также не рассмотрен.

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

5.6.4 Если в конкретной ситуации валидация показала несоответствие данных выбранной модели, т.е. \chi_{obs}^{2} превышает квантиль \chi_{v}^{2} уровня 95% (см. 5.6.2), вычисленные стандартные неопределенности u(a) и u(b) и ковариацию cov(a, b) (см. 5.5.2) не следует использовать для расчета неопределенности прогнозируемых значений (см. 5.7).

5.7. Использование калибровочной функции

5.7.1 Калибровочную функцию, как правило, используют для прогноза (обратная оценка), когда по заданному значению Y и соответствующей ему стандартной неопределенности, определяют значение X и соответствующую ему стандартную неопределенность. При определении оценки стандартной неопределенности X используют стандартные неопределенности оценок a и b, а также их ковариацию (см. 11.1).

5.7.2 Иногда при определении оценки X и ее неопределенности, соответствующей значению Y с соответствующей стандартной неопределенностью, необходима предварительная оценка, например, при сопоставлении данных калибровки с набором аналогичных методов (см. 11.2).

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

5.8. Определение наилучшей прямой

5.8.1 Наилучшей прямой, соответствующей исходным данным согласно методу наименьших квадратов, является прямая с коэффициентами a и b (оценками параметров A и B), которые минимизируют сумму

\sum_{i = 1}^{m} (y_{i} - A - Bx_{i})^{2} . (2)

Эти значения удовлетворяют уравнениям, полученным приравниванием к нулю частных производных первого порядка по A и B выражения (2).

5.8.2 Значения оценок a и b могут быть вычислены при выполнении следующих действий:

1. вычисление x_{0} = \frac{1}{m} \sum_{i = 1}^{m} x_{i} и y_{0} = \frac{1}{m} \sum_{i = 1}^{m} y_{i} ;

2. вычисление \tilde{x}_{i} = x_{i} - x_{0} и \tilde{y}_{i} = y_{i} - y_{0} , i = 1,..., m;

3. вычисление b = \frac{\sum_{i = 1}^{m} x_{i}\tilde{y}_{i}}{\sum_{i = 1}^{m} \tilde{x}_{i}^{2}} и a = y0 — bx0.

5.8.3 Значения x0 и y0 таковы, что наилучшая прямая, проведенная через точку (\tilde{x}_{i},\tilde{y}_{i}) , проходит через начало координат и имеет тот же угол наклона, что и наилучшая прямая для исходных данных (xi, yi).

Примечание — Математически, наилучшие параметры определяют, решая систему из двух линейных уравнений, использующих матрицу размерности 2 x 2. Для преобразованных значений эта матрица является диагональной, позволяя легко определить параметры решения. Преобразование данных также позволяет достичь более высокой точности при использовании компьютера (см. [4, страница 33]).

5.8.4 Методы, описанные в разделах 6 - 10, представляют собой расширения вычислений, представленных в 5.8.2 с учетом информации о неопределенности.

6. Модель, учитывающая неопределенность yi

6.1. Общие положения

6.1.1 В данном разделе рассмотрена ситуация 5.3.2 a), а именно, когда имеется следующая информация для i = 1,..., m:

a) результаты измерений (xi, yi);

b) стандартная неопределенность u(yi), соответствующая yi.

В приложении D приведено руководство по получению неопределенности. Все другие неопределенности и ковариации, соответствующие данным, предполагаются несущественными.

6.1.2 Ситуация 5.3.2 a) соответствует статистической модели

yi = A* + B*xi + ei, i = 1,..., m, (3)

где ei — реализации независимых случайных величин с нулевым математическим ожиданием и дисперсиями u2(yi) (см. [9, страница 1]). A* и B* — (неизвестные) значения параметров калибровочной функции для измерительной системы, на которой получены результаты измерений. Эту модель (без учета неопределенности xi) называют функциональной моделью.

6.1.3 В качестве оценок a и b определяют значения, минимизирующие по A и B взвешенную сумму квадратов

\sum_{i = 1}^{m} R_{i}^{2} \equiv \sum_{i = 1}^{m} w_{i}^{2}(y_{i} - A - Bx_{i})^{2} , (4)

wi = 1/u(yi), i = 1,..., m.

Эта задача (4) является методом взвешенных наименьших квадратов. Искомые оценки определяют из уравнений, полученных приравниванием к нулю частных производных первого порядка выражения по A и B (4).

6.2. Оценки параметров калибровки и соответствующие стандартные неопределенности и ковариации

Оценки a и b определяют, выполняя вычисления 1 - 5, затем вычисляют стандартные неопределенности u(a) и u(b) и ковариацию cov(a, b) (вычисление 6):

1. w_{i} = \frac{1}{u(y_{i})} , i = 1,..., m, и F^{2} = \sum_{i = 1}^{m} w_{i}^{2} ;

2. g_{0} = \frac{1}{F^{2}} \sum_{i = 1}^{m} w_{i}^{2}x_{i} и h_{0} = \frac{1}{F^{2}} \sum_{i = 1}^{m} w_{i}^{2}y_{i} ;

3. gi = wi(xi — g0) и hi = wi(yi — h0), i = 1,..., m;

4. G^{2} = \sum_{i = 1}^{m} g_{i}^{2} ;

5. b = \frac{1}{G^{2}} \sum_{i = 1}^{m} g_{i}h_{i} и a = h0 — bg0;

6. u^{2}(a) = \frac{1}{F^{2}} + \frac{g_{0}^{2}}{G^{2}} , u^{2}(b) = \frac{1}{G^{2}} и cov(a,b) = - \frac{g_{0}}{G^{2}} .

Примечание 1 - Вычисления 1 - 5 эквивалентны следующим:

1. w_{i} = \frac{1}{u(y_{i})} , i = 1,..., m, x_{0} = \frac{\sum_{i - 1}^{m} w_{i}^{2}x_{i}}{\sum_{i - 1}^{m} w_{i}^{2}} и y_{0} = \frac{\sum_{i - 1}^{m} w_{i}^{2}y_{i}}{\sum_{i - 1}^{m} w_{i}^{2}} ;

2. \tilde{x}_{i} = x_{i} - x_{0} и \tilde{y}_{i} = y_{i} - y_{0} , i = 1,..., m;

3. b = \frac{\sum_{i - 1}^{m} w_{i}^{2}x_{i}\tilde{y}_{i}}{\sum_{i - 1}^{m} w_{i}^{2}\tilde{x}_{i}^{2}} и a = y0 — bx0.

Примечание 2 — В процессе вычислений 1 - 5 определяют решение системы уравнений (метод наименьших квадратов)

wia + wixib = wiyi, i = 1,..., m.

Примечание 3 — Если все u(yi) идентичны так, что все wi идентичны, то a и b те же, что в 5.8.2.

Примечание 4 — Значения u2(a), u2(b) и cov(a, b) в вычислении 6 получены на основе применения закона распространения неопределенности в соответствии с Руководством ИСО/МЭК 98-3:2008 к a и b в соответствии с вычислениями 1 - 5.

6.2.2 Оценки a и b, полученные в соответствии с 6.1.3, обладают следующими свойствами (см. [15]) для данных yi и модели (3):

1. Оценки a и b являются линейной комбинацией данных yi.

2. Оценки a и b можно рассматривать, как реализации случайных величин с математическими ожиданиями A* и B* соответственно.

3. Матрица ковариации для случайных величин в перечислении 2) имеет элементы u2(a), u2(b) и cov(a, b), вычисленные в соответствии с 6.2.1. Свойство перечисления 1) состоит в том, что оценки a и b получены методом линейной оценки. В соответствии с перечислением 2) оценки являются несмещенными. В соответствии с перечислениями 2) и 3) оценки являются состоятельными, т.е. с увеличением m оценки a и b сходятся к A* и B* соответственно.

Метод, установленный в 6.1.3, обладает следующим оптимальным свойством для данных yi и модели (3):

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

Свойство перечисления 4) можно интерпретировать следующим образом. Для констант c и d, стандартная неопределенность u(ca + db) , соответствующая линейной комбинации оценок a и b , полученных любым несмещенным методом линейной оценки, является, по крайней мере, столь же большой, как u(ca + db). Свойства перечислений 1) - 4) оправдывают использование методов наименьших квадратов для данных, совместимых с моделью (3). Необходимо заметить, что все утверждения относятся только к математическим ожиданиям и дисперсиям ei, соответствующие распределения далее не определены. Если сделано дополнительное предположение о том, что ei являются реализацией нормально распределенных случайных величин, то могут быть сделаны утверждения о следующих свойствах, связанных с методом взвешенных наименьших квадратов:

5. Случайные величины в перечислении 2) характеризуются двумерным нормальным распределением со средними A* и B* и матрицей ковариаций с элементами u2(a), u2(b) и cov(a, b).

6. Оценки a и b являются оценками максимального правдоподобия, соответствующими наиболее вероятным значениям A и B, которые, возможно, могут быть определены по наблюдениям (результатам измерений yi).

7. С позиции Байесовского анализа распределение знаний об A и B с учетом наблюдаемых результатов измерений yi является двумерным нормальным распределением, со средними a и b и ковариационной матрицей с элементами u2(a), u2(b) и cov(a, b).

6.3. Валидация модели

Если m > 2, соответствие модели исходным данным может быть проверено с использованием взвешенных остатков ri (продолжение 6.2.1). Для этого выполняют следующие действия:

8. Формирование ri = wi(yi — a - bxi), i = 1,..., m.

9. Вычисление наблюдаемого значения хи-квадрат \chi_{obs}^{2} = \sum_{i = 1}^{m} r_{i}^{2} и определение числа степеней свободы v = m — 2.

10. Сопоставление \chi_{obs}^{2} с квантилем \chi_{v}^{2} распределения уровня 95%. Если значение x_{obs}^{2} больше этого квантиля, линейную модель отклоняют.

Примечание — Критерий \chi^{2} основан на предположении, что ei в модели (3) являются реализацией независимых нормальных случайных величин.

6.4. Организация вычислений

Вычисления в 6.2.1 и 6.3 могут быть выполнены в одной или двух таблицах при использовании электронных таблиц, в соответствии с таблицами 2 и 3, которые могут быть объединены в одну таблицу.

Таблица 2

Данные для определения линейной калибровочной функции методом взвешенных наименьших квадратов

xi	yi	u(yi)
x1	y1	u(y1)
x2	y2	u(y2)
\vdots	\vdots	\vdots
xm	ym	u(ym)

Таблица 3

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

							6.2.1, вычисление 5
				6.2.1, вычисление 2, 3			6.3, вычисление 7	6.3, вычисление 8
				g0	h0			a
w1	w_{1}^{2}	w_{1}^{2}x_{1}	w_{1}^{2}y_{1}	g1	h1	g_{1}^{2}	g1h1	r1	r_{1}^{2}
w2	w_{2}^{2}	w_{2}^{2}x_{2}	w_{2}^{2}y_{2}	g2	h2	g_{2}^{2}	g2h2	r2	r_{2}^{2}
\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots
wm	w_{m}^{2}	w_{m}^{2}x_{m}	w_{m}^{2}y_{m}	gm	hm	g_{m}^{2}	gmhm	rm	r_{m}^{2}
	F^{2} = \sum w_{i}^{2}	\sum w_{i}^{2}x_{i}	\sum w_{i}^{2}y_{i}			G^{2} = \sum g_{i}^{2}	\sum g_{i}h_{i}	b	\chi_{obs}^{2} = \sum r_{i}^{2}

Пример - (равные веса). В таблице 4 приведено шесть значений и соответствующие им значения стандартной неопределенности. Результаты измерений xi являются точными, а стандартная неопределенность yi равна u(yi) = 0,5. Поэтому wi = 1/u(yi) = 2,0, i = 1,..., 6.

Таблица 4

Данные, представляющие результаты шести измерений с равными весами

xi	yi	u(yi)
1,0	3,3	0,5
2,0	5,6	0,5
3,0	7,1	0,5
4,0	9,3	0,5
5,0	10,7	0,5
6,0	12,1	0,5

Результаты вычислений приведены в таблице 5. В соответствии с таблицей 5 g0 = 84,000/24,000 = 3,500, h0 = 192,400/24,000 = 8,017, b = 123,000/70,000 = 1,757 и a = 8,017 — (1,757)(3,500) = 1,867.

Таблица 5

Вычисления на основе данных таблицы 4

wi	w_{i}^{2}	w_{i}^{2}x_{i}	w_{i}^{2}y_{i}	gi	hi	g_{i}^{2}	gihi	ri	r_{i}^{2}
				3,500	8,017			a = 1,867
2,000	4,000	4,000	13,200	-5,000	-9,433	25,000	47,167	-0,648	0,419
2,000	4,000	8,000	22,400	-3,000	-4,833	9,000	14,500	0,438	0,192
2,000	4,000	12,000	28,400	-1,000	-1,833	1,000	1,833	-0,076	0,006
2,000	4,000	16,000	37,200	1,000	2,567	1,000	2,567	0,810	0,655
2,000	4,000	20,000	42,800	3,000	5,367	9,000	16,100	0,095	0,009
2,000	4,000	24,000	48,400	5,000	8,167	25,000	40,833	-0,619	0,383
	24,000	84,000	192,400			70,000	123,000	b = 1,757	1,665

Стандартная неопределенность и ковариация, соответствующие параметрам прямой, могут быть вычислены на основе формулы, приведенной в 6.2.1 и данных таблицы 5:

u2(a) = 1/24,000 + (3,500)2/70,000, так что u(a) = 0,465;

u2(b) = 1/70,000, так что u(b) = 0,120;

cov(a, b) = -3,500/70,000 = -0,050.

Наблюдаемое значение \chi_{obs}^{2} = 1,665 с v = 4. Так как \chi_{obs}^{2} не превышает квантиль \chi_{v}^{2} уровня 95%, а именно, 9,488, можно считать, что данные соответствуют модели.

Данные и полученная линейная калибровочная функция показаны на рисунке 2. Стандартная неопределенность yi показана вертикальными отрезками, охватывающими yi, конечные точки которых равны соответственно (yi — u(yi)) и (yi + u(yi)). Взвешенные остатки показаны на рисунке 3.

Рисунок 2 — Данные таблицы 4 и линейная калибровочная функция, полученная в таблице 5

Рисунок 3 — Взвешенные остатки ri, полученные в таблице 5

Пример - (неравные веса). В таблице 6 приведено шесть значений и соответствующие им стандартные неопределенности. Значения xi измерены точно. Значения yi получены с помощью двух настроек прибора так, что для больших значений X yi являются менее точными.

Таблица 6

Данные, представляющие шесть результатов измерений (неравные веса)

xi	yi	u(yi)
1,0	3,2	0,5
2,0	4,3	0,5
3,0	7,6	0,5
4,0	8,6	1,0
5,0	11,7	1,0
6,0	12,8	1,0

Результаты вычислений приведены в таблице 7. В соответствии с таблицей 7 g0 = 39,000/15,000 = 2,600, h0 = 93,500/15,000 = 6,233, b = 65,000/31,600 = 2,057 и a = 6,233 — (2,057)(2,600) = 0,885.

Таблица 7

Вычисления для данных таблицы 6

wi	w_{i}^{2}	w_{i}^{2}x_{i}	w_{i}^{2}y_{i}	gi	hi	g_{i}^{2}	gihi	ri	r_{i}^{2}
				2,600	6,233			a = 0,885
2,000	4,000	4,000	12,800	-3,200	-6,067	10,240	19,413	0,516	0,266
2,000	4,000	8,000	17,200	-1,200	-3,867	1,440	4,640	-1,398	1,955
2,000	4,000	12,000	30,400	0,800	2,733	0,640	2,187	1,088	1,183
1,000	1,000	4,000	8,600	1,400	2,367	1,960	3,313	-0,513	0,263
1,000	1,000	5,000	11,700	2,400	5,467	5,760	13,120	0,530	0,281
1,000	1,000	6,000	12,800	3,400	6,567	11,560	22,327	-0,427	0,182
	15,000	39,000	93,500			31,600	65,000	b = 2,057	4,131

Стандартная неопределенность и ковариация, соответствующие параметрам прямой, могут быть вычислены, используя формулу, приведенную в 6.2.1 и данные таблицы 7;

u2(a) = 1/15,000 + (2,600)2/31,600, так чтобы u(a) = 0,530;

u2(b) = 1/31,600, так чтобы u(b) = 0,178;

cov(a, b) = -2,600/31,600 = -0,082.

Наблюдаемое значение \chi_{obs}^{2} = 4,131 с v = 4 степенями свободы. Так как \chi_{obs}^{2} не превышает квантиль \chi_{v}^{2} уровня 95%, а именно 9,488, можно считать, что данные соответствуют линейной модели.

Данные и полученная линейная калибровочная функция показаны на рисунке 4. Взвешенные остатки показаны на рисунке 5.

Рисунок 4 — Данные таблицы 6 и линейная калибровочная функция, полученная в таблице 7

Рисунок 5 — Взвешенные остатки ri, полученные в таблице 7

7. Модель, учитывающая неопределенности xi и yi

7.1. Общие положения

7.1.1 В данном разделе рассмотрена ситуация 5.3.2 b), когда имеется следующая информация для i = 1,..., m:

a) результаты измерений (xi, yi);

b) стандартная неопределенность u(xi), соответствующая xi;

c) стандартная неопределенность u(yi), соответствующая yi.

В приложении D приведено руководство по определению неопределенности. Все ковариации, соответствующие данным, считаются несущественными.

7.1.2 Ситуации 5.3.2 b) соответствует статистическая модель

x_{i} = X_{i}^{*} + d_{i} , y_{i} = Y_{i}^{*} + e_{i} , Y_{i}^{*} = A^{*} + B^{*}X_{i}^{*} , i = 1,..., m, (5)

где di и ei — реализации независимых случайных величин с нулевым математическим ожиданием и дисперсиями u2(xi) и u2(yi), соответственно. Эту модель называют структурной моделью. В модели (xi, yi) представляют измеренные координаты точек ( X_{i}^{*} , Y_{i}^{*} ), лежащих на линии Y = A* + B*X.

7.1.3 Поскольку xi (в дополнение к yi — см. раздел 6) соответствуют неопределенности, это необходимо учитывать при определении линейной калибровочной функции. Задача определения a и b в этом случае является одной из задач взвешенной ортогональной регрессии (см. [3]) или обобщенной регрессии (GDR) (см. [2])). В статистической литературе ее называют моделью ошибок в переменных (см. [7], [9], [17]). Оценки a и b обеспечивают минимум по A, B, и Xi, i = 1,..., m с весами \upsilon_{i} = 1/u(x_{i}) и wi = 1/u(yi) сумме квадратов

\sum_{i = 1}^{m} [\upsilon_{i}^{2}(x_{i} - X_{i})^{2} + w_{i}^{2}(y_{i} - A - BX_{i})^{2}] . (6)

Каждая оценка x_{i}^{*} вместе с a и b определяет оценку (x_{i}^{*},y_{i}^{*}) , y_{i}^{*} = a + bx_{i}^{*} для (X_{i}^{*},Y_{i}^{*}) модели (5).

7.1.4 Данные A и B и значения x_{i}^{*} , минимизирующие сумму квадратов (6) относительно Xi, удовлетворяют соотношению

x_{i}^{*} = x_{i}^{*}(A,B) = [u^{2}(y_{i})x_{i} + (y_{i} - A)Bu^{2}(x_{i})]T_{i} , T_{i} = \frac{1}{u^{2}(y_{i}) + B^{2}u^{2}(x_{i})} . (7)

Используя выражение (7) и выполняя замену Xi в выражении (6) на x_{i}^{*}(A,B) , можно записать задачу оптимизации для параметров A и B.

\sum_{i = 1}^{m} {_{iiiiii}^{\upsilon^{2}[x - x^{*}(A,B)]^{2} + w^{2}[y - y^{*}(A,B)]^{2}}} , y_{i}^{*}(A,B) = A + Bx_{i}^{*}(A,B) . (8)

7.1.5 Если

{}_{ii}^{R = R(A,B) =}{_{iiii}^{- B[x - x^{2}(A,B)] + [y - y^{*}(A,B)]}}_{i}^{T^{12}} , (9)

то сумма квадратов (8) эквивалентна

\sum_{i = 1}^{m} R_{i}^{2} .

Для Ri существует следующая геометрическая интерпретация. Вектор, перпендикулярный к прямой Y = A + BX, имеет вид (-B, 1)T/(1 + B2)1/2, Ri — весовой коэффициент соответствующего компонента (x_{i} - x_{i}^{2}(A,B),y_{i} - y_{i}^{*}(A,B))^{T} в направлении этого вектора (с учетом его значения).

Примечание 1 — В обычных методах наименьших квадратов (см. 5.8) и взвешенных наименьших квадратов (см. раздел 6), расстояние до линии измеряют по вертикали, т.е. в направлении оси Y, отражая тот факт, что отклонение (xi, yi) от линии может быть вычислено с помощью ошибки ei, связанной с yi, так как xi предполагается точно известным. Метод взвешенной ортогональной регрессии применяют в случае, когда существует неопределенность xi.

Примечание 2 - Выражение (7) получено приравниванием к нулю частных производных первого порядка выражения (6) по A, B и Xi.

Примечание 3 — Если ui(xi) = 0, то в выражениях (7) x_{i}^{*}(A,B) = x_{i} (т.е. y_{i}^{*}(A,B) = A + Bx_{i} ) и T_{i} = 1/u^{2}(y_{i}) = w_{i}^{2} . Следовательно, Ri в выражении (9) принимает вид Ri = wi(yi — A - Bxi). Таким образом, если u(xi) = 0, Ri оценивают аналогично (4) в 6.1.3.

Примечание 4 — Если u(xi) = u(yi) = ui, то x_{i}^{*}(A,B) определяет точку на линии Y = A + BX самую близкую к (xi, yi)

T_{i} = \frac{1}{u_{i}^{2}} \frac{1}{1 + B^{2}} , R_{i} = \frac{1}{u_{i}} \frac{1}{(1 + B^{2})^{12}} {_{iiii}^{- B[x - x^{*}(A,B)] + [y - y^{*}(A,B)]}} .

Так как (-B, 1)T/(1 + B2)1/2 — вектор, перпендикулярный к линии, Ri — взвешенное расстояние от точки (xi, yi) до линии Y = A + BX.

7.1.6 В 7.1.3 A, B и Xi, i = 1,..., m использованы как переменные при минимизации. В 7.2.1 приведены вычисления по этой минимизации в процессе двухэтапной итерации (см. [2]):

1. по приближениям a и b определяют соответствующее оптимальное x_{i}^{*} ;

2. на основе x_{i}^{*} определяют новые приближения a и b, которые уменьшают сумму квадратов (6).

Примечание - В рекомендациях не использованы различные обозначения x_{i}^{*} для итераций и окончательного результата.

7.2. Оценки параметров, соответствующие стандартные неопределенности и ковариации

7.2.1 Оценки a и b определяют, выполняя вычисления 1 - 6, используя схему итерации, описанную в 7.1.6; стандартные неопределенности u(a) и u(b) и ковариацию cov(a, b) определяют, выполняя вычисления 7 (см. приложение B):

1. определение начального приближения \tilde{a} и \tilde{b} оценок a и b, например, определение методом взвешенных наименьших квадратов наилучшей линии (см. 6.2.1 вычисления 1 - 5), игнорируя наличие неопределенности xi;

2. t_{i} = \frac{1}{u^{2}(y_{i}) + \tilde{b}^{2}u^{2}(x_{i})} , x_{i}^{*} = [x_{i}u^{2}(y_{i}) + (y_{i} - \tilde{a})\tilde{b}u^{2}(x_{i})]t_{i} и z_{i} = y_{i} - \tilde{a} - \tilde{b}x_{i} , i = 1,..., m;

3. f_{i} = t_{i}^{12} , g_{i} = f_{i}x_{i}^{*} и hi = fizi, i = 1,..., m;

4. определение решения методом (невзвешенных) наименьших квадратов \delta a и \delta b для системы уравнений (\delta A)f_{i} + (\delta B)g_{i} = h_{i} , i = 1,..., m:

a) F^{2} = \sum_{i = 1}^{m} f_{i}^{2} ;

b) g_{0} = \frac{1}{F^{2}} \sum_{i = 1}^{m} f_{i}g_{i} и h_{0} = \frac{1}{F^{2}} \sum_{i = 1}^{m} f_{i}h_{i} ;

c) \tilde{g}_{i} = g_{i} - g_{0}f_{i} и h_{i} = h_{i} - h_{0}f_{i} , i = 1,..., m;

d) \tilde{G}^{2} = \sum_{i = 1}^{m} \tilde{g}_{i}^{2} ;

e) \delta b = \frac{1}{\tilde{G}^{2}} \sum_{i = 1}^{m} \tilde{g}_{i}h_{i} и \delta a = h_{0} - (\delta b)g_{0} ;

5. обновление параметров и остатков: \tilde{a}: = \tilde{a} + \delta a , \tilde{b}: = \tilde{b} + \delta b , r_{i} = h_{i} - (\delta b)\tilde{g}_{i} , i = 1,..., m;

6. повторение вычислений 2 - 5, до тех пор, пока не будет достигнута необходимая сходимость. Присвоение a = \tilde{a} , b = \tilde{b} ;

7. вычисление u^{2}(a) = \frac{1}{F^{2}} + \frac{g_{0}^{2}}{G^{2}} , u^{2}(b) = \frac{1}{\tilde{G}^{2}} и cov(a,b) = - \frac{g_{0}}{\tilde{G}^{2}} , где g0, h0 и т.д. — значения, полученные при выполнении вычислений 4.

Примечание 1 - Вычисления 4 аналогичны вычислениям 1 - 5 в 6.2.1.

Примечание 2 — При выполнении вычислений 2 значению x_{i}^{*} соответствует точка (x_{i}^{*},a + b,x_{i}^{*}) текущего приближения наилучшей линейной калибровочной функции, наиболее близкой к точке результатов измерений (xi, yi) (с учетом взвешенного расстояния).

Примечание 3 — При выполнении вычислений 3 значение hi представляет значение обобщенного расстояния Ri в выражении (9) от i-й точки до текущей оценки линейной калибровочной функции. Алгоритм минимизирует сумму квадратов таких расстояний.

Примечание 4 — При выполнении вычислений 4 значения \delta a и \delta b уменьшаются в одно и тоже количество раз от итерации до итерации. Коэффициент уменьшения зависит в значительной степени от неопределенности, соответствующей данным: чем меньше эта неопределенность, тем больше уменьшение. Итерации могут быть закончены, когда величины \delta a и \delta b становятся несущественными.

Примечание 5 — Остатки, вычисленные в соответствии с 5, связаны с решением системы уравнений, решаемой ранее, при выполнении вычислений 4. Конвергенция ri, определенная при выполнении вычислений 5, совпадает с hi, значение которой определено при выполнении вычислений 3.

Примечание 6 — Строго говоря, остатки (см. вычисление 5) необходимы только на последующей итерации. Однако, в формате таблицы (таблица 9 в 7.4) остатки вычисляют на каждой итерации.

Примечание 7 — Значения u2(a), u2(b) и cov(a, b) (вычисление 7) получены с применением закона распространения неопределенности по Руководству ИСО/МЭК 98-3:2008 к a и b в соответствии с вычислениями 1 - 6.

7.2.2 Несмотря на то, что свойства оценки, полученной методом взвешенных наименьших квадратов можно определить в 6.2.2, оценки a и b соответствуют минимуму суммы квадратов (6) и нелинейно зависят от данных xi и yi. Это означает, что соответствующие свойства оценок взвешенной ортогональной регрессии не могут быть установлены прямо. Оценки a и b, определенные в 7.1.3, обладают следующими свойствами для данных xi и yi, соответствующих модели (5):

1. оценки a и b являются нелинейными функциями данных xi и yi.

2. оценки a и b можно рассматривать как реализацию случайных величин, математические ожидания которых составляют приблизительно A* и B*, соответственно.

3. элементы ковариационной матрицы для случайных величин в 2) близки к u2(a), u2(b) и cov(a, b), вычисленным в 7.2.1.

Приближения в 2) и 3) являются более точными для данных с меньшей неопределенностью. Однако метод оценки обладает следующими свойствами:

4. для данных, удовлетворяющих модели (5) с увеличением m, оценки a и b сходятся к A* и B*, соответственно (см. [16]).

Метод взвешенных наименьших квадратов недооценивает угловой коэффициент (см. [5]) для данных, соответствующих модели (5).

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

5. случайные переменные в 2 подчиняются приближенно двумерному нормальному распределению со средними A* и B* и ковариационной матрицей с элементами u2(a), u2(b) и cov(a, b);

6. оценки a и b являются оценками максимального правдоподобия, соответствующими наиболее вероятным значениям A и B, которым соответствуют наблюдаемые результаты измерений xi и yi;

7. в соответствии с Байесовским анализом распределение, характеризующее знания об A и B с учетом наблюдаемых результатов измерений xi и yi, является приближенно двумерным нормальным распределением со средними a и b и ковариационной матрицей с элементами u2(a), u2(b) и cov(a, b).

7.3. Валидация модели

Если m > 2, соответствие данных модели может быть частично проверено с использованием взвешенных остатков ri, (вычисление 5 в 7.2.1) при их сходимости в процессе итераций (продолжение 7.2.1):

8. определение наблюдаемого значения \chi_{obs}^{2} = \sum_{i = 1}^{m} r_{i}^{2} с числом степеней свободы v = m — 2;

9. сопоставление \chi_{obs}^{2} с квантилем распределения \chi_{v}^{2} уровня 95%. Если \chi_{obs}^{2} превышает этот квантиль, модель не соответствует исходным данным.

Примечание — Тест \chi^{2} основан на предположении, что di и ei в модели (5) представляют собой реализации независимых случайных величин, подчиняющихся нормальному распределению на первом этапе итерации.

7.4. Организация вычислений

Вычисления в 7.2.1 и 7.3 могут быть выполнены в двух последовательно дополняемых таблицах, подходящих для использования электронных таблиц. В первой таблице (таблица 8) даны приближения \tilde{a} и \tilde{b} (см. 7.2.1 вычисление 1), вычисление fi, gi и hi (см. 7.2.1 вычисление 3). Во второй таблице (таблица 9) использованы значения fi, gi и hi для вычисления поправок \delta a и \delta b (см. 7.2.1 вычисление 4).

Таблица 8

Вычисления для определения параметров a и b линейной калибровочной функции по данным приближениям \tilde{a} и \tilde{b}

				7.2.1, шаги 1 - 5
					7.2.1, шаг 2	7.2.1, шаг 2	7.2.1, шаг 3
				\tilde{a}	\tilde{b}
x1	u(x1)	y1	u(y1)	t1	x_{1}^{*}	z1	f1	g1	h1
x2	u(x2)	y2	u(y2)	t2	x_{2}^{*}	z2	f1	g2	h2
\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots
xm	u(xm)	ym	u(ym)	tm	x_{m}^{*}	zm	fm	gm	hm

Таблица 9

Организация вычислений для определения поправок \delta a и \delta b для GDR линейной калибровочной функции

			7.2.1, шаги 4 a), 4 b)			7.2.1, шаги 4 e), 5		7.3, шаг 8
			g0	h0			\delta a
f_{1}^{2}	f1g1	f1h1	\tilde{g}_{1}	h_{1}	\tilde{g}_{1}^{2}	\tilde{g}_{1}h_{1}	r1	r_{1}^{2}
f_{1}^{2}	f2g2	f2h2	\tilde{g}_{2}	\tilde{h}_{2}	\tilde{g}_{2}^{2}	\tilde{g}_{2}\tilde{h}_{2}	r2	r_{2}^{2}
\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots
f_{m}^{2}	fmgm	fmhm	\tilde{g}_{m}	\tilde{h}_{m}	\tilde{g}_{m}^{2}	\tilde{g}_{m}\tilde{h}_{m}	rm	r_{m}^{2}
F^{2} = \sum f_{i}^{2}	\sum f_{i}g_{i}	\sum f_{i}h_{i}			\tilde{G}^{2} = \sum \tilde{g}_{i}^{2}	\sum \tilde{g}_{i}h_{i}	\delta b	\sum r_{i}^{2}

Пример - В таблице 10 приведены шесть результатов измерений и соответствующие им стандартные неопределенности.

Таблица 10

Шесть результатов измерений с соответствующими неопределенностями

xi	u(xi)	yi	u(yi)
1,2	0,2	3,4	0,2
1,9	0,2	4,4	0,2
2,9	0,2	7,2	0,2
4,0	0,2	8,5	0,4
4,7	0,2	10,8	0,4
5,9	0,2	13,5	0,4

Для определения начальных приближений \tilde{a} и \tilde{b} (7.2.1 вычисления 1) используют метод взвешенных наименьших квадратов и определяют параметры линейной калибровочной функции. После схемы, описанной в 6.2, получают таблицы 11 и 12.

Таблица 11

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

xi	yi	u(yi)
1,2	3,4	0,2
1,9	4,4	0,2
2,9	7,2	0,2
4,0	8,5	0,4
4,7	10,8	0,4
5,9	13,5	0,4

Таблица 12

Вычисление начальных аппроксимаций \tilde{a} и \tilde{b} на основе данных таблицы 11

wi	w_{i}^{2}	w_{i}^{2}x_{i}	w_{i}^{2}y_{i}	gi	hi	g_{i}^{2}	gihi	ri	r_{i}^{2}
				2,5733	6,1867			a = 0,6583
5,0000	25,0000	30,0000	85,0000	-6,8667	-13,9333	47,1511	95,6756	0,8186	0,6701
5,0000	25,0000	47,5000	110,0000	-3,3667	-8,9333	11,3344	30,0756	-1,7006	2,8920
5,0000	25,0000	72,5000	180,0000	1,6333	5,0667	2,6678	8,2756	1,5577	2,4264
2,5000	6,2500	25,0000	53,1250	3,5667	5,7833	12,7211	20,6272	-1,8791	3,5310
2,5000	6,2500	29,3750	67,5000	5,3167	11,5333	28,2669	61,3189	0,1113	0,0124
2,5000	6,2500	36,8750	84,3750	8,3167	18,2833	69,1669	152,0564	0,4163	0,1733
	93,7500	241,2500	580,0000			171,3083	368,0292	b = 2,1483	9,7052

Начальные приближения \tilde{a} = 0,6583 и \tilde{b} = 2,1483 . На основе приближений вычисляют fi, gi и hi (таблица 13). Затем вычисляют поправки (таблица 14) \delta a = - 0,0784 и \delta b = 0,0111 (7.2.1 вычисления 4). В конце итерации приближения \tilde{a} и \tilde{b} обновляют (7.2.1 вычисления 5):

\tilde{a}: = \tilde{a} + \delta a = 0,6583 - 0,0784 = 0,5799 ;

\tilde{b}: = \tilde{b} + \delta b = 2,1483 + 0,0111 = 2,1594 ;

Таблица 13

Первая итерация при определении fi, gi и h на основе \tilde{a} и \tilde{b}

xi	u(xi)	yi	u(yi)	ti	x_{i}^{*}	zi	fi	gi	hi
				0,6583	2,1483
1,2000	0,2000	3,4000	0,2000	4,4522	1,2626	0,1637	2,1100	2,6642	0,3455
1,9000	0,2000	4,4000	0,2000	4,4522	1,7699	-0,3401	2,1100	3,7345	-0,7176
2,9000	0,2000	7,2000	0,2000	4,4522	3,0192	0,3116	2,1100	6,3706	0,6575
4,0000	0,2000	8,5000	0,4000	2,9019	3,8126	-0,7515	1,7035	6,4947	-1,2802
4,7000	0,2000	10,8000	0,4000	2,9019	4,7111	0,0447	1,7035	8,0253	0,0761
5,9000	0,2000	13,5000	0,4000	2,9019	5,9416	0,1667	1,7035	10,1214	0,2840

Таблица 14

Первая итерация при определении поправок \delta a и \delta b на основе fi, gi и hi

f_{i}^{2}	fi	fihi	\tilde{g}_{i}	h_{i}	\tilde{g}_{i}^{2}	\tilde{g}_{i}h_{i}	ri	r_{i}^{2}
			3,1239	-0,0437			\delta a = - 0,0784
4,4522	5,6216	0,7290	-3,9273	0,4378	15,4236	-1,7193	0,4814	0,2318
4,4522	7,8799	-1,5141	-2,8570	-0,6253	8,1623	1,7864	-0,5935	0,3523
4,4522	13,4422	1,3874	-0,2209	0,7498	0,0488	-0,1656	0,7523	0,5659
2,9019	11,0636	-2,1807	1,1732	-1,2057	1,3764	-1,4145	-1,2187	1,4852
2,9019	13,6710	0,1297	2,7038	0,1506	7,3108	0,4073	0,1206	0,0145
2,9019	17,2416	0,4838	4,7999	0,3585	23,0387	1,7208	0,3052	0,0931
22,0622	68,9199	-0,9648			55,3606	0,6152	\delta b = 0,0111	2,7429

По этим обновленным значениям \tilde{a} и \tilde{b} снова выполняют вычисления (таблицы 15 и 16) и определяют поправки \delta a = - 0,0010 и \delta b = 0,0002 . Процесс повторяют в третий раз (таблицы 17 и 18). В этом случае значения поправок менее 0,0005 можно считать незначительными, а итерацию с оценками параметров a = 0,5788 и b = 2,1597 заключительной.

Таблица 15

Вторая итерация (таблица, аналогичная таблице 13)

xi	u(xi)	yi	u(yi)	ti	x_{i}^{*}	zi	fi	gi	hi
				0,5799	2,1594
1,2000	0,2000	3,4000	0,2000	4,4146	1,2873	0,2289	2,1011	2,7047	0,4808
1,9000	0,2000	4,4000	0,2000	4,4146	1,7922	-0,2827	2,1011	3,7655	-0,5941
2,9000	0,2000	7,2000	0,2000	4,4146	3,0365	0,3579	2,1011	6,3799	0,7519
4,0000	0,2000	8,5000	0,4000	2,8858	3,8212	-0,7175	1,6988	6,4913	-1,2189
4,7000	0,2000	10,8000	0,4000	2,8858	4,7177	0,0709	1,6988	8,0142	0,1205
5,9000	0,2000	13,5000	0,4000	2,8858	5,9448	0,1796	1,6988	10,0988	0,3051

Таблица 16

Вторая итерация (таблица, аналогичная таблице 14)

f_{i}^{2}	fi	fihi	\tilde{g}_{i}	h_{i}	\tilde{g}_{i}^{2}	\tilde{g}_{i}h_{i}	ri	r_{i}^{2}
			3,1412	-0,0003			\delta a = - 0,0010
4,4146	5,6827	1,0103	-3,8953	0,4814	15,1734	-1,8751	0,4823	0,2326
4,4146	7,9117	-1,2482	-2,8344	-0,5935	8,0339	1,6822	-0,5928	0,3514
4,4146	13,4047	1,5798	-0,2201	0,7524	0,0484	-0,1656	0,7525	0,5662
2,8858	11,0271	-2,0706	1,1551	-1,2184	1,3342	-1,4074	-1,2187	1,4852
2,8858	13,6143	0,2046	2,6781	0,1209	7,1720	0,3238	0,1203	0,0145
2,8858	17,1555	0,5183	4,7626	0,3056	22,6824	1,4553	0,3044	0,0927
21,9012	68,7961	-0,0057			54,4443	0,0132	\delta b = 0,0002	2,7427

Таблица 17

Третья итерация (таблица, аналогичная таблице 13)

xi	u(xi)	yi	u(yi)	ti	x_{i}^{*}	zi	fi	gi	hi
				0,5788	2,1597
1,2000	0,2000	3,4000	0,2000	4,4138	1,2875	0,2296	2,1009	2,7050	0,4823
1,9000	0,2000	4,4000	0,2000	4,4138	1,7924	-0,2822	2,1009	3,7657	-0,5928
2,9000	0,2000	7,2000	0,2000	4,4138	3,0366	0,3582	2,1009	6,3795	0,7525
4,0000	0,2000	8,5000	0,4000	2,8855	3,8212	-0,7174	1,6987	6,4909	-1,2187
4,7000	0,2000	10,8000	0,4000	2,8855	4,7176	0,0708	1,6987	8,0137	0,1203
5,9000	0,2000	13,5000	0,4000	2,8855	5,9447	0,1792	1,6987	10,0980	0,3044

Таблица 18

Вторая итерация (таблица, аналогичная таблице 14)

f_{i}^{2}	figi	fihi	\tilde{g}_{i}	h_{i}	\tilde{g}_{i}^{2}	\tilde{g}_{i}h_{i}	ri	r_{i}^{2}
			3,1414	0,0000			\delta a = 0,0000
4,4138	5,6829	1,0133	-3,8947	0,4823	15,1685	-1,8785	0,4823	0,2327
4,4138	7,9113	-1,2454	-2,8340	-0,5928	8,0315	1,6800	-0,5928	0,3514
4,4138	13,4027	1,5809	-0,2202	0,7525	0,0485	-0,1657	0,7525	0,5662
2,8855	11,0258	-2,0702	1,1548	-1,2187	1,3335	-1,4073	-1,2187	1,4852
2,8855	13,6126	0,2043	2,6776	0,1203	7,1695	0,3220	0,1203	0,0145
2,8855	17,1531	0,5171	4,7619	0,3044	22,6756	1,4496	0,3044	0,0927
21,8977	68,7884	0,0000			54,4271	0,0001	\delta b = 0,0000	2,7427

Стандартная неопределенность и ковариация (7.2.1 вычисления 7), соответствующие этим параметрам, могут также быть оценены по данным таблицы 18:

u2(a) = 1/21,8977 + (3,1414)2/54,4271 так, что u(a) = 0,4764;

u2(b) = 1/54,4271, так, что u(b) = 0,1355;

cov(a, b) = -3,1414/54,4271 = -0,0577.

Наблюдаемое значение \chi_{obs}^{2} = 2,743 с v = 4 степенями свободы. Так как \chi_{obs}^{2} не превышает \chi_{v}^{2} уровня 95%, а именно 9,488, можно считать, что модель достаточно хорошо соответствует исходным данным.

Исходные данные и линия, полученная методом взвешенной ортогональной регрессии представлены на рисунке 6. На графике также для каждого i показано положение точки (x_{i}^{*},y_{i}^{*}) на полученной прямой и точки (xi, yi). Взвешенные остатки показаны на рисунке 7.

Рисунок 6 — Данные таблицы 10 и полученная линейная калибровочная функция (см. таблицы 11 - 18)

Рисунок 7 — Взвешенные остатки, полученные в таблице 18

8. Модель, учитывающая неопределенности xi и yi и ковариации, соответствующие парам (xi, yi)

8.1. Общие положения

8.1.1 В данном разделе рассмотрена ситуация 5.3.2 c), когда имеется следующая информация для i = 1,..., m:

a) результаты измерений (xi, yi);

b) стандартная неопределенность u(xi), соответствующая xi;

c) стандартная неопределенность u(yi), соответствующая yi;

d) ковариация cov(xi, yi), соответствующая xi и yi;

В приложении D приведено руководство по определению неопределенностей и ковариаций. Все другие ковариации, соответствующие данным, считаются несущественными.

8.1.2 Ситуации 5.3.2 c) соответствует статистическая модель

x_{i} = X_{i}^{*} + d_{i} , y_{i} = Y_{i}^{*} + e_{i} , Y_{i}^{*} = A^{*} + B^{*}X_{i}^{*} , i = 1,..., m, (10)

где каждая пара (di, ei) является реализацией двумерной случайной величины с математическим ожиданием (0,0)T и ковариационной матрицей, имеющей диагональные элементы u2(xi) и u2(yi), а недиагональные элементы cov(xi, yi) = cov(yi, xi), т.е.

.

Матрица не зависит от других случайных величин.

Примечание — Предположение о том, что (di, ei) являются реализациями двумерных нормальных случайных величин, необходимо только для валидации модели (10).

8.2. Оценки параметров калибровочной функции и соответствующие стандартные неопределенности и ковариация

8.2.1 Алгоритмически данный случай является расширением (см. приложение B) обработки, приведенной в разделе 7. Вычисления в данном случае идентичны приведенным в разделе 7, кроме того, что вычисления 2) в 7.2.1 должны быть заменены на следующие:

2. t_{i} = \frac{1}{u^{2}(y_{i}) - 2\tilde{b}cov(x_{i},y_{i}) + \tilde{b}^{2}u^{2}(x_{i})} ,

${}{i}^{x^{*} =}{{(i)iiiii(i)i}^{[u^{2}y — \tilde{b}cov(x,y)]x — [cov(x,y) — \tilde{b}u^{2}x](y — \tilde{a})}}_{i}^{t}$

\zeta_{i} = y_{i} - \tilde{a} - \tilde{b}x_{i} , i = 1,..., m.

8.2.2 Все свойства, указанные в 7.2.2, применимы к данным, полученным в соответствии с моделью (10), остальную часть раздела 7 выполняют аналогично.

9. Модель, учитывающая неопределенности и ковариации, соответствующие yi

9.1. Общие положения

9.1.1 В данном разделе рассмотрена ситуация 5.3.2 d), когда имеется следующая информация для i = 1,..., m:

a) результаты измерений (xi, yi);

b) стандартная неопределенность u(yi), соответствующая yi;

c) ковариации cov(yi, yj), соответствующие парам (yi, yj), j = 1,..., m, j =/ i.

9.1.2 Квадраты стандартных неопределенностей и ковариации составляют ковариационную матрицу Uy размерности m x m для вектора результатов измерений y = (y1,..., ym)T.

.

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

9.1.3 Ситуация 5.3.2 d) соответствует статистической модели

yi = A* + B*xi + ei, i = 1,..., m, (11)

где e = (ei,..., em)T является реализацией многомерной случайной величины с вектором математического ожидания равным нулевому вектору размерности m x 1 и ковариационной матрицей Uy размерности m x m (см. [21]).

9.1.4 Оценки a и b минимизируют обобщенную сумму квадратов относительно A и B (см. [8])

, (12)

где e = y — A1 — Bx. Задача определения a и b в этом случае называется регрессионной задачей Гаусса-Маркова (см. [2]).

Примечание — В случае, когда Uy является диагональной матрицей, обобщенная сумма квадратов (12) упрощается до выражения (4) в 6.1.3 и задача сводится к задаче взвешенных наименьших квадратов.

9.2. Оценки параметров калибровочной функции, соответствующих стандартной неопределенности и ковариации

9.2.1 Если Uy положительно определенная матрица, такая, что нижняя треугольная матрица (фактор Холецкого) Ly размерности m x m существует и Uy = LyLyT (см. [10], также см. A.4), оценки a и b параметров A и B могут быть вычислены непосредственно с использованием общей схемы, приведенной в 6.2.1, после некоторых предварительных вычислений с применением матрично-векторных операций. В противном случае необходимо применение более сложных вычислений. Эти операции преобразовывают обобщенную сумму квадратов (12) в обычную сумму квадратов (2) (см. в 5.8.1), т.е. задача сводится к задаче метода невзвешенных наименьших квадратов без ковариации.

9.2.2 Оценки параметров a и b определяют в соответствии с вычислениями ₁ - ₇, а стандартные неопределенности u(a) и u(b) и ковариацию cov(a, b) определяют в соответствии с вычислением 8:

1. вычисление фактора Холецкого Ly размерности m x m, для которого Uy = LyLyT (см. A.4.1);

2. решение трех нижних треугольных систем уравнений Lyf = 1, Lyg = x и Lyh = y, где f = (f1,...fm)T и т.д. для f, g и h (см. A.4.3). 1 — вектор размерности m x 1.

3. вычисление F^{2} = \sum_{i = 1}^{m} f_{i}^{2} ;

4. вычисление g_{0} = \frac{1}{F^{2}} \sum_{i = 1}^{m} f_{i}g_{i} , h_{0} = \frac{1}{F^{2}} \sum_{i = 1}^{m} f_{i}h_{i} ;

5. вычисление \tilde{g}_{i} = g_{i} - g_{0}f_{i} , h_{i} = h_{i} - h_{0}f_{i} ; i = 1,..., m;

6. вычисление G_{2} = \sum_{i = 1}^{m} \tilde{g}_{i}^{2} ;

7. вычисление b = \frac{1}{G_{2}} \sum_{i = 1}^{m} \tilde{g}_{i}h_{i} и a = h0 — bg0;

8. вычисление u^{2}(a) = \frac{1}{F^{2}} + \frac{g_{0}^{2}}{G_{2}} , u^{2}(b) = \frac{1}{G_{2}} и cov(a,b) = \frac{g_{0}}{G_{2}} .

9.2.3 Оценки a и b, определенные в соответствии с 9.1.4, обладают следующими свойствами (см. [15]) для данных yi, соответствующих модели (11):

1. оценки a и b являются линейной комбинацией данных yi.

2. оценки a и b можно рассматривать как реализации случайных величин, математические ожидания которых равны A* и B* соответственно.

3. ковариационная матрица для случайных величин в перечислении 2) включает элементы u2(a), u2(b) и cov(a, b), вычисленные в соответствии с 9.2.2.

Свойство перечисления 1) означает, что a и b получены методом линейной оценки. Свойство перечисления 2) означает, что полученные оценки являются несмещенными. Свойства перечислений 2) и 3) показывают сходимость полученных оценок, т.е. при увеличении m, оценки a и b стремятся к A* и B* соответственно.

Метод оценки в соответствии с 9.1.4 обладает следующими оптимальными свойствами для данных yi, согласованными с моделью (11):

4. оценки \tilde{a} и \tilde{b} , полученные любым несмещенным методом линейной оценки, можно рассматривать как реализацию случайных величин, дисперсии которых не меньше, чем у оценок, полученных методом регрессионной оценки Гаусса-Маркова.

Свойство перечисления 4) может быть интерпретировано следующим образом. Для констант c и d, стандартная неопределенность u(c\tilde{a} + d\tilde{b}) , соответствующая линейной комбинации оценок \tilde{a} и \tilde{b} , полученных любым несмещенным методом линейной оценки, является не менее u(ca + db). Свойства перечислений 1) - 4) обосновывают использование метода наименьших квадратов для данных, согласующихся с моделью (11). Необходимо отметить, что эти утверждения при их использовании относятся только к математическим ожиданиям и дисперсиям ei. Соответствующие распределения далее не определяют. Если сделаны дополнительные предположения о том, что ei является реализацией случайных величин, характеризующихся многомерным нормальным распределением, то могут быть выделены следующие свойства, связанные с методом оценки Гаусса-Маркова:

5. случайные величины, указанные в перечислении 2), характеризуются двумерным нормальным распределением со средними A* и B*, ковариационной матрицей с элементами u2(a), u2(b) и cov(a, b).

6. оценки a и b являются оценками максимального правдоподобия, соответствующими наиболее вероятным значениям A и B, которые соответствуют наблюдаемым результатам измерений yi.

7. в соответствии с Байесовским анализом распределение, характеризующее знания об A и B с учетом наблюдаемых результатов измерений yi, является двумерным нормальным распределением со средними a и b и ковариационной матрицей с элементами u2(a), u2(b) и cov(a, b).

Примечание 1 — Приведенные выше свойства относятся также к методу оценки взвешенных наименьших квадратов 6.1.3 для данных, соответствующих модели (3).

Примечание 2 — Значения u2(a), u2(b) и cov(a, b), полученные при выполнении вычислений 8, определены на основе применения закона распространения неопределенности (см. Руководство ИСО/МЭК 98-3:2008) к a и b в соответствии с вычислениями 1 - 7.

9.3. Валидация модели

Если m > 2, соответствие модели данным может быть проверено при использовании взвешенных остатков ri (продолжение 9.2.2):

9. определение r_{i} = h_{i} - b\tilde{g}_{i} , i = 1,..., m;
10. определение наблюдаемого значения \chi_{obs}^{2} = \sum_{i = 1}^{m} r_{i}^{2} и числа степеней свободы v = m — 2;
11. сопоставление \chi_{obs}^{2} с квантилем \chi_{v}^{2} уровня 95%. Если \chi_{obs}^{2} превышает квантиль, то линейная модель не соответствует исходным данным.

Примечание — Тест \chi^{2} основан на предположении, что ei в модели (11) является реализацией случайных величин, характеризуемых многомерным нормальным распределением.

9.4. Организация вычислений

Вычисления в соответствии с 9.2.2 и 9.3 могут быть выполнены с применением таблиц 19 - 21. В таблице 20 приведены fi, gi и hi, вычисленные в соответствии с 1 и 2 в 9.2.2 с учетом разложения ковариационной матрицы Uy на множители Холецкого Ly. В таблице 21 значения fi, gi и hi использованы для вычисления оценок a и b параметров линейной калибровочной функции.

Таблица 19

Данные для линейной калибровочной функции Гаусса-Маркова

xi	yi
x1	y1
x2	y2
\vdots	\vdots
xm	ym

Таблица 20

Предварительные вычисления для применения метода Гаусса-Маркова

9.2.2 вычисления 1, 2
fi	gi	hi
f1	g1	h1
f2	g2	h2
\vdots	\vdots	\vdots
fm	gm	hm

Таблица 21

Организация вычислений для определения параметров линейной калибровочной функции методом Гаусса-Маркова

							9.2.2, вычисления 7;
			9.2.2 вычисления 4, 5				9.3, вычисления 9;	9.3, вычисления 10
			g0	h0			a
f_{1}^{2}	f1g1	f1h1	\tilde{g}_{1}	h_{1}	\tilde{g}_{1}^{2}	\tilde{g}_{1}h_{1}	r1	r_{1}^{2}
f_{2}^{2}	f2g2	f2h2	\tilde{g}_{2}	\tilde{h}_{2}	\tilde{g}_{2}^{2}	\tilde{g}_{2}\tilde{h}_{2}	r2	r_{2}^{2}
\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots	\vdots
f_{m}^{2}	fmgm	fmhm	\tilde{g}_{m}	\tilde{h}_{m}	\tilde{g}_{m}^{2}	\tilde{g}_{m}\tilde{h}_{m}	rm	r_{m}^{2}
F^{2} = \sum f_{i}^{2}	\sum f_{i}g_{i}	\sum f_{i}h_{i}			\tilde{G}^{2} = \sum \tilde{g}_{i}^{2}	\sum \tilde{g}_{i}h_{i}	b	\chi_{obs}^{2} = \sum_{i = 1}^{m} r_{i}^{2}

Пример - В таблице 22 приведено десять результатов измерений (xi, yi) и соответствующая матрица стандартных неопределенностей yi.

Данные получены с использованием модели, описанной в D.2.2 с uR = 1,0, uS,1 = 1,0 и uS,2 = 2,0.

Таблица 22

Десять результатов измерений yi и соответствующая ковариационная матрица

xi	yi
1,0	1,3
2,0	4,1
3,0	6,9
4,0	7,5
5,0	10,2
6,0	12,0
7,0	14,5
8,0	17,1
9,0	19,5
10,0	21,0

Ковариационная матрица Uy размерности 10 x 10, соответствующая yi

Фактор Холецкого Ly размерности 10 x 10 для Uy = LyLyT, вычисленный с использованием любого из алгоритмов, описанных в A.4.1

Векторы f, g и h в таблице 23 получены в соответствии с вычислениями 2 в 9.2.2.

Таблица 23

Таблица предварительных вычислений, соответствующих данным таблицы 22

fi	gi	hi
0,7071	0,7071	0,9192
0,4082	1,2247	2,8169
0,2887	1,7321	4,4167
0,2236	2,2361	3,9578
0,1826	2,7386	5,6963
0,4472	2,6833	5,3666
0,1491	1,6398	3,6522
0,0925	1,8490	4,4284
0,0673	2,2198	5,3208
0,0529	2,6463	5,5360

Параметры лучшей прямой, приведенные в таблице 24, вычислены в соответствии с таблицей 21. В соответствии с таблицей 24 g0 = 4,4048/1,0714 = 4,1111, h0 = 9,0048/1,0714 = 8,4044, b = 54,2185/24,6296 = 2,2014 и a = 8,4044 — (2,2014)(4,1111) = -0,6456.

Таблица 24

Таблица вычислений для данных таблицы 22

f_{i}^{2}	figi	fihi	\tilde{g}_{i}	h_{i}	\tilde{g}_{i}^{2}	\tilde{g}_{i}h_{i}	ri	r_{i}^{2}
			4,1111	8,4044			a = -0,6456
0,5000	0,5000	0,6500	-2,1999	-5,0236	4,8395	11,0514	-0,1809	0,0327
0,1667	0,5000	1,1500	-0,4536	-0,6142	0,2058	0,2786	0,3844	0,1477
0,0833	0,5000	1,2750	0,5453	1,9906	0,2973	1,0854	0,7902	0,6245
0,0500	0,5000	0,8850	1,3168	2,0785	1,7340	2,7370	-0,8202	0,6727
0,0333	0,5000	1,0400	1,9880	4,1619	3,9523	8,2739	-0,2145	0,0460
0,2000	1,2000	2,4000	0,8447	1,6080	0,7136	1,3583	-0,2516	0,0633
0,0222	0,2444	0,5444	1,0269	2,3994	1,0546	2,4640	0,1387	0,0192
0,0085	0,1709	0,4094	1,4689	3,6514	2,1578	5,3636	0,4177	0,1745
0,0045	0,1493	0,3579	1,9433	4,7555	3,7763	9,2412	0,4777	0,2282
0,0028	0,1401	0,2930	2,4287	5,0912	5,8986	12,3650	-0,2552	0,0651
1,0714	4,4048	9,0048			24,6296	54,2185	b = 2,2014	2,0740

Стандартные неопределенности и ковариация, соответствующие a и b, определены по данным таблицы 24 в соответствии с вычислениями 8 в 9.2.2:

u2(a) = 1/1,0714 + (4,1111)2/24,6296, так, чтобы u(a) = 1,2726;

u2(b) = 1/24,6296, так чтобы u(b) = 0,2015;

cov(a, b) = -4,1111/24,6296 = -0,1669.

Наблюдаемое значение \chi_{obs}^{2} = 2,074 с 8 степенями свободы определено по данным таблицы 24 с использованием 9.3. Так как \chi_{obs}^{2} не превышает квантиль \chi_{v}^{2} уровня 95%, а именно 15,507, может быть принято решение о соответствии линейной модели и данных.

Данные и полученная линейная калибровочная функция приведены на рисунке 8. Взвешенные остатки приведены на рисунке 9.

Рисунок 8 — Данные таблицы 22 и полученная линейная калибровочная функция

Рисунок 9 — Взвешенные остатки ri

10. Модель, учитывающая неопределенности и ковариации, соответствующие xi и yi

10.1. Общие положения

10.1.1 В данном разделе рассмотрена ситуация 5.3.2 e), т.е. наиболее общий случай, когда все результаты измерений имеют соответствующие неопределенности и ковариации. В приложении D приведено руководство по получению неопределенностей и ковариаций.

10.1.2 Стандартные неопределенности и ковариации являются элементами ковариационной матрицы

Размерности 2m x 2m, соответствующей вектору (x1,..., xm, y1,..., ym)T результатов измерений размерности 2m x 1.

10.1.3 Ситуация 5.3.2 e) соответствует статистической модели

x_{i} = X_{i}^{*} + d_{i} , y_{i} = Y_{i}^{*} + e_{i} , Y_{i}^{*} = A^{*} + B^{*}X_{i}^{*} , i = 1,..., m, (13)

где вектор (d1,..., dm, e1,..., em)T размерности 2m x 1 является реализацией многомерной случайной величины с нулевым вектором математического ожидания размерности 2m x 1 и ковариационной матрицей U размерности 2m x 2m (см. [21]).

10.1.4 Оценки a и b минимизируют обобщенную сумму квадратов

, (14)

где d = x — X, e = y — A1 — Bx в соответствии с A, B и Xi, i = 1,..., m. Задача определения оценок a и b в данном случае является задачей обобщенной регрессии Гаусса-Маркова (см. [2]).

10.2. Оценки параметров калибровочной функции и соответствующие стандартные неопределенности и ковариации

10.2.1 Если U положительно определенная и существует нижняя треугольная матрица L (фактор Холецкого) размерности 2m x 2m, такая что U = LLT (см. [10] и A.4), оценки a и b параметров A и B могут быть вычислены по итеративной схеме, используя матричные и векторные операции. В противном случае требуется применение более сложных методов вычислений. Эти операции преобразовывают обобщенную сумму квадратов (14) в обычную сумму квадратов (2) (см. 5.8.1), т.е. задача сводится к задаче взвешенных наименьших квадратов без ковариации. В итеративной схеме использованы приближения x_{i}^{*} , соответствующие точке на линии ( x_{i}^{*} , A + Bx_{i}^{*} ), самой близкой к точке результатов измерений (xi, yi), если близость определяется в виде взвешенного расстояния с учетом неопределенности, указанной в матрице U.

10.2.2 Оценки a и b определяют в соответствии с вычислениями 1 - 10, используя итеративную схему, аналогичную приведенной в 6.2.1, стандартные неопределенности u(a) и u(b) и ковариацию cov(a, b), выполняя вычисления 11:

1. Определяют начальные приближения \tilde{t} = (\tilde{x}_{1},...,\tilde{x}_{m},\tilde{a},\tilde{b})^{T} ;

2. Вычисляют вектор размерности 2m x 1

и матрицу размерности 2m x (m + 2) (якобиан)

,

где \tilde{x} = (\tilde{x}_{1},...,\tilde{x}_{m})^{T} , а значения \tilde{a} и \tilde{b} получены по текущему значению \tilde{t} вектора параметров;

3. Вычисляют фактор Холецкого L размерности 2m x 2m, для которого U = LLT (см. A.4.1 и [10]);

4. Определяют решение нижних треугольных систем

L\tilde{f} = f и L\tilde{J} = J ,

для определения преобразованного вектора \tilde{f} размерности 2m x 1 и преобразованной матрицы \tilde{J} размерности 2m x (m + 2) (см. A.4.3);

5. Формируют вектор g = \tilde{J}^{T}\tilde{f} размерности (m + 2) x 1 и матрицу H = \tilde{J}^{T}\tilde{J} размерности (m + 2) x (m + 2);

6. Определяют фактор Холецкого M, нижнюю треугольную матрицу размерности (m + 2) x (m + 2), для которой H = MMT (см. A.4.1);

7. Определяют решение нижней треугольной системы Mq = -g для определения вектора q размерности (m + 2) x 1 (см. A.4.3);

8. Определяют решение верхней треугольной системы M^{T}\delta t = q для определения вектора поправок \delta t размерности (m + 2) x 1 (см. A.4.4);

9. Обновление текущего приближения оценок параметров: \tilde{t}: = \tilde{t} + \delta t ;

10. Повторение вычислений 2 - 9 до тех пор, пока не достигнута необходимая сходимость. Установление a = \tilde{a} и b = \tilde{b} (элементы m + 1 и m + 2 вектора \tilde{t} );

11. Представление матрицы M, полученной при выполнении вычислений в виде:

,

,

где M22 — нижняя правая треугольная матрица размерности 2 x 2 матрицы M. Тогда

u^{2}(a) = \frac{m_{22}^{2} + m_{21}^{2}}{m_{11}^{2}m_{22}^{2}} , u^{2}(b) = \frac{m_{11}^{2}}{m_{11}^{2}m_{22}^{2}} = \frac{1}{m_{22}^{2}} , cov(a,b) = \frac{m_{21}}{m_{11}m_{22}^{2}} .

Примечание 1 — При выполнении вычислений 1 начальные приближения соответствуют вектору \tilde{t} = (x_{1},...,x_{m},a_{0},b_{0})^{T} , где a0 и b0 — значения параметров прямой, определенные с помощью метода взвешенных наименьших квадратов (см. 6.2.1).

Примечание 2 — При выполнении вычислений 8 вектор поправок \delta t уменьшается по величине при каждой итерации приблизительно в одно и то же число раз. Коэффициент уменьшения зависит от неопределенности данных: чем меньше неопределенность, тем больше сокращение. Итерации могут быть закончены, когда величина поправки станет несущественной.

Примечание 3 — При выполнении вычислений 8 поправки \delta t определяют, решая методом наименьших квадратов матричное уравнение

\tilde{J}dt = \tilde{f} .

Решение этого матричного уравнения находят из нормальных уравнений

H = \tilde{J}^{T}\tilde{J}dt = - J\tilde{f} = - g .

Примечание 4 — При выполнении вычислений 5 - 8 находят решения нормальных уравнений, используя факторизацию Холецкого. В цифровой форме более устойчивым подходом является использование факторизации QR для \tilde{J} (см. A.5.1 и [10]). В схеме, описанной в C.2, использована QR-факторизация без вычисления обратной матрицы L, как в вычислениях 4.

Примечание 5 — В матричной форме ковариационная матрица, соответствующая оценкам a и b, имеет вид:

U_{a} = M_{22}^{- T}M_{22}^{- 1} .

Примечание 6 — Более общий и численно более устойчивый подход к решению обобщенной регрессионной задачи Гаусса-Маркова в общих чертах описан в C.2. Этот подход предполагает, что матрица U положительно определенная и не включает высоких корреляций.

Примечание 7 — Значения u2(a), u2(b) и cov(a, b) получены в вычислении 11 с применением закона распространения неопределенности в соответствии с Руководством ИСО/МЭК 98-3:2008 к оценкам a и b, полученным в соответствии с вычислениями 1 - 10.

10.2.3 Поскольку оценки a и b определены в результате минимизации суммы квадратов (14), они нелинейно зависят от данных xi и yi, а следовательно, свойства для оценок обобщенного метода Гаусса-Маркова не могут быть установлены непосредственно. Оценки a и b, определенные в соответствии с 10.1.4, обладают следующими свойствами для данных xi и yi, согласованных с моделью (13):

1. Оценки a и b являются нелинейными функциями xi и yi.

2. Оценки a и b можно рассматривать как реализации случайных величин, математические ожидания которых приближенно равны A* и B*, соответственно.

3. Элементы ковариационной матрицы для случайных величин в перечислении 2) приближенно равны u2(a), u2(b) и cov(a, b), вычисленным в соответствии с 10.2.2.

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

4. Для данных, согласованных с моделью (13), при увеличении m оценки a и b сходятся к A* и B*, соответственно (см. [16]).

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

5. Случайные величины в перечислении 2) подчиняются приближенно двумерному нормальному распределению со средними A* и B* и ковариационной матрицей с элементами u2(a), u2(b) и cov(a, b).

6. Оценки a и b являются оценками максимального правдоподобия, соответствующими наибольшей вероятности значений A и B, которые соответствуют наблюдаемым результатам измерений xi и yi.

7. С позиции Байесовского анализа распределение, характеризующее знания об A и B с учетом наблюдаемых результатов измерений xi и yi, приближенно является двумерным нормальным распределением со средними a и b и ковариационной матрицей с элементами u2(a), u2(b) и cov(a, b).

10.3. Валидация модели

Если m > 2, соответствие модели исходным данным может быть частично проверено с использованием взвешенных остатков f_{i} (продолжение 10.2.2): 12) определение наблюдаемого значения \chi_{obs}^{2} = \sum_{i = 1}^{2m} f_{i}^{2} с числом степеней свободы v = m — 2; 13) сопоставление \chi_{obs}^{2} с квантилем \chi_{v}^{2} уровня 95%. Если \chi_{obs}^{2} превышает этот квантиль, то предположение о линейной модели отклоняют.

Примечание — Критерий основан на предположении, что di и ei в модели (13) в первом приближении являются реализациями случайных величин, характеризующихся многомерным нормальным распределением.

Пример — В таблице 25 приведено семь результатов измерений (xi, yi), полученных с использованием моделей измерений, описанных в D.2 и D.4.

Таблица 25

Данные семи результатов измерений xi и yi с соответствующей ковариационной матрицей

xi	yi
50,4	52,3
99,0	97,8
149,9	149,7
200,4	200,1
248,5	250,4
299,7	300,9
349,1	349,2

Ковариационная матрица, соответствующая yi, получена с использованием модели измерений (D.1) для uS = 2,0 и uR = 1,0.

Данные xi и соответствующая ковариационная матрица получены с использованием модели измерений (D.2) для z1 = 50, z2 = 100, z3 = 200, u(z1) = 0,5, u(z2) = u(z3) = 1,0, и uD,i = 0,5.

Ковариационная матрица Ux размерности 7 x 7, соответствующая xi

Фактор Холецкого Lx размерности 7 x 7 Ux = LXLXT, вычисленный с использованием любого алгоритма, описанного в A.4.1.

Ковариационная матрица Uy размерности 7 x 7, соответствующая yi

Фактор Холецкого Ly размерности 7 x 7 Uy = LyLyT, вычисленный с использованием любого алгоритма, описанного в A.4.1, имеет вид

Ковариационная матрица U размерности 14 x 14

.

Примечание - В данном примере существует корреляция, соответствующая каждой паре xi и xj и каждой паре yi и yj, отсутствует корреляция, соответствующая парам xi и yj, т.е. cov(xi, yj) = 0 для всех i и j.

Фактор Холецкого L размерности 14 x 14 для U = LLT:

.

Применение метода взвешенных наименьших квадратов к данным (6.2.1 вычисления 1 - 5) дает \tilde{a} = 0,2707 и \tilde{b} = 1,0011 . Итеративную схему начинают с \tilde{t} = (x_{1},...,x_{7},\tilde{a},\tilde{b})^{T} .

В таблице 26 приведены начальный вектор t_{0} , поправки \delta t_{k} для k-го повторения, k = 1..., 4 и заключительную оценку \tilde{t} = t_{4} .

Таблица 26

Изменения вектора \tilde{t}

t_{0}	\delta t_{1} \times 10^{- 2}	\delta t_{2} \times 10^{- 4}	\delta t_{3} \times 10^{- 6}	\delta t_{4} \times 10^{- 8}	t_{4}
50,4000	17,2531	1,2580	3,0782	0,2904	50,5727
99,0000	-43,1501	-3,2145	-6,3201	-0,7101	98,5682
149,9000	-29,1641	-3,9604	-3,8889	-0,7564	149,6080
200,4000	2,9677	-10,7629	-0,6024	-1,7165	200,4286
248,5000	24,0394	-11,4064	3,2378	-1,7064	248,7393
299,7000	-22,2510	-15,7767	-3,3581	-2,6110	299,4759
349,1000	-20,6192	-16,6217	-3,3805	-2,7429	348,8921
0,2707	7,5040	-33,3957	0,1006	-5,3019	0,3424
1,0011	0,0110	0,2113	0,0076	0,0337	1,0012

Наилучшими оценками A и B являются a = 0,3424 и b = 1,0012.

В заключительной итерации матрица M размерности 9 x 9 имеет вид:

где M22 (10.2.2, вычисление 11) имеет вид:

Стандартные неопределенности и ковариация, соответствующие a и b (см. вычисления 11 в 10.2.2):

u^{2}(a) = \frac{(110,9677)^{2} + (107,2706)^{2}}{(0,6762)^{2}(110,9677)^{2}} , u(a) = ₂,0569;

u^{2}(b) = \frac{1}{(110,9677)^{2}} , u(b) = 0,0090;

cov(a,b) = \frac{(107,2706)^{2}}{(0,6762)^{2}(110,9677)^{2}} = - 0,0129 .

Наблюдаемое значение \chi_{obs}^{2} = 1,772 с v = 5 степенями свободы получено в соответствии с вычислениями 12 в 10.3. Так как \chi_{obs}^{2} не превышает квантиль \chi_{v}^{2} уровня 95%, а именно, 11,070, это не противоречит предположению о линейности модели.

11. Использование калибровочной функции

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

11.1. Прогноз

11.1.1 Предположим, что в соответствии с применением одного из разделов 6 - 10 установлено следующее:

a) оценки параметров прямой a и b и соответствующие им стандартные неопределенности u(a) и u(b) и ковариация cov(a, b);

b) результат измерения y величины Y и соответствующая стандартная неопределенность u(y).

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

11.1.2 Оценка x величины X, соответствующая y, имеет вид:

x = \frac{y - a}{b} . (15)

11.1.3 Стандартную неопределенность u(x), соответствующую x, определяют следующим образом:

c(a) = - \frac{1}{b} , c(b) = - \frac{y - a}{b^{2}} , c(y) = \frac{1}{b} ,

u2(x) = c2(a)u2(a) + c2(b)u2(b) +

• 2c(a)c(b)cov(a, b) + c2(y)u2(y).

Примечание 1 — Формула для u2(x) получена на основе закона распространения неопределенности в соответствии с Руководством ИСО/МЭК 98-3:2008. Это приближенная формула, основанная на линеаризации формулы (15), где c(a), c(b) и c(y) — коэффициенты чувствительности.

Примечание 2 — Для вычислительных целей может быть удобно матричное представление:

, .

Примечание 3 — В случае b = 0, когда наилучшей прямой является y = a (недопустимая калибровочная функция), прогноз невозможен.

Примечание 4 — Валидация стандартной неопределенности u(x) зависит от выполнения соответствующего критерия \chi^{2} , приведенного в разделах 6 - 10.

Пример 1 — В соответствии с примером метода взвешенных наименьших квадратов с известными равными весовыми коэффициентами, приведенными в разделе 6, параметры наилучшей прямой и их стандартные неопределенности и ковариация имеют вид:

a = 1,867, b = 1,757, u(a) = 0,465,

u(b) = 0,120, cov(a, b) = -0,050.

Пусть y = 10,5 — результат дополнительного измерения Y, а u(y) = 0,5 — соответствующая стандартная неопределенность.

В соответствии с 11.1.2 оценка x величины X, соответствующей y;

x = (10,5 — 1,867)/₁,757 = 4,913.

В соответствии с 11.1.3 вычисления для определения стандартной неопределенности u(x) дают

c(a) = -1/1,867 = -0,569,

c(b) = -(10,5 — 1,867)/(1,757)2 = -2,796,

c(y) = 1/1,757 = 0,569,

u2(x) = (-0,569)2(0,217) + (-2,796)2(0,014) +

• (2)(-0,569)(-2,796)(-0,050) + (-0,569)2(0,5)2 = 0,104.

Таким образом, u(x) = 0,322.

Пример 2 — В соответствии с примером метода взвешенных наименьших квадратов с известными неравными весовыми коэффициентами, описанными в разделе 6, параметры наилучшей прямой и их стандартные неопределенности и ковариация имеют вид:

a = 0,885, b = 2,057, u(a) = 0,530,

u(b) = 0,178, cov(a, b) = -0,082.

Пусть y = 10,5 — результат дополнительного измерения Y, а u(y) = 1,0 — его стандартная неопределенность.

Из 11.1.2 оценка значений x величины X, соответствующей y, имеет вид

x = (10,5 — 0,885)/2,057 = 4,674.

В соответствии с 11.1.3 вычисления для определения стандартной неопределенности u(x) дают

c(a) = -1/0,885 = -0,486,

c(b) = -(10,5 — 0,885)/(2,057)2 = -2,272,

c(y) = 1/2,057 = 0,486,

u2(x) = (-0,486)2(0,281) + (-2,272)2(0,032) + (2)(-0,486)(-2,272)(-0,082) + (-0,486)2(1,0)2 = 0,284.

Таким образом, u(x) = 0,533.

В этом примере и примере 1 в 11.1 заметно влияние различных неопределенностей y на неопределенность x.

11.2. Предварительная оценка

Предположим, что в соответствии с применением одного из разделов 6 - 10 установлено следующее:

a) оценки параметров прямой a и b, их стандартные неопределенности u(a) и u(b) и соответствующая им ковариация cov(a, b);

b) результат измерений x величины X и его стандартная неопределенность u(x).

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

11.2.1 Оценка y величины Y, соответствующая значению x,

y = a + bx. (16)

11.2.2 Стандартную неопределенность u(y), соответствующую y, определяют при выполнении следующих вычислений:

c(a) = 1, c(b) = x, c(x) = b,

u2(y) = c2(a)u2(a) + c2(b)u2(b) +

• 2c(a)c(b)cov(a, b) + c2(x)u2(x).

Примечание 1 — Формула для u2(y) установлена с использованием закона распространения неопределенности в соответствии с Руководством ИСО/МЭК 98-3:2008. Эта аппроксимация основана на линеаризации формулы (16). Величины c(a), c(b) и c(y) представляют собой коэффициенты чувствительности.

Примечание 2 — В вычислительных целях может быть полезна матричная форма:

, .

Примечание 3 — Валидация стандартной неопределенности u(y) зависит от выполнения критерия \chi^{2} в разделах 6 - 10.

Пример — В соответствии с примером метода взвешенных наименьших квадратов с известными равными весовыми коэффициентами, приведенными в разделе 6, параметры наилучшей прямой, их стандартная неопределенность и ковариация имеют вид:

a = 1,867, b = 1,757, u(a) = 0,465,

u(b) = 0,120, cov(a, b) = -0,050.

Пусть x = 3,5 — результат дополнительного измерения X, а u(x) = 0,2 — его стандартная неопределенность и пусть cov(x, a) = cov(x, b) = 0, т.е. отсутствует корреляция x c a и x c b.

В соответствии с 11.2.1 оценка y величины Y, соответствующей x, имеет вид:

y = 1,867 + (1,757)(3,5) = 8,017.

В соответствии с 11.2.2 стандартная неопределенность u(y) имеет вид:

u2(y) = 0,217 + (3,5)2(0,014) + (2)(3,5)(-0,050) +

• (1,757)2(0,2)2 = 0,165.

Таким образом, u(y) = 0,406.

Приложение A (справочное)

ОПЕРАЦИИ С МАТРИЦАМИ

A.1 Общие положения

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

A.2 Элементарные операции

Далее используются следующие обозначения:

A — матрица результатов измерений размерности m x n с элементом A(i, j) = aij в i-й строке и j-м столбце. B — матрица размерности n x k, C — (квадратная) матрица размерности m x m, d — вектор результатов измерений размерности n x 1 с j-м элементом dj.

A.2.1 Умножение матрицы на вектор

Произведение матрицы на вектор Ad представляет собой вектор e размерности m x 1 с i-м элементом ei

e_{i} = \sum_{j = 1}^{n} a_{ij}d_{j} = a_{i1}d_{1} + a_{i2}d_{2} + ... + a_{in}d_{n} .

A.2.2 Операция умножения матрицы на матрицу

Произведение двух матриц AB представляет собой матрицу размерности m x k, j-й столбец которой является произведением матрицы A на j-й столбец B.

A.2.3 Транспонирование матрицы

Результатом AT транспонирования матрицы A является матрица размерности n x m с элементом A(j, i) = aji в j-й строке и i-м столбце.

A.2.4 Единичная матрица

Единичной матрицей порядка m является матрица l размерности m x m, у которой l(j, j) = 1, для j = 1,..., m, а все другие элементы равны нулю.

A.2.5 Инверсия квадратной матрицы

Инверсией матрицы C, если она существует, является такая матрица C-1 <*> размерности m x m, что

CC-1 = C-1C = l.

---

<*> Матрицу C-1 также называют обратной матрицей по отношению к матрице C.

Транспонирование C-1 эквивалентно инверсии CT и дает C-T.

A.3 Элементарные определения

Далее использованы следующие определения: C — (квадратная) матрица размерности m x m с элементом C(i, j) = cij в i-ой строке и j-ом столбце.

A.3.1 Симметричная матрица

Матрица C является симметричной, если cij = cji, i = 1,..., m, j = 1,..., m, т.е., C = CT.

A.3.2 Обратимая матрица

Матрица C является обратимой, если ее обратная матрица C-1 (см. 2.5), существует.

A.3.3 Нижняя треугольная и верхняя треугольная матрица

Матрица C является нижней треугольной матрицей, если cij = 0, i < j и верхней треугольной матрицей, если cij = 0, i > j.

A.3.4 Ортогональная матрица

Матрица C является ортогональной, если CTC = l.

A.4 Факторизация (разложение на множители) Холецкого

Факторизация Холецкого симметричной положительно определенной матрицы U размерности m x m — это представление матрицы в виде U = LLT (см. [10]), где L — нижняя треугольная матрица размерности m x m.

A.4.1 Алгоритмы факторизации Холецкого

A.4.1.1 Следующий алгоритм позволяет вычислить нижнюю треугольную матрицу L, такую, что U = LLT.

Инициализация

For k = 1: m

For j = k: m

L(j, k):= U(j, k)

end

end

for k = 2: m

for j = 1: k — 1

L(j, k): = 0

end

end

Факторизация

For k = 1: m

L(k,k): = L(k,k): = \sqrt{L(k,k)}

For j = k + 1: m

L(j, k):= L(j, k)/L(k, k)

end

for j = k + 1: m

for l = j: m

L(l, _j):= L(l, j) — L(l, k) L(j, k)

end

end

end

Примечание — Чтобы переписать элементы U(i, j), i >= j нижней треугольной матрицы U с разложением Холецкого, выполняют только действия стадии Разложение на множители алгоритма, приведенного в A.4.1.1, используя U вместо L.

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

Инициализация

For j = 1: m

L(j, 1: j): = U(j, 1: j)

end

for j = 1: m — 1

L(j, j + 1: m): = 0

end

Разложение на множители

для j = 1: m,

if j > 1

L(j: mj): = L(j: mj) — L(j: m, 1: j — 1) L(j, 1: j — 1))T

end

L(j:m,j) = L(j:m,j)/ \sqrt{L(j,j)} .

end

Примечание — Для получения элементов U(i, j), i >= j нижней треугольной матрицы U, имеющей разложение Холецкого, выполняют только действия этапа «Факторизация», приведенного алгоритма в A.4.1.2, используя U вместо L.

A.4.2 Интерпретация разложения Холецкого ковариационной матрицы

A.4.2.1 Пусть Ei, i = 1,..., m, — m независимых случайных величин с нулевым математическим ожиданием и единичной дисперсией каждая, а ei — реализация Ei.

y1 = l11e1,

y2 = l21e₁ + l22e₂.

Тогда u^{2}(y_{1}) = l_{11}^{2} и u^{2}(y_{2}) = l_{21}^{2} + l_{22}^{2} . Зависимость y1 и y2 от e1 означает, что у y1 и y2 существует корреляция с ковариацией cov(y1, y2) = l11l21. Далее продолжим, что

y3 = l31e1 + l32e2 + l33e3,

\vdots

ym = lm1e1 + lm2e2 + ... lmmem.

A.4.2.2 В матричной форме y = Le с нижней треугольной матрицей L. Общая зависимость y1 и y3 от e1 означает, что существует корреляция между y1 и y3. Аналогично общая зависимость y2 и y3 от e1 и e2 означает, что существует корреляция между y2 и y3 и так далее.

A.4.2.3 Для данной ковариационной матрицы U, соответствующей данным yi, разложение Холецкого U = LLT позволяет вычислить элементы матрицы таким образом, что ковариационную матрицу можно рассматривать, предполагая, что yi определены в соответствии с A.4.2.1 как реализации линейных комбинаций значений lij независимых случайных переменных Ei. На практике ковариационные матрицы часто определяют с помощью разложения на множители U = BBT, что дает U бесконечно много вариантов матриц B, которые можно использовать для построения U. Разложение на факторы Холецкого, в котором линейные комбинации представлены нижней треугольной матрицей, является однозначным с точностью до знака колонок L.

A.4.3 Решение нижней треугольной системы

A.4.3.1 Если L нижняя треугольная матрица размерности m x m такая, что L(j, j) =/ 0, j = 1,..., m, и x — вектор размерности m x 1, следующий алгоритм позволяет вычислить вектор y, где Ly = x, т.е., y = L-1x.

Инициализация

For j = 1: m

y(j): = x(j)

end

Решение

y(1): = y(1)/L(1, 1)

for j = 2: m

for k = 1: j — 1

y(j) = y(j) — L(j, k)y(k)

end

y(j): = y(j) — L(j, j)

end

Примечание — Для определения вектора x, соответствующего решению y, выполнять только действия этапа «Решение» алгоритма, приведенного в A.4.3.1, используя x вместо y.

A.4.3.2 Алгоритм, приведенный в A.4.3.1, может быть применен к решению матричного уравнения LY = X, последовательно применяя его к каждому столбцу X. Решение имеет вид Y = L-1X.

A.4.4 Решение верхней треугольной системы

A.4.4.1 Решение верхней треугольной системы может быть определено с помощью транспонирования нижней треугольной матрицы. Если L — нижняя треугольная матрица размерности m x m такая, что L(j, j) =/ 0, j = 1,..., m и x — вектор размерности m x 1, следующий алгоритм позволяет определить элемент вектора y, где y удовлетворяет уравнению LTy = x, т.е. y = L-Tx.

Инициализация

For j = 1: m

y(j): = x(j)

end

Решение

y(m): = y(m)/L(m, m)

for j = j = m — 1: — 1: 1

for k = j + 1: m

y(j): = y(j) — L(k, j) y(k)

end

y(j): = y(j)/L(j, j)

end

Примечание — Для определения вектора x, соответствующего решению y, выполняют только действия этапа «Решение» алгоритма, приведенного в A.4.4.1, используя x вместо y.

A.4.4.2 Алгоритм, приведенный в A.4.4.1, может быть применен для решения матричного уравнения LTY = X, последовательно применяя его к каждой колонке X. Решение имеет вид Y = L-TX.

A.5 Ортогональная факторизация

Ортогональные матрицы являются комбинациями вращений и отображений и имеют свойство, состоящее в том, что умножение вектора на ортогональную матрицу не изменяет длины вектора (квадратный корень из суммы квадратов элементов вектора). Столбцы ортогональной матрицы можно рассматривать как системы ортогональных осей. Важность методов ортогонального разложения состоит в том, что они позволяют решать матричные уравнения в цифровой форме устойчивым методом. Алгоритмы вычисления ортогонального разложения матрицы описаны в [1, 10, 20].

A.5.1 QR-факторизация

QR-факторизация матрицы A размерности m x n, с m >= _n, имеет вид:

,

где Q = [Q1, Q2] является ортогональной матрицей размерности m x m, Q1 — матрица, состоящая из первых n столбцов матрицы Q, Q_{1}^{T}Q_{1} = l , R1 — верхняя треугольная матрица размерности n x n.

Примечание - Также может быть получена QR-факторизация матрицы A размерности m x n, с m < n. Так как для матрицы в настоящих рекомендациях, для которых требуется QR-факторизация, необходимо выполнение неравенства m >= n, разложение не существует.

A.5.2 RQ-факторизация

A.5.2.1 RQ-факторизация матрицы B размерности m x n, с m >= n, имеет вид:

,

где Z — ортогональная матрица; T2 — верхняя треугольная матрица.

A.5.2.2 RQ-факторизация матрицы B размерности m x n, с m < n, имеет вид:

,

где Z — ортогональная матрица; T2 — верхняя треугольная матрица.

Приложение B (справочное)

ПРИМЕНЕНИЕ АЛГОРИТМА ГАУССА-НЬЮТОНА К ОБОБЩЕННОЙ РЕГРЕССИИ

B.1 В данном приложении приведены алгоритмы в соответствии с 7.2.1 и 8.2.1 с использованием алгоритма Гаусса-Ньютона.

B.2 Алгоритмы, приведенные в 7.2.1 и 8.2.1, являются частным случаем итеративного алгоритма Гаусса-Ньютона (см. [10]) для минимизации суммы квадратов нелинейных функций:

F(A) = \sum_{i = 1}^{m} f_{i}^{2}(A) , A = (A1,..., An)T, m >= n.

B.3 Пусть \tilde{a} — приближение искомого параметра a и

, ,

где f и j — соответственно вектор размерности m x 1 значений функции и якобиан размерности m x n частных производных первого порядка по параметрам, оцениваемым по приближению \tilde{a} к параметрам.

B.4 Пусть p — решение уравнения

JTJp = -JTf. (B.1)

Тогда обновленная оценка искомых параметров имеет вид \tilde{a}: = \tilde{a} + p .

B.5 Для алгоритмов в 7.2.1 и 8.2.1 A = (A, B)T и функция fi(A) является мерой обобщенного расстояния от i-й точки (xi, yi) до линии y = A + Bx.

B.6 Пусть Ui — ковариационная матрица i-й точки,

,

x_{i}^{*} \equiv x_{i}^{*}(A,B) — решение уравнения

, (B.2)

как функция A и B.

B.7 Если f_{i}^{2}(A,B) определяется равенством

f_{i}^{2}(A,B) = d_{i}^{2}(x_{i}^{*}(A,B),A,B) ,

т.е. d_{i}^{2}(x,A,B) оценивают в точке x_{i}^{*} , то значения A и B минимизируют

F(A,B) = \sum_{i = 1}^{m} f_{i}^{2}(A,B) .

Определяют лучшую линию обобщенной регрессии. Выполнение алгоритма Гаусса-Ньютона требует определения частных производных первого порядка от fi(A, B) по A и B в форме якобиана J.

B.8 Пусть n = (-B, 1)T — вектор ортогональный к линии y = A + Bx, x_{i}^{*} — решение задачи (B.2). Если xi = (xi, yi)T, x_{i}^{*} = (x_{i}^{*},A + Bx_{i}^{*})^{T} , то

ti = nTUin, (B.3)

f_{i}(A,B) = t_{i}^{- 1/2}n^{T}(x_{i} - x_{i}^{*}) , (B.4)

, . (B.5)

B.9 Решение x_{i}^{*} задачи (B.2) имеет вид:

, x_{i}^{*} = \frac{- q_{i}x_{i} + p_{i}(y_{i} - A)}{- q_{i} + p_{i}B} . (B.6)

Примечание — В выражениях (B.3), (B.4), (B.5) и (B.6) использованы Ui, а не U_{i}^{- 1} . Не требуется существования матрицы обратной к Ui, но nTUin должно быть отличным от нуля.

B.10 Алгоритмы в 7.2.1 и 8.2.1 представляют собой алгоритм Гаусса-Ньютона. В них использованы явные выражения для fi(A, B), \partial f_{i}/\partial A и \partial f_{i}/\partial B . Решение для обновления p в выражении (B.1) сформулировано как проблема определения методом взвешенных наименьших квадратов наилучшей прямой (см. 6.2.1 вычисления 1 - 5) для преобразованных данных, полученных по результатам измерений (xi, yi), соответствующих ковариационной матрице Ui и текущим аппроксимациям A и B.

Приложение C (справочное)

ПРИМЕНЕНИЕ ОРТОГОНАЛЬНОЙ ФАКТОРИЗАЦИИ К РЕШЕНИЮ ОБОБЩЕННОЙ ЗАДАЧИ ГАУССА-МАРКОВА

C.1 Общие положения

В итеративном алгоритме, описанном в 10.2.2, использовано предположение о том, что ковариационная матрица U размерности 2m x 2m является положительно определенной и, следовательно, имеет обратную матрицу. В частности, свойство обратимости предполагает, что u(xi) > 0 и u(yi) > 0. В данном приложении описан общий алгоритм, который подходит для всех действительных (симметричных положительно полуопределенных) ковариационных матриц U. Необходимо, чтобы ковариационная матрица могла быть представлена в виде U = BBT, где B — матрица размерности 2m x p(p >= m). Часто ковариационную матрицу задают в виде такого разложения на множители. Если U является обратимой, B может быть ее фактором Холецкого. Алгоритм аналогичен, описанному в 10.2.2, и требует вычисления остатков f и якобиана J, но поправки \delta t определяют, используя две ортогональных факторизации. Математически \delta t минимизирует c^{T}c = \sum_{i = 1}^{p} c_{i}^{2} при условии f = - J\delta t + Bc .

C.2 Оценки параметров калибровочной функции, соответствующие стандартной неопределенности и ковариации

Оценки a и b вычисляют в соответствии с вычислениями 1 - 9, приведенными ниже; стандартные неопределенности u(a) и u(b) оценивают в соответствии с вычислениями 10:

1. определяют начальные приближения \tilde{t} = (\tilde{x}_{1},...,\tilde{x}_{m},\tilde{a},\tilde{b})^{T} к параметрам;

2. вычисляют вектор размерности 2m x 1,

,

и якобиан размерности 2m x (m + 2),

,

где \tilde{x} = (\tilde{x}_{1},...,\tilde{x}_{m})^{T} , \tilde{a} и \tilde{b} определяют на основе текущей оценки \tilde{t} вектора параметров;

3. определяют разложение на множители QR-факторизацию матрицы J:

,

где Q — ортогональная матрица размерности 2m x 2m, R1 — верхняя треугольная матрица размерности (m + 2) x (m + 2) (см. A.5.1);

4. формируют матричное произведение QTB и находят RQ-факторизацию

QTB = TZ,

где T — матрица размерности 2m x p и Z — ортогональная матрица размерности p x p (см. A.5.2);

5. определяют \tilde{f} = Q^{T}f , \tilde{f} и T:

, ,

где f_{1} — вектор размерности (m + 2) x 1, f_{2} — вектор размерности (m — 2) x 1, T11 — матрица размерности (m + 2) x (p — m + 2), T12 — матрица размерности (m + 2) x (m — 2), и T22 — верхняя треугольная матрица размерности (m — 2) x (m — 2);

6. решают верхнюю треугольную систему T_{22}\tilde{e}_{2} = f_{2} , для определения вектора \tilde{e}_{2} = (\tilde{e}_{2,1},...,\tilde{e}_{2,m - 2})^{T} размерности (m — 2) x 1 (см. A.4.4);

7. решают верхнюю треугольную систему R_{1}\delta t = T_{12}\tilde{e}_{1} - f_{1} для определения поправок \delta t (см. A.4.4);

8. обновляют текущие приближения параметров: \tilde{t}: = \tilde{t} + \delta t ;

9. повторяют вычисления 2 - 8 до тех пор, пока не будет достигнута необходимая сходимость. Устанавливают a = \tilde{a} и b = \tilde{b} (элементы m + 1 и m + 2 из \tilde{t} );

10. пусть Ra — нижний правый минор размерности 2 x 2 из R1, Ta — нижний правый минор размерности 2 x 2 из T11. Решают верхнюю треугольную систему

RaKa = Ta,

для верхней треугольной матрицы Ka размерности 2 x 2 (см. 4.4), устанавливают U_{a} = K_{a}K_{a}^{T} , затем вычисляют

u2(a) = Ua(1, 1), u2(b) = Ua(2, 2) и cov(a, b) = Ua(1, 2).

Примечание 1 — Подход, описанный в C.2, представляет собой общее решение при определении параметров линейной калибровочной функции на основе метода наименьших квадратов. Все другие подходы, описанные в настоящих рекомендациях, являются частными случаями этого подхода.

Примечание 2 - Вычисления 1, 2, 8 и 9 в C.2 идентичны соответственно вычислениям 1, 2, 9 и 10 в 10.2.2.

C.3 Валидация модели

Если m > 2, соответствие модели данным может быть частично проверено с использованием элементов вектора \tilde{e}_{2} (продолжение C.2): 11) определяют наблюдаемое значение \chi_{obs}^{2} = \sum_{i = 1}^{m - 2} \tilde{e}_{2,i}^{2} с числом степеней свободы v = m — 2; 12) сопоставляют \chi_{obs}^{2} с квантилем \chi_{v}^{2} уровня 95%. Если \chi_{obs}^{2} превышает этот квантиль, линейную модель отклоняют.

Примечание — Критерий \chi^{2} основан на предположении, что di и ei в модели (13) представляют реализацию случайных величин, характеризуемых в первом приближении многомерным нормальным распределением. В условиях этого предположения вектор \tilde{e}_{2} размерности (m — 2) x 1 подчиняется многомерному нормальному распределению ковариационной матрицы, равной матрице идентичности размерности (m — 2) x (m — 2) так, что \chi_{obs}^{2} соответствует распределение с m — 2 степенями свободы.

Пример 1 — Подход QR-факторизации может быть применен к числовому примеру, описанному в разделе 10.

Ковариационная матрица Ux в факторизованной форме имеет вид (см. D.4)

U_{x} = B_{x}B_{x}^{T} ,

Ковариационная матрица Uy в факторизованной форме имеет вид: U_{y} = B_{y}B_{y}^{T}

Полная ковариационная матрица U размерности 14 x 14 имеет вид U = BBT, где B — матрица размерности 14 x 18

.

Для данного примера алгоритм, приведенный в C.2, математически эквивалентен алгоритму, приведенному в 10.2.2. Оба подхода дают очень близкие числовые результаты.

Пример 2 — В таблице C.1 приведено семь результатов измерений (xi, yi), полученных с использованием моделей измерений, описанных в D.2 и D.5.

Таблица C.1

Данные семи результатов измерений xi и yi

xi	yi
50,5	47,1
99,7	98,4
150,2	153,7
199,5	194,0
249,9	251,9
299,2	297,5
349,7	349,0

Ковариационная матрица, соответствующая yi, с использованием модели (D.1) с uS =2,0 и uR = 1,0 является такой же, как в примере 1 приложения C.

Данные xi и соответствующая ковариационная матрица получены с использованием модели (D.3) с z1 = 50, z2 = 100, z3 = 200, u(zi) = 0,5 и u(z2) = u(z3) = 1,0, так, что

U_{x} = B_{x}B_{x}^{T} ,

Полная ковариационная матрица U размерности 14 x 14 может быть представлена в виде U = BBT, где B — матрица размерности 14 x 11

.

Для этого примера не может быть применен алгоритм, описанный в 10.2.2, так как U не является положительно определенной. Вместо него может быть использован алгоритм, описанный в C.2.

В таблице C.2 приведены начальный вектор t_{0} , поправки \delta t_{k} для k-ой итерации k = 1,..., 5 и заключительная оценка t_{0} = t_{5} .

Таблица C.2

Изменение параметров вектора \tilde{t}

t_{0}	\delta t_{1} \times 10^{- 2}	\delta t_{2} \times 10^{- 4}	\delta t_{3} \times 10^{- 6}	\delta t_{4} \times 10^{- 8}	\delta t_{5} \times 10^{- 10}	t5
50,5000	30,8229	3,1874	23,2957	8,1124	16,4231	50,8086
99,7000	55,8313	-13,8365	26,1136	-0,2063	15,6770	100,2570
150,2000	86,6542	-10,6491	49,4093	7,9061	32,1002	151,0655
199,5000	-59,0711	-8,5976	-49,5849	-44,7904	-43,0470	198,9044
249,9000	-28,2482	-45,4102	-26,2891	-36,6780	-26,6237	249,6130
299,2000	-3,2398	-62,4341	-23,4713	-44,9967	-27,3698	299,1613
349,7000	27,5831	-59,2467	-0,1755	-36,8843	-10,9468	349,9699
-1,8528	-50,6203	-140,0856	-63,9316	-100,9432	-68,1345	-2,3731
1,0042	0,1738	0,8571	0,3217	0,6108	0,3722	1,0060

Приложение D (справочное)

ПРЕДСТАВЛЕНИЕ НЕОПРЕДЕЛЕННОСТЕЙ И КОВАРИАЦИЙ РЕЗУЛЬТАТОВ ИЗМЕРЕНИЙ x И y

D.1 Общие положения

В данном приложении установлено, как могут быть получены неопределенности и ковариации, соответствующие результатам измерений и моделируемым значениям. Подход основан на использовании модели процесса измерений, лежащей в основе определения отклика и моделируемых данных и применения закона распространения неопределенности в соответствии с Руководством ИСО/МЭК 98-3:2008. С этой целью использованы иллюстративные примеры.

D.2 Данные отклика 1

D.2.1 Общие положения

D.2.1.1 Предполагается, что величина Y, представляющая отклик средства измерений, может быть описана моделью измерений

Y = Y0 + E, (D.1)

где Y0 — величина, реализованная указанным откликом, E — величина, представляющая систематическое воздействие. Предположим, что знания об Y0 можно описать распределением со стандартным отклонением uR. Это распределение обычно основано на анализе большого количества повторений Y. Y0 оценивают средним наблюдаемых значений, uR — стандартная неопределенность, соответствующая этой оценке. Предположим, что знания о E таковы, что E обладает нулевым математическим ожиданием (т.е. были применены все необходимые корректировки) и дисперсией и u_{S}^{2} (полученной на основе знаний о средстве измерений).

D.2.1.2 Из выражения (D.1) следует, что применяя закон распространения неопределенности в соответствии с Руководством ИСО/МЭК 98-3:2008, стандартная неопределенность u(yi) измеренного значения yi величины Y имеет вид:

u^{2}(y_{i}) = u_{S}^{2} + u_{R}^{2} .

Кроме того, ковариация измеренных значений yi и yj величины Y имеет вид:

cov(y_{i},y_{j}) = u_{S}^{2} .

D.2.1.3 Таким образом ковариационная матрица в этом случае

.

D.2.2 Модель измерений для неопределенности и ковариации, соответствующих yi.

D.2.2.1 Данные, используемые в примере раздела 9, получены для измерительной системы, на которой выполнены две группы измерений. Каждая группа измерений подвергалась различным воздействиям системы. Эти воздействия являются некоррелированными, т.е.

,

где Y0,i — величина i-го отклика, E1 и E2 — величины, характеризующие воздействия системы. Предположим, что знания об Y0,i таковы, что Y0,i имеет дисперсию u_{R}^{2} , а знания об Ek таковы, что позволяют считать, что Ek обладает нулевым математическим ожиданием и дисперсией u_{S,k}^{2} для k = 1,2.

D.2.2.2 Стандартная неопределенность u(yi) результата измерений yi величины Yi имеет вид:

Ковариация значений yi и yj имеет вид:

D.2.2.3 Ковариационная матрица в этом случае имеет вид:

.

D.3 Данные отклика 2

D.3.1 Модель измерений идентична выражению (D.1) за исключением того, что воздействия системы E являются абсолютными; D — относительное воздействие системы.

Y = Y0(1 + D).

D.3.2 Обработка аналогична, проведенной в D.2, за исключением того, что теперь, используют uD для обозначения стандартной неопределенности оценки Y0,

u^{2}(y_{i}) = y_{i}^{2}u_{D}^{2} + u_{R}^{2} ,

cov(y_{i},y_{j}) = y_{i} \cdot y_{j} \cdot u_{D}^{2} .

D.3.3 Ковариационная матрица в этом случае имеет вид:

.

D.4 Данные моделирования

D.4.1 Данные, использованные в примере раздела 10, получены на основе модели измерений, в которой в соответствии с практикой, применяемой в метрологии, используют большое количество стандартных образцов для формирования большого количества значений xi при калибровке. В соответствии с моделью xi представляют собой реализации случайных величин Xi, i = 1,..., 7, определенных с помощью случайных величин Zk, k = 1, 2, 3 и Di, i = 1,..., 7:

X1 = Z1 + D1,

X2 = Z2 + D2,

X3 = Z1 + Z2 + D3,

X4 = Z3 + D4,

X5 = Z1 + Z3 + D5,

X6 = Z2 + Z3 + D6,

X7 = Z1 + Z2 + Z3D3. (D.2)

Случайные величины Zk, k = 1, 2, 3 обладают математическим ожиданием zk и дисперсией u2(zk), случайные величины Di обладают нулевым математическим ожиданием и дисперсией u_{D,i}^{2} . (При калибровке масс значения zk являются калиброванными значениями масс, а u(zk) — неопределенностями).

D.4.2 Неопределенности u(zk) и u(di) распространяют с помощью модели измерений, приведенной в D.4.1, для определения неопределенностей оценок xi, величин Xi с использованием закона распространения неопределенности в соответствии с Руководством ИСО/МЭК 98-3:2008. Общая зависимость Xi от Zk означает, что некоторые из ковариаций являются отличными от нуля. Распространение неопределенности наиболее легко может быть описано в матричной форме. Пусть

C = [CD CZ]

представляет собой матрицу чувствительности размерности 7 x 10, где CD = l — матрица идентичности размерности 7 x 7.

.

D.4.3 Пусть SD — диагональная матрица размерности 7 x 7 с диагональными элементами SD(i, i) = uD,i, i = 1,..., 7, SZ — диагональная матрица размерности 3 x 3 с диагональными элементами SZ(k, k) = u(zk), k = 1, 2, 3. Пусть

.

D.4.4 Наилучшая оценка X имеет вид x = CZz, где x — вектор размерности 7 x 1, с соответствующей матрицей ковариации размерности 7 x 7

U_{x} = B_{x}B_{x}^{T} = S_{D}^{2} + C_{Z}S_{Z}^{2}C_{Z}^{T} .

Элемент S_{D}^{2} является вкладом в дисперсию, вызванным Di, а второй элемент — вкладом в дисперсию, вызванным Zk.

D.5 Входные данные 2

D.5.1 Входные данные, используемые в примере 2 приложения C, получены на основе следующей модели измерений, связанной с описанной в D.4. Исходные данные xi являются реализацией случайных величин Xi, i = 1,..., 7, определенных как функции случайных величин Zk, k = 1, 2, 3:

X1 = Z1,

X2 = Z2,

X3 = Z1 + Z2,

X4 = Z3,

X5 = Z1 + Z3,

X6 = Z2 + Z3,

X7 = Z1 + Z2 + Z3. (D.3)

Случайные переменные Zk обладают математическими ожиданиями zk и дисперсиями u2(zk), k = 1, 2, 3.

D.5.2 Неопределенность u(zk) на основе модели измерений, приведенной в D.5.1, распространена на оценки xi величины Xi с помощью закона распространения неопределенности (см. Руководство ИСО/МЭК 98-3:2008). В соответствии с D.4 наилучшая оценка Xi имеет вид x = CZz, где x — вектор размерности 7 x 1 с соответствующей ковариационной матрицей размерности 7 x 7.

U_{x} = B_{x}B_{x}^{T} = C_{Z}S_{Z}^{2}C_{Z}^{T} , B_{x} = C_{Z}S_{Z} .

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

D.6 Исходные и наблюдаемые данные

D.6.1 Корреляции, т.е. ковариации отличные от нуля, соответствующие результатам измерений xi и yi, возникают в результате воздействий, влияющих на величины (xi и yi).

D.6.2 Предположим, что X и Y могут быть описаны моделью измерений

X = X0 + T, Y = Y0 + T, (D.4)

где X0 и Y0 и T — независимые случайные величины, математические ожидания которых равны x0, y0 и нуль и дисперсиями u2(x0), u2(y0) и u2(t) соответственно.

D.6.3 Из выражения (D.4) следует, что применяя закон распространения неопределенности в соответствии с Руководством ИСО/МЭК 98-3:2008, стандартные неопределенность u(xi) и u(yi) результатов измерений xi и yi величин X и Y имеют вид:

u2(xi) = u2(x0) + u2(t), u2(yi) = u2(y0) + u2(t).

Кроме того, ковариация xi и yi имеет вид:

cov(xi, yi) = u2(t).

D.6.4 Если вместо X и Y может быть применена модель измерений

X = X0 + T, Y = Y0 — T,

ковариация xi и yi имеет вид

cov(xi, yi) = -u2(t).

Приложение E (справочное)

НЕОПРЕДЕЛЕННОСТЬ, ИЗВЕСТНАЯ С ТОЧНОСТЬЮ ДО ПОСТОЯННОГО МНОЖИТЕЛЯ

E.1 В данном приложении описан метод оценки неопределенности результатов измерений в случае, если неопределенность известна с точностью до постоянного коэффициента.

E.2 В настоящих рекомендациях предполагается, что величины (переменные) характеризуются согласно принципам Руководства ИСО/МЭК 98-3:2008 и Руководства ИСО/МЭК 98-3:2008/Дополнение 1:2008 (см. [13]) распределением вероятностей. Результатам измерений соответствует математическое ожидание и дисперсия соответствующей случайной величины.

Примечание - Для измеренного значения y, как реализации случайной величины, характеризуемой t-распределением с параметром масштаба s и v степенями свободы (v > 2), стандартная неопределенность имеет вид u(y) = [v/(v — 2)]1/2, где s — стандартное отклонение распределения.

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

В настоящих рекомендациях применен этот подход в ситуации, когда неопределенность известна с точностью до постоянного коэффициента. В самом общем случае (раздел 9) предполагается, что измеренные значения y имеют идентичные неопределенности, но их общая стандартная неопределенность \sigma неизвестна. (Это — более общий случай, когда ковариационная матрица имеет вид U = \sigma^{2}U_{0} , где U0 — задано, \sigma неизвестно.) Если m > 2, можно получить оценку \hat{\sigma} на основе разброса исходных данных вокруг подобранной линии калибровочной функции. Эта оценка известна как апостериорная оценка \sigma , квалификация апостериорного отношения к данным после определения наилучшей линии калибровки.

E.4 Апостериорные оценки определяют используя те же понятия, что и при валидации модели. Делая предположение о том, что исходные данные являются реализацией величины, характеризуемой многомерным нормальным распределением, апостериорную оценку \hat{\sigma} выбирают так, чтобы \chi_{obs}^{2} был равен (m — 2) — (математическому ожиданию \chi^{2} - \text{распределения} с (m — 2) степенями свободы). В этом случае валидация модели не может быть выполнена, так как апостериорная оценка выбрана так, чтобы критерий валидации был выполнен.

E.5 Эти методы должны быть использованы с большой осторожностью. Например, если график данных указывает, что они не соответствуют линейной калибровочной функции, метод не должен быть использован.

E.6 Оценки параметров a и b не зависят от коэффициента пропорциональности \sigma . Оценка \sigma необходима только для определения стандартных неопределенностей u(a) и u(b) и ковариации cov(a, b). Для случая, когда U полностью известна, u(a), u(b) и cov(a, b) могут быть оценены по данным и U. Предположения о распределениях данных не нужны. В предположении о нормальности оценки параметров могут рассматриваться как реализации случайных величин, характеризуемых двумерным распределением.

E.7 В случае, если данные могут быть рассмотрены, как реализации многомерного нормального распределения с известной ковариационной матрицей U, двумерное распределение в E.6 является нормальным с ковариационной матрицей Ua и элементами u2(a) и u2(b) и cov(a, b) как в выражении (1).

E.8 Для случая E.4, когда многомерному нормальному распределению соответствует ковариационная матрица U = \sigma^{2}U_{0} , где U0 известно, а \sigma неизвестно, U0 используют вместо U в вычислениях. Ковариационная матрица

оценок параметров линейной калибровочной функции может быть вычислена. Если m > 2, наблюдаемое значение \chi^{2} (см. 6.3) может быть использовано для определения апостериорной оценки коэффициента пропорциональности, связанного с исходной неопределенностью. Пусть \chi_{obs}^{2} определено в соответствии с вычислениями 8 в 6.3

\hat{\sigma}^{2} = \chi_{obs}^{2}/(m - 2) .

E.9 Скорректированная по масштабу ковариационная матрица

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

\widehat{U}_{a} = \hat{\sigma}^{2}U_{a,0} .

т.е. скорректированные по масштабу стандартные неопределенности \hat{u}(a) и \hat{u}(b) и ковариация \widehat{c}ov(a,b) полученных оценок параметров имеют вид:

\hat{u}^{2}(a) = \hat{\sigma}^{2}u_{0}^{2}(a) , \hat{u}^{2}(b) = \hat{\sigma}^{2}u_{0}^{2}(b) , \widehat{c}ov(a,b) = \hat{u}^{2}cov_{0}(a,b) . (E.1)

E.10 Оценки (E.1) основаны на конечном количестве m исходных данных. Для небольших m метод занижает значение дисперсии распределения для оценок параметров. Для m > 4 лучшую оценку определяют (см. [19, глава 8]) используя

\tilde{\sigma}^{2} = \frac{m - 2}{m - 4} \frac{\chi_{obs}^{2}}{m - 2} = \frac{\chi_{obs}^{2}}{m - 4} .

Примечание - В случае предположения о нормальности распределения оценки параметров рассматривают как реализацию случайной величины с двумерным t-распределением с матрицей параметров масштаба U_{a} и (m — 2) степенями свободы. Для m > 4 ковариационная матрица этого распределения имеет вид:

, (E.2)

где коэффициент (m — 2)/(m — 4) учитывает то, что \sigma является оценкой, а не известным значением.

Пример - (неизвестные весовые коэффициенты). В данном примере xi определены точно, а yi имеют равные, но неизвестные стандартные неопределенности, апостериорные оценки неопределенностей подобранных параметров оценены на основе полученных остатков. Аппроксимацию проводят, выбирая весовые коэффициенты, равные единице (это значит, что стандартная неопределенность u(yi) также равна единице). Данные приведены в таблице E.1.

Таблица E.1

Данные шести результатов измерений

xi	yi	u(yi)
1,000	3,014	1
2,000	5,225	1
3,000	7,004	1
4,000	9,061	1
5,000	11,201	1
6,000	12,762	1

В таблице E.2 приведены результаты вычислений параметров наилучшей прямой. В соответствии с этой таблицей g0 = 21,000/6,000 = 3,500, h0 = 48,267/6,000 = 8,044, b = 34,363/17,500 = 1,964, a = 8,044 — (1,964)(3,500) = 1,172.

Таблица E.2

Вычисления на основе данных таблицы E.1

wi	w_{i}^{2}	w_{i}^{2}x_{i}	w_{i}^{2}y_{i}	gi	hi	g_{i}^{2}	gihi	ri	r_{i}^{2}
				3,500	8,044			a = 1,172
1,000	1,000	1,000	3,014	-2,500	-5,031	6,250	12,576	-0,122	0,015
1,000	1,000	2,000	5,225	-1,500	-2,819	2,250	4,229	0,126	0,016
1,000	1,000	3,000	7,004	-0,500	-1,040	0,250	0,520	-0,059	0,003
1,000	1,000	4,000	9,061	0,500	1,017	0,250	0,508	0,035	0,001
1,000	1,000	5,000	11,201	1,500	3,157	2,250	4,735	0,211	0,045
1,000	1,000	6,000	12,762	2,500	4,718	6,250	11,794	-0,191	0,037
	6,000	21,000	48,267			17,500	34,363	b = 1,964	0,116

Данные и полученная линейная калибровочная функция приведены на рисунке E.1. Взвешенные остатки показаны на рисунке E.2. Поскольку всем u(yi) присвоено значение 1, в этом случае неопределенности сильно превышают остатки по величине.

Рисунок E.1 — Данные таблицы E.1 и полученная линейная калибровочная функция (см. таблицу E.2)

Рисунок E.2 — Взвешенные остатки для подобранной линейной калибровочной функции в соответствии с таблицей E.2

Если априорно известно, что u(yi) = 1, i = 1,..., m, то неопределенность полученных параметров, вычисленная по данным таблицы E.2 имеет вид:

u2(a) = 1/6,000 + (3,500)2/17,500, так, что u(a) = 0,931;

u2(b) = 1/17,500, так, что u(b) = 0,239;

cov(a, b) = -3,500/17,500 = -0,200.

Поскольку эти вычисления основаны на произвольном присвоении u(yi) = 1, необходима апостериорная оценка \hat{\sigma} по u(yi) для определения неопределенностей параметров полученной функции. В соответствии с таблицей

\hat{\sigma}^{2} = \frac{\chi_{obs}^{2}}{m - 2} = \frac{0,116}{4} = 0,029 , \hat{\sigma} = 0,171 .

Значение \hat{\sigma} представляет собой оценку стандартной неопределенности u(yi), соответствующей yi, основанную на наблюдаемом значении \chi_{obs}^{2} . Учитывая эту апостериорную оценку, вычисления могут быть повторены с u(yi) = 0,171. Оценки a и b при этом не изменились, но наблюдаемое значение \chi_{obs}^{2} и неопределенности определены следующим образом:

r_{i} = \frac{r_{i}}{\hat{\sigma}} ,

так, что \frac{\chi_{obs}^{2}}{\hat{\sigma}^{2}} = m - 2 = 4 — математическое ожидание \chi^{2} - \text{распределения} с 4 степенями свободы. Используя формулу (E.1) можно записать

\hat{u}^{2}(a) = \hat{\sigma}^{2}u_{0}^{2}(a) = 0,867\hat{\sigma}^{2} = 0,025 , так, что \hat{u}(a) = 0,931\hat{\sigma} = 0,159 ;

\hat{u}(b) = \hat{\sigma}^{2}u_{0}^{2}(b) = 0,057\hat{\sigma}^{2} = 0,002 , так, что \hat{u}(b) = 0,239\tilde{\sigma} = 0,041 ;

\widehat{c}ov(a,b) = \hat{\sigma}^{2}cov_{0}(a,b) = - 0,200\hat{\sigma}^{2} = - 0,006 .

Элементы матрицы \widehat{U}_{a} затем оценивают, если априорно известно, что u(y_{i}) = \hat{\sigma} . Однако, \hat{\sigma} является оценкой стандартной неопределенности yi. Для m > 4 коэффициент (m — 2)/(m — 4) может быть включен в ковариационную матрицу для учета дополнительной неопределенности, которая является результатом того, что оценки \hat{\sigma} получены по m значениям. Используя формулу (E.2), можно записать

,

\tilde{u}(a) = (0,050)^{1/2} = 0,225 , \tilde{u}(b) = (0,003)^{1/2} = 0,058 и \widetilde{c}ov(a,b) = - 0,012 .

Приложение F (справочное)

РАЗРАБОТКА ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ ДЛЯ ОПИСАННЫХ АЛГОРИТМОВ

F.1 Программное обеспечение, реализующее алгоритмы, описанные в настоящих рекомендациях для определения и использования линейных калибровочных функций, разработано Национальной Физической Лабораторией (NPL) Соединенного Королевства. Программное обеспечение доступно как архивированная папка ZIP с веб-сайтов NPL www.npl.co.uk/mathematics-scientific-computing/software-support-for-metrology/software-downloads-(ssfm) http://www.npl.co.uk/mathematics-scientific-computing/software-support-for-metrology/software-downloads-(ssfm) и ИСО standards.iso.org/iso/ts/28037/.

F.2 Программное обеспечение разработано на языке программирования MATLAB (см. [18]), в форме M-файлов и файлов html и использует Версии 7.10.0 MATLAB (R2010a). Для пользователей MATLAB M-файлы могут быть выполнены непосредственно и также изменены для выполнения алгоритмов обработки различных данных. Для пользователей, не имеющих доступа к MATLAB, программное обеспечение более всего подходит для использования как файлы html. Программное обеспечение может быть использовано в качестве основы для подготовки выполнения алгоритмов на других языках программирования. В пределах файлов использованы обращения ко многим функциям MATLAB, которые также включены в программное обеспечение. Например, функция algm_gdr1_steps_2_to_5 выполняет вычисления 2 - 5 процедуры вычислений ситуации 5.3.2 b) (неопределенность соответствует xi и yi, а все ковариации, соответствующие данным, являются несущественными), установленной в 7.2.1. Кроме того, некоторое использование встроенных функций MATLAB предусмотрено для разложения Холецкого. Скрипты MATLAB (имеющие расширение '.m') и файлы html ('.html') обеспечены следующим:

• TS28037_WLS1 (выполняет числовой пример метода взвешенных наименьших квадратов с известными равными весовыми коэффициентами, описанный в разделе 6 и выполняет прогноз, описанный в 11.1, пример 1 и предварительную оценку, описанную в 11.2);
• TS28037_WLS2 (выполняет числовой пример метода взвешенных наименьших квадратов с известными неравными весовыми коэффициентами, описанный в разделе 6, и выполняет прогноз, описанный в 11.1, пример 2);
• TS28037_WLS3 (выполняет числовой пример метода взвешенных наименьших квадратов с неизвестными равными весовыми коэффициентами, описанный в приложении E);
• TS28037_GDR1 (выполняет числовой пример обобщенного регрессионного анализа расстояний, описанный в разделе 7);
• TS28037_GDR2 (выполняет числовой пример, иллюстрирующий алгоритм для обобщенной регрессии расстояний, описанный в разделе 8);
• TS28037_GMR (выполняет числовой пример регрессии Гаусса-Маркова (GMR), описанный в разделе 9);
• TS28037_GGMR1 (выполняет числовой пример обобщенной регрессии Гаусса-Маркова, описанный в разделе 10);
• TS28037_GGMR2 (выполняет числовой пример обобщенной регрессии Гаусса-Маркова, описанный в разделе 10 и приложении C, Пример 1, с использованием ортогонального разложения, описанного в C.2);
• TS28037_GGMR3 (выполняет числовой пример, описанный в приложении C, Пример 2, с использованием ортогонального разложения, описанного в C.2).

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

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

F.4 Предоставлено соглашение о лицензии на программное обеспечение, имеется лицензионное соглашение (REF:MSC/L/10/001) и использование программного обеспечения должно соответствовать правовым требованиям этого соглашения. Используя MATLAB, пользователь принимает условия соглашения. Запросы на программное обеспечение следует направлять в NPL по адресу enquiries@npl.co.uk.

Приложение G (справочное)

ПЕРЕЧЕНЬ ОСНОВНЫХ УСЛОВНЫХ ОБОЗНАЧЕНИЙ

A	-	точка пересечения линейной калибровочной функции с осью абсцисс;
A*	-	неизвестное значение параметра A для конкретной измерительной системы;
a	-	оценка параметра A;
a	-	вектор (a, b)T оценок параметров калибровочной функции;
B	-	угловой коэффициент линейной калибровочной функции;
B*	-	неизвестное значение параметра B для конкретной измерительной системы;
b	-	оценка параметра B;
cov(a, b)	-	ковариация оценок a и b;
di	-	разность (x_{i} - X_{i}^{*}) , представляющая собой реализацию случайной величины с нулевым математическим ожиданием и дисперсией u2(xi);
ei	-	разность (y_{i} - Y_{i}^{*}) , представляющая собой реализацию случайной величины с нулевым математическим ожиданием и дисперсией u2(yi);
L	-	нижняя треугольная матрица;
m	-	количество точек результатов измерений;
ri	-	взвешенный остаток или взвешенное расстояние i-й точки для оценок a и b;
Ri	-	взвешенный остаток или взвешенное расстояние i-й точки для параметров A и B;
U	-	матрица ковариации размерности 2m x 2m, соответствующая результатам измерений (xi, yi), i = 1,..., m;
Ua	-	ковариационная матрица a размерности 2 x 2;
Ux	-	ковариационная матрица размерности m x m, соответствующая результатам измерений xi, i = 1,..., m;
Uy	-	ковариационная матрица размерности m x m, соответствующая результатам измерений yi, i = 1,..., m;
uR	-	стандартное отклонение случайной величины с распределением, отражающим знания о случайных воздействиях;
uS	-	стандартное отклонение случайной величины с распределением, отражающим знания о влиянии системы;
u(z)	-	стандартная неопределенность z, в качестве z могут быть использованы a, b, xi, yi и т.д.;
\upsilon_{i}	-	величина, обратная u(xi);
wi	-	величина, обратная u(yi);
X	-	независимая величина (переменная);
Xi	-	i-я независимая величина (переменная);
X_{i}^{*}	-	неизвестное значение i-й независимой величины, обеспечиваемое измерительной системой;
x	-	оценка X (в случае прогноза) или результат измерений X (предварительная оценка);
xi	-	i-й результат измерений X;
x_{i}^{*}	-	оценка i-й независимой величины переменной;
Y	-	зависимая величина (переменная);
Yi	-	i-я зависимая величина;
Y_{i}^{*}	-	неизвестное значение i-й зависимой величины, измерения которой обеспечивает измерительная система;
y	-	результат измерений величины Y (в случае прогноза) или оценка Y (предварительная оценка);
yi	-	i-й результат измерений величины Y;
y_{i}^{*}	-	оценка i-й зависимой величины;
v	-	число степеней свободы для модели \chi^{2} - \text{распределения} или t-распределения;
\sigma	-	стандартное отклонение случайной величины, характеризуемой распределением вероятностей;
\hat{\sigma}	-	апостериорная оценка \sigma ;
\chi_{obs}^{2}	-	наблюдаемое значение \chi^{2} ;
\chi_{v}^{2}	-	случайная величина, подчиняющаяся \chi^{2} - \text{распределению} с v степенями свободы.

Приложение ДА (справочное)

СВЕДЕНИЯ О СООТВЕТСТВИИ ССЫЛОЧНЫХ МЕЖДУНАРОДНЫХ СТАНДАРТОВ ССЫЛОЧНЫМ НАЦИОНАЛЬНЫМ СТАНДАРТАМ РОССИЙСКОЙ ФЕДЕРАЦИИ

Таблица ДА.1

Обозначение ссылочного международного документа	Степень соответствия	Обозначение и наименование соответствующего национального стандарта
Руководство ИСО/МЭК 99:2007	-	<*>
Руководство ИСО/МЭК 98-3:2008	IDT	ГОСТ Р 54500.3-2011/Руководство ИСО/МЭК 98-3:2008 «Неопределенность измерения. Часть 3. Руководство по выражению неопределенности измерения»
Руководство ИСО/МЭК 98-3:2008/, Дополнение 1:2008	IDT	ГОСТ Р 54500.3.1-2011/Руководство ИСО/МЭК 98-3:2008/Дополнение 1:2008 «Неопределенность измерения. Часть 3. Руководство по выражению неопределенности измерения по выражению неопределенности измерения. Дополнение 1. Трансформирование распределений с использованием метода Монте-Карло»
<*> Соответствующий национальный стандарт отсутствует. До его утверждения рекомендуется использовать перевод на русский язык данного международного стандарта. Перевод данного международного стандарта находится в Федеральном информационном фонде технических регламентов и стандартов. Примечание — В настоящей таблице использовано следующее условное обозначение степени соответствия стандартов: - IDT — идентичные стандарты.

БИБЛИОГРАФИЯ

[1]	Anderson, E., Bai, Z., Bischof, C.H., Blackford, S., Demmel, J., Dongarra, J.J., Croz, J.D.A. Greenbaum, Hammarling, S., McKenney, A., and Sorensen, D.C. LAPACK Users' Guide, 3rd ed., SIAM, Philadelphia, PA, 1999, http://www.netlib.org/lapack/lug/
[2]	Bartholomew-Biggs, M., Butler, B.P., and Forbes, A.B. Optimisation algorithms for generalised regression on metrology, In Advanced Mathematical and Computational Tools in Metrology IV (Singapore, 2000), P. Ciarlini, A.B. Forbes, F. Pavese, and D. Richter, Eds., World Scientific, pp. 21 — 31
[3]	Boggs, P.T., Byrd, R.H., and Schnabel, R.B. A stable and efficient algorithm for nonlinear orthogonal distance regression, SIAM J. Sci. Stat. Comput. 8,6 (1987), 1052 — 1078
[4]	Butler, B.P., Cox, M.G., Ellison, S.L.R., and Hardcastle, W.A., Eds., Statistics Software Qualification: Reference Data Sets, Royal Society of Chemistry, Cambridge, 1996
[5]	Carroll, R.J., Ruppert, D., and Stefanski, L. a. Measurement error in nonlinear models, Chapman&Hall/CRC, Boca Raton, 1995
[6]	Cox, M.G., Forbes, A.B., Harris, P.M., and Smith, I.M. The classification and solution of regression problems for calibration, Tech. Rep. CMSC 24/03, National Physical Laboratory, Teddington, UK, 2003
[7]	Draper, N.R., and Smith, H. Applied Regression Analysis, Wiley, New York, 1998, Third edition
[8]	Forbes, A.B., Harris, P.M., and Smith, I.M. Generalised Gauss-Markov regression, In Algorithms for Approximation IV (Huddersfield, UK, 2002), J. Levesley, I. Anderson, and J.C. Mason, Eds., University of Huddersfield, pp. 270 — 277
[9]	Fuller, W. a. Measurement Error Models, Wiley, New York, 1987
[10]	Golub, G.H., and Van Loan, C.F. Matrix Computations, North Oxford Academic, Oxford, 1983
[11]	ISO 3534-1:2006	Statistics — Vocabulary and symbols — Part 1: General statistical terms and terms used in probability
[12]	ISO 3534-2:2006	Statistics — Vocabulary and symbols — Part 2: Applied statistics
[13]	ISO/IEC Guide 98-3/Suppl. 1, Uncertainty of measurement — Part 3: Guide to the expression of uncertainty in measurement (GUM: 1995) — Supplement 1: Propagation of distributions using a Monte Carlo method
[14]	ISO 11095:1996	Linear calibration using reference materials
[15]	Kendall, M.G., and Stuart, a. The Advanced Theory of Statistics, Volume 2: Inference and Relationship. Charles Griffin, London, 1961
[16]	Kukush, A., and Van Huffel, S. Consistency of elementwise-weighted total least squares estimator in a multivariate errors-in-variables model AX-B, Metrika 59, 1 (February 2004), 75 — 97
[17]	Mardia, K.V., Kent, J.T., and Bibby, J.M. Multivariate Analysis, Academic Press, London, 1979
[18]	MATLAB http://www.mathworks.com/products/matlab/
[19]	Migon, H.S. and Gamerman, D. Statistical Inference: An Integrated Approach, Arnold, London, 1999
[20]	Paige, C.C. Fast numerically stable computations for generalized least squares problems, SIAM J. Numer. Anal. 16 (1979), 165 — 171
[21]	Strang, G., and Borre, K. Linear Algebra, Geodesy and GPS, Wiley, Wellesley-Cambridge Press, 1997

---