Вы также можете заказать у нас разработку методик оценивания неопределенности измерений
Введение
При обработке результатов измерений, получаемых при испытаниях, калибровке или поверке средств измерений, при проведении анализа образцов продукции, в качестве характеристики качества результата измерений используют неопределенность измерений.
В настоящее время выделяют три надежных подхода по количественной оценке и расчету неопределенности измерений:
– метод моделирования, с применением закона распределения неопределенности;
– метод моделирования Монте-Карло;
– эмпирические методы, основанные на использовании оценок повторяемости, воспроизводимости и правильности.
Метод Монте-Карло (ММК) используется как метод трансформирования распределений на основе моделирования случайных выборок из этих распределений. Этот метод может быть применен к любым моделям, имеющим единственную выходную величину, в которых входные величины характеризуются любыми заданными функциями распределения вероятностей.
Согласно ГОСТ 34100.3.1 для оценивания неопределенности измерений метод Монте-Карло используется, когда:
– вклад разных составляющих неопределенности существенно неодинаков;
– распределение выходной величины не является нормальным или масштабированным смещенным t - распределением;
– трудно найти частные производные от функции измерения, как того требует закон трансформирования неопределенностей;
– модель достаточно сложная;
– плотности распределения вероятностей входных величин асимметричны;
– оценка выходной величины и соответствующая стандартная неопределенность имеют приблизительно равные значения.
Применение метода Монте-Карло
Этот метод позволяет на основе аналитического описания объекта моделировать поведение выходной величины с учётом случайного изменения входных величин. Для этого с помощью генератора случайных задаются значения входных величин Xi для первого испытания, в зависимости от их функции плотности вероятности (ФПВ), и вычисляется значение выходной величины Y. Затем аналогично для второго испытания и т.д. В ходе проведения расчета неопределенности результатов измерений происходит трансформирование ФПВ входных величин сквозь разработанную модель, связывающая Y и Xi, в результате получают совокупность значений выходной величины, которая описывает ФПВ для Y. С помощью статистического анализа оценивается закон распределения и параметры выходной величины.
Поскольку ММК требует проведения большого числа испытаний, его часто называют методом статистических испытаний.
Оценка входной величины
Следует выбрать M — количество оцениваний модели, которое необходимо произвести. Число попыток, необходимых для получения этих результатов с предписанной числовой точностью, зависит от «формы» ФПВ выходной величины и требуемой вероятности охвата.
Лучше всего выбирать достаточно большие значение М (например, превышающие в 104 раз по сравнению с 1/(1-p)). Тогда можно ожидать, что функция распределения G обеспечит приемлемое дискретное представление GY(η) (функции распределения для Y) вблизи границ 100 p %-ного интервала охвата для Y.
Так как точно не известно, что это или любое заранее предписанное число будет достаточным, может быть использована процедура адаптивного выбора M, т. е. по мере выполнения попыток.
Для применения метода формируют М векторов xr, r=1, …, M в соответствии с плотностями распределения вероятностей gXi(ξi) для N входных величин Xi или, если это необходимо, из совместной плотности распределения gX(о). Модель оценивается для каждого из М извлечений из функции плотности вероятности (ФПВ) для N входных величин.
Конкретнее, необходимо обозначить M извлечений через x1, …, xM, где r-е значение xr состоит из случайных значений x1,r, …, xN,r, а xi,r – случайное значение из ФПВ для Xi. Тогда значения модели можно представить в виде:
yr = f(xr), r = 1, …, M.
Дискретное представление G функции распределения GY(η) для выходной величины Y может быть получено следующим образом:
а) необходимо рассортировать значения модели yr, r = 1, …, M, полученные по ММК, в неубывающем порядке. Обозначить рассортированные значения модели как y(r), r = 1, …, M;
б) если необходимо, создать возмущения для любых дублирующих значений модели y(r) так, чтобы конечный полный набор y(r), r = 1, …, M, формировал строго возрастающую последовательность;
в) полученная последовательность y(r), r = 1, …, M определяет G.
Если выходная величина Y будет рассматриваться как входная величина при оценивании неопределенности другого измерения, то выборку из ее распределения легко получить случайным (равновероятным) выбором значений из y(r), r = 1, …, M.
Последовательность y(r) (или yr) может быть скомпонована в гистограмму (при подходящей ширине ячеек) и представлять собой частотное распределение, которое, при условии, что оно нормировано, чтобы иметь единичную площадь, обеспечивает приближение к ФПВ gY(η) для Y. Вычисления характеристик распределения обычно проводятся не в терминах этой гистограммы, разрешение которой зависит от выбора ширины ячеек, а в терминах G. Тем не менее, гистограмма может быть полезна для понимания природы ФПВ, например степени ее асимметрии. В ряде случаев полезна аппроксимация GY(η) непрерывной функцией.
Выбор распределения вероятностей
Приписывание, в некоторых общих случаях, ФПВ входным величинам Xi на этапе постановки задачи оценивания неопределенности может быть основано на теореме Байеса или на принципе максимума энтропии.
Чаще всего выбор вида плотности распределения вероятностей осуществляет метролог, составляющий модель измерений, который в конечном счете несет ответственность за качество конечных результатов.
Плотность распределения вероятностей случайных величин можно выбрать на основе имеющейся о них информации. В таблице 1 описана информация и ФПВ, которая приписывается на основании данной информации.
|
Информация о величине |
Приписанная ФПВ |
Изображение ФПВ |
|---|---|---|
|
Нижняя и верхняя границы a, b |
Прямоугольное: R(a, b) |
 |
|
Неточные нижняя и верхняя границы a ± d, b ± d |
Криволинейная трапеция: CTrap(a, b, d) |
 |
|
Сумма двух величин с приписанными прямоугольными распределениями с нижними и верхними границами a1, b1 и a2, b2 |
Трапецеидальное: Trap(a, b, β) a = a1 + a2, b = b1 + b2, β = |(b1 – – a1) – (b2 – a2)|/(b − a) |
 |
| Сумма двух величин с приписанными прямоугольными распределениями с нижними и верхними границами a1, b1 и a2, b2 и одинаковой полушириной (b1 – a1) = (b2 – a2) |
Треугольное: T(a, b) a = a1 + + a2, b = b1 + b2 |
|
| Синусоидальное циклическое изменение между нижней и верхней границами a, b | Арксинусное (U — образное): U(a, b) |
|
| Наилучшая оценка x и связанная с ней стандартная неопределенность u(x) | Гаусса: N(x, u2 (x)) |
|
| Наилучшая оценка x векторной величины и соответствующая матрица неопределенностей Ux | Многомерное Гаусса: N(x,Ux) | - |
| Серия выбранных независимо наблюдений x1, …, xn величины, имеющей распределение Гаусса с неизвестным ожиданием и неизвестной дисперсией |
|
|
| Наилучшая оценка x, расширенная неопределенность Up, коэффициент охвата kp и эффективные степени свободы veff | Масштабированное и сдвинутое t: tVeff (x ,(Up /kp )2) |
|
| Наилучшая оценка x неотрицательной величины | Экспоненциальное: Ex(1/x) |
|
| Количество q подсчитываемых объектов | Гамма: G(q + 1,1) |
|
Далее будут разобраны наиболее часто применяемые основные законы распределения:
– равномерное (прямоугольное);
– треугольное;
– нормальное (Гаусса) и др.
Прямоугольное распределение
Если единственной доступной информацией о величине X являются нижняя граница a и верхняя граница b, a < b, тогда, согласно принципу максимума энтропии, X будет приписано прямоугольное распределение R(a,b) на интервале [a,b].
ФПВ для X следующая
(1)
X имеет ожидание и дисперсию
(2)
Треугольное распределение
Предположим, что величина X определена как сумма двух независимых величин, каждой из которых приписано прямоугольное распределение, но с равными полуширинами, т. е. b1 – a1 = b2 – a2. Распределение для X — это трапецеидальное распределение Trap(a, b, 0), которое редуцируется до треугольного распределения T(a, b) на интервале [a, b].
Функцией плотности вероятности для Х является
где x = (a + b)/2 и ω = (b – a)/2.
Формула (3) может быть выражена как
(4)
*для обработки на компьютере.
Величина X имеет ожидание и дисперсию
(5)
Распределение Гаусса
Если наилучшая оценка x и соответствующая стандартная неопределенность u(x) являются единственной доступной информацией, касающейся величины X, тогда, согласно принципу максимума энтропии, величине X должно быть приписано Гауссовское распределение вероятностей N(x, u2 (x)).
Для Х функция плотности вероятности имеет вид
(6)
Величина X имеет следующее ожидание и дисперсию
(7)
Представление результатов
Обычно после использования трансформирования распределений представляют следующее:
– оценку y выходной величины Y;
– стандартную неопределенность u(y), связанную с y;
– оговоренную вероятность охвата 100p % (например, 95 %);
– границы выбранного 100p % интервала охвата (например, 95 % интервала охвата) для Y;
– любую другую значимую информацию, такую как: интервал охвата — это вероятностно симметричный интервал или наименьший интервал охвата.
Оценка выходной величины
Оценкой выходной величины Y является выборочное среднее
(8)
В качестве оценки стандартной неопределенности Y используют выборочное стандартное отклонение
(9)
При некоторых особых обстоятельствах, когда одной из входных величин приписано t-распределение с числом степеней свободы менее трех, математическое ожидание и стандартное отклонение Y, описываемой ФПВ gY(η), могут не существовать. Формулы (8) и (9) могут не обеспечить получение содержательных результатов. Однако плотность распределения вероятностей может быть использована для получения интервала охвата для выходной величины, т. к. G имеет смысл и может быть определена.
Интервал охвата
Пусть q = pM, если pM - целое число. В противном случае в качестве q можно выбрать целую часть (pM+1/2). Тогда [ylow, yhigh] – это 100p % интервал охвата для Y, где для любого r из ряда r = 1, …, (M – q), ylow= y(r) и yhigh= y(r+q).
Вероятностно симметричный 100p % интервал охвата можно получить, взяв
r = (M − q)/2, если (M − q)/2 — целое число, r=int[(M− q + 1)/2] - в противном случае. Для нахождения наименьшего 100p % интервала охвата задается определение r*, чтобы, что для r = 1, …, M – q выполнялось неравенство
y(r*+ q) – y(r*) ≤ y(r+q) – y(r).
Из-за стохастичности в ММК некоторые из протяженностей этих M-q интервалов могут быть короче, чем они будут в среднем, а некоторые - длиннее. Таким образом, при выборе такой наименьшей протяженности (аппроксимация) к наикратчайшему 100p % интервалу охвата имеет тенденцию быть, по крайней мере, меньше, чем при вычислении из GY(η), приводя к тому, что типичная вероятность охвата меньше, чем 100p %. Однако для больших М этим отличием можно пренебречь.Адаптивная реализация метода Монте-Карло
Получение установившихся числовых оценок статистических характеристик – суть адаптивной процедуры ММК. Численный результат считается установившимся, если соответствующее ему удвоенное стандартное отклонение станет меньше заданной точности вычисления стандартной неопределенности u(y).
Пусть ndig - число десятичных знаков, рассматриваемых как значащие в числовом представлении величины z. Предел погрешности вычисления δ, в таком случае, определяется следующим образом:
– выражают значение z в виде
z=c⸱10l , (10)
где с – это целое от ndig десятичных знаков,
l – целое число;
– δ определяют по формуле
(11)
Практический подход, включающий проведение последовательных применений ММК заключается в следующем:
– необходимо задать ndig;
б) задают M = max(J, 104), где J – наименьшее целое, которое больше или равно 100/(1–p);
в) задают h = 1 (счетчик итераций ММК);
г) проводят M испытаний методом Монте-Карло;
д) используют M полученных на выходе модели значений y1, …, yM, для вычислений h-й оценки y(h) величины Y, ее стандартной неопределенности u[y(h)], левой
и правой 
границ 100p % интервала охвата;
е) если h = 1, то увеличивают h на единицу и возвращаются к шагу г);
ж) вычисляют выборочное стандартное отклонение sy, связанное со средним оценок y(1), …, y(h) величины Y, задающееся как
(12)
(13)
и) аналогичным образом вычисляют выборочное стандартное отклонение для средних значений оценок u(y), 
к) используют далее все hM доступные значения модели для вычисления u(y);
л) определяют предел погрешности δ, связанной с u(y);
м) если хотя бы одно из значений 2sy, 2su(y), 
превышает δ, необходимо увеличить h на единицу и вернуться к шагу г);
н) если возврата к этапу г) не произошло, и значения всех вычисляемых значений можно считать установившимися, то на основе hM полученных значений выходной величины вычисляют y, u(y) и 100p% интервал охвата.
Обычно на этапе а) выбирается ndig=1 или ndig=2. На этапе ж) y может рассматриваться как реализация случайной величины со стандартным отклонением sy. В ситуациях, где интервал охвата не требуется, проверка стабилизации вычислений на этапе м) может быть основана вместо этого только на 2sy и 2su(y).
Альтернативный неадаптивный подход для симметричного интервала охвата при вероятности 95%, основанный на использовании статистик биноминального распределения, состоит в следующем. Выбирают М = 105 или М= 106. Формируется интервал [y(r), y(s)], где для М = 105, r= 2420 и s=97581, а для М = 106, r= 24747 и s= 975254. Этот интервал – 95 % статистический интервал охвата при уровне доверия 0,99, т. е. вероятность охвата будет не менее 95 % в 99 % применений ММК. Средняя вероятность охвата такого интервала будет (s − r)/(M + 1), что превышает 95 % на величину, которая становится меньше при возрастании M, а именно, 95,16 % для М = 105 и 95,05 % для М=106.
В результате применения адаптивной процедуры должны быть определены:
– оценка y величины Y,
– соответствующая стандартная неопределенность u(y);
– границы ylow и yhigh интервала охвата для Y, соответствующего заданной вероятности охвата.Преимущества метода Монте-Карло
ММК имеет следующие преимущества:
– возможна оценка любых статистических характеристик выходной величины Y, а не только стандартного отклонения;
– возможность поэтапного оценивания неопределенности, когда выход одной модели служит входом другой модели (допускается любое количество таких этапов);
– применимость для любых типов математических моделей результатов измерений, как линейных, так и нелинейных;
– улучшаются оценки выходной величины для нелинейных моделей;
– программное обеспечение метода доступно и относительно недорого;
– нет необходимости оценивать коэффициенты чувствительности.
Пример оценивания неопределённости измерений методом Монте-Карло
В качестве примера для оценки неопределённости измерений ММК был выбран метод определения количества удаляемого воздуха из помещения по ГОСТ 12.3.018.
Количество удаляемого воздуха из помещения L, м3/ч, определяют по формуле
L=Sк ∙ v∙ 3600, (14)
где Sк – площадь отверстия вентиляционного канала, м2;
v – значение скорости движения воздушного потока в вентиляционных каналах, м/с;
3600 – количество секунд в 1 ч.
Площадь отверстия вентиляционного канала круглого сечения Sк, м2, определяют по формуле
(15)
где dк – диаметр вентиляционного канала, м;
π = 3,14 – математическая константа.
Для упрощения расчётов объединим формулы (14), (15) и получим
Входные величины:
1. Диаметр вентиляционного канала:
– измеренное значение dк = 0,200 м;
– распределение вероятностей: прямоугольное (равномерное);
– интервал входной величины: ±Δ=±0,001 м;
2. Скорость движения воздушного потока в вентиляционных каналах
– измеренное значение v = 1,58 м;
– распределение вероятностей: нормальное;
– интервал входной величины: ±Δ=0,04 м/с.
С помощью программы Excel генерируем 2 массива случайных чисел (для входных величин dк и v) объемом M=106, подчиняющихся прямоугольному и нормальному закону распределения соответственно. Массивы данных входных величин и построенные гистограммы для каждого из них с целью проверки распределения сгенерированных данных представлены на рисунке 1.
Рисунок 1 – Массив данных входных величин
Примечание. Чтобы создать случайные данные выберите пункт «Анализ данных» – «Генерация случайного числа». Если на панели инструментов нет данного пункта, то его нужно добавить через окно «Параметры» – «Настройки» – «Пакет анализа», после чего во вкладке «Данные» появится пункт «Анализ данных».
Далее по каждому вектору вычисляется значение выходной величины, в результате получаем массив значений Y.
После чего массив значений выходной величины сортируется в неубывающем порядке и строится гистограмма, которая даст представление о распределении Y. На рисунке 2 представлена гистограмма распределения выходной величины.
Рисунок 2 – Гистограмма распределения выходной величины
Примечание. Для сортировки в неубывающем порядке выходной величины следует скопировать массив Y и вставить его в соседний столбец через параметр «Специальная вставка» – «Вставить значения».
По формулам (8), (9) вычисляем оценки параметров полученного распределения
Расчёт расширенной неопределённости осуществляется по формуле
(16)
Далее находим коэффициент охвата по формуле
(17)
Результат измерений:
Количество удаляемого воздуха из помещения составило L=178,79 ± 10,70 м3/ч.
