Эта публикация - продолжение базовой статьи Кригинг F-фактора или кормить, любить и никогда не покидать — «достаточно, но не чрезмерно», размещенной ранее на Хабре. В ней мы попробуем решить стандартную для технолога разработчика задачу: "Подбор времени выдержки продукта при температуре стерилизации" с опорой на ГОСТ. Статья адресована двум категориям читателей:

  • Технологам. Для них - простая пошаговая процедура решения задачи с использованием таблиц ГОСТа.

  • Специалистам в области контроля и диагностики. С ними мы попробуем в коде сделать то, что "зашито" в таблицах ГОСТа (ссылка на блокнот с  кодом в конце статьи).

План действий:

  1. Коротко о толерантных интервалах и ГОСТ Р 50779.29-2017 (ИСО 16269-6:2014).

  2. Детализация задачи "Подбор времени выдержки продукта при температуре стерилизации" и ее табличное решение для технологов.

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

Коротко о толерантных интервалах и ГОСТ Р 50779.29-2017 (ИСО 16269-6:2014)

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

Область применения: оценка средних эффектов, долей дефектов, коэффициентов регрессии, разностей средних и т. п. в  исследованиях и аналитике, где главная цель — выводы о параметрах популяции.

Толерантные интервалы (употребляемый синоним — «интервалы охвата/coverage intervals») характеризуют не параметры, а само распределение и определяются по выборке так, что с заданным уровнем доверия интервал накрывает долю совокупности не менее заданной.

Области применения: инженерные и регуляторные задачи, где нужно гарантировать долю изделий/наблюдений внутри границ с заданным доверием (размеры деталей, время нагрева, выход годных).

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

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

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

ГОСТ Р 50779.29-2017 (ИСО 16269-6:2014) это модифицированный перевод на русский язык международного стандарта ИСО 16269-6:2014.

По ГОСТу:
Статистический толерантный интервал - интервал, определяемый по выборке, относительно которого можно утверждать с уровнем доверия 1-\alpha, например 0,95, что он содержит не менее указанной доли p совокупности. Границы статистического толерантного интервала называют статистическими толерантными границами. Уровень доверия 1-\alpha - это вероятность того, что построенный установленным способом толерантный интервал содержит не менее заданной доли совокупности p.

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

Методика ГОСТа позволяет задав, четыре значения,:

  • два внешних требования (уровень доверия \gamma=1-\alpha и уровень покрытия p) и

  • два значения, относящихся к анализируемой выборке (порядковые статистики t_{(v)}и t_{(n-w+1)}, соответствующие наименьшему и наибольшему наблюдениям)

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

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

  • [0.9, 0.95, 0.99, 0.999] для уровня доверия -1-\alpha

  • [0.9, 0.95, 0.99] для уровня покрытия -p.

Границы интервала оценивают в узлах сетки по таблицам из Приложения Е (обязательное). Таблицы Е.1 и Е.2. Страницы 38, 39.

Детализация задачи "Подбор времени выдержки продукта при температуре стерилизации" и ее табличное решение

В базовой статье мы писали о том, что процесс автоклавирования может управляться:

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

  • по целевому значению показателя летальности (F‑фактор, таргет).

Во втором случае в продукт помещается датчик и на каждом временном шаге по формуле Флауменбаума рассчитывается показатель летальности. Когда его величина достигает целевого значения этап выдержки при температуре стерилизации завершают и переходят к этапу охлаждения.

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

Рисунок 1. Время выдержки при температуре стерилизации до достижения целевого значения показателя летальности.
Рисунок 1. Время выдержки при температуре стерилизации до достижения целевого значения показателя летальности.

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

Различия времен выдержки для разных датчиков связаны с неравномерностью теплового поля и, как следствие, поля F-фактора (их изображения на заглавной картинке). Основные причины различий описаны в базовой статье.

Пример плана эксперимента (создается технологом):

  • целевое значение показателя летальности - 20 единиц

  • суммарное количество варок с датчиками - 80

  • количество датчиков - 5

  • схема размещения датчиков в автоклаве - должна быть единой для экспериментальных варок и учитывать размещение датчиков в самых "холодных точках" (для современных горизонтальных автоклавов в режимах с душированием это могут быть середины корзин)

  • количество автоклавов на которых будет обрабатываться однотипный продукт - 4.

Проведению эксперимента предшествуют подготовительные работы:

  • валидация участвующих автоклавов

  • валидация датчиков

  • планирование графика экспериментальных варок с учетом гарантированной доступности электрических и тепловых мощностей в интервалах выдержки.  

Цель подготовительных работ - минимизировать влияние неслучайных факторов на экспериментальные варки.

Реализовав план технолог сохраняет результаты измерений - 400 значений интервала времени (80 * 5) и агрегирует их,  выбирая по каждой варке показания "самого холодного" датчика, что соответствует самому продолжительному времени выдержки. Это наиболее консервативный подход, обычно применяемый в предметной области, где цена ошибки недовара  выше цены перевара. В итоге получается 80 значений интервалов выдержки.

Для статьи генерируются синтетические данные с учетом дискретности измерений (ссылка на блокнот Colab с кодом генерации в конце статьи). Код записывает созданные данные в Excel файл synthetic_heating_30s.xlsx и сохраняет его в сессионном хранилище. Лист Sensors содержит сырые данные. Лист Aggregates - агрегированные. В листе Y_sorted агрегированные данные отсортированы по возрастанию времени выдержки и пронумерованы от 1 до 80 (суммарное количество варок с датчиками).

Рисунок 2. Гистограммы времени выдержки при температуре стерилизации.
Рисунок 2. Гистограммы времени выдержки при температуре стерилизации.

Вывод по результатам визуализации:

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

Процедура табличного решения для технологов

Определимся с прикладной задачей: необходимо найти такое значение времени выдержки при температуре стерилизации, чтобы с одной стороны обеспечивать достижение заданного таргета по фактору летальности, а с другой — не перегревать продукт чрезмерно, поскольку это приводит к деструкции витаминов, повышению кислотного числа и содержания амино-аммиачного азота.  Таргет по F‑фактору 20 единиц выбран с двукратным технологическим запасом, поэтому целесообразно уточнить правую границу (самые длительные выдержки) с точки зрения ее обоснованности.

Выборочное распределение агрегированных времён (Aggregates) показывает, что большая часть «гарантирующих» интервалов лежит в области примерно до  38:30 минут, а наблюдения в районе 39:00-39:30 выглядят как относительно редкие выбросы.  


С учётом того, что целевое значение F‑фактора задано с запасом, можно использовать нормализованный подход: по процедурам ГОСТ определить односторонний толерантный интервал и по его верхней границе выбрать предельную длительность выдержки, сохраняя приемлемый риск недостижения таргета в контролируемом  числе варок при заданном уровне доверия.

Замечание.   На практике, оценить достаточность запаса по таргету для отсекаемых по границе редких варок можно по автоклавным лог-журналам эксперимента. В них содержатся фактические значения F-фактора по всем датчикам, в том числе и на момент отсечения.

Первый шаг

С учётом описанной ранее схемы эксперимента результаты измерений заносятся в лист Sensors.  Затем для каждой варки выполняется агрегация по «самой холодной» точке, и эти значения записываются в лист Aggregates.

Далее агрегированные значения сортируются по возрастанию, и каждому присваивается ранг от 1 до 80. Создается лист Y_sorted с тремя колонками: Y_sorted_sec (время в секундах), Y_sorted_mmss (формат ММ:СС) и rank (порядковый номер).  Таким образом, Y_sorted — это вариационный ряд\{t_{(1)},\dots,t_{(n)}\} сn=80, гдеt_{(j)}j-я по величине наблюдаемая выдержка для «самого холодного» датчика варки.

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

Это даёт основания интерпретировать упорядоченные значения (вариационный ряд) как реализацию порядковых статистик из одного непрерывного распределения неизвестного вида, для которого с учетом сформулированной прикладной задачи применим верхний односторонний непараметрический метод оценивания границы толерантного интервала по ГОСТ Р 50779.29‑2017 (ИСО 16269‑6:2014).

Второй шаг

Задаём входные требования к толерантному интервалу:

  • Уровень доверия1-\alpha из перечня [0.9, 0.95, 0.99, 0.999] — вероятность того, что построенный интервал будет накрывать не менее долиp генеральной совокупности времён выдержки.

  • Уровень покрытияp из перечня [0.9, 0.95, 0.99] — минимальная доля всех потенциальных варок (при неизменной схеме процесса), выдержки которых должны попадать в интервал.

Чем выше значенияp и1-\alpha, тем более консервативным становится толерантный интервал: требуется больший объём выборки и границы сдвигаются дальше от центра выборочного распределения (для двусторонней постановки, при верхней односторонней - вправо, при нижней - влево), приближаясь к экстремальным наблюдениям.  

Для примера возьмём сочетание (p,1-\alpha)=(0.9, 0.9): оно означает, что требуется интервал, который будет накрывать не менее 90 % всех потенциальных варок,  с вероятностью не ниже 0.9.  Это наименее жесткие требования по сетке стандарта.

Третий шаг

Для непараметрических толерантных интервалов используются таблицы, связывающие объём выборкиn, заданныеp и1-\alpha с  параметрамиv иw, определяющими номера порядковых статистик, задающих границы интервала.  

Рисунок 3. Таблица Е.1 (фрагмент).
Рисунок 3. Таблица Е.1 (фрагмент).
  • Открываем  ГОСТ Р 50779.29-2017 и находим Приложение Е (обязательное), таблицы Е.1 и Е.2.

  • В таблице Е.1 выбираем столбец, соответствующий выбранным уровням (p,1-\alpha)=(0.9, 0.9).

  • В этом столбце находим наибольшее табличное значение объёма выборки, не превосходящее фактическоеn=80; в нашем примере это значение78, ближайшее снизу к 80 (зеленая рамка на рисунке 3).

  • В строке, соответствующей этому значению78, считываем табличную величинуv+w (в ГОСТ она задаётся одной колонкой) — в примере это5.

Для верхнего одностороннего интервала по стандарту параметрv (число «отброшенных» наблюдений слева) принимается равным нулю:v=0.  Это отражает тот факт, что нижняя граница интервала не фиксируется и  весь контроль ведется по  правой, а табличная величинаv+w целиком относится к правому хвосту:w=5.

Граница верхнего одностороннего интервала задаётся через порядковую статистикуt_{(n-w+1)}, поэтому граничный рангk равен: k = n - w + 1.

Для нашего примера: k = 80 - 5 + 1 = 76.

Обращаемся к листу Y_sorted и по рангуk=76 берём значение Y_sorted_mmss, в примере оно соответствует выдержке в 38:00 минут.  

Это значение является оценкой верхней односторонней статистической толерантной границы при (p,1-\alpha)=(0.9 ,0.9) в непараметрической постановке (распределение непрерывно, вид неизвестен).

Рисунок 4. Лист Y_sorted.
Рисунок 4. Лист Y_sorted.

По аналогии таблицы Е.1 и Е.2 позволяют получить границы и для других комбинацийp и1-\alpha (желтые рамки на рисунке 3), что даёт возможность сравнить несколько альтернативных границ для времени выдержки в зависимости от степени консервативности требований.

Вариант интерпретации результата (для технолога)

Полученный в соответствии с методикой ГОСТ Р 50779.29-2017 интервал[t_{min}, t_{(76)}] с верхней границейt_{(76)}=38{:}00 мин — это односторонний непараметрический толерантный интервал для распределения времён выдержки по «самым холодным» датчикам варок.  По определению толерантного интервала это означает, что при фиксированной в эксперименте схеме технологического процесса и с уровнем доверия не ниже 1-\alpha=0.9 можно утверждать: не менее доли p=0.9 все варки будут иметь время достижения таргета по F‑фактору в наиболее холодной точке не больше 38:00 минут.

При этом единичные варки с временем достижения таргета в районе 38:30-39:30 (всего 4 значения - по рисунку 4) трактуются как редкие реализации верхнего хвоста распределения времени выдержки. С учётом того, что целевой F‑фактор задан с двукратным запасом, такой уровень «остаточного» риска можно считать приемлемым (по данным автоклавных лог-журналов эксперимента фактическое значение F-фактора по этим четырем измерениям на момент выдержки 38:00 минут составляло не менее 16 единиц, что соответствует критерию промышленной стерильности для выпускаемой продукции).

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

Обращение к читателю: насколько понятным и насколько полезным показался табличный подраздел?

Не табличные решения для оценивания границ непараметрических толерантных интервалов

Кратко о методе Уилкса

В ГОСТ Р 50779.29-2017 (ИСО 16269-6:2014) имя Сэмюэла Уилкса прямо не упоминается, но приложение G и таблицы E.1–E.2 основаны на   подходе, предложенном им в  статье Determination of Sample Sizes for Setting Tolerance Limits, опубликованной в 1941 году. Выражение G.1 из приложения G — частный случай его подхода для центрального интервала с поиском точек в хвостах распределений.

Метод Уилкса предлагает строить непараметрические толерантные интервалы на основе биномиальной модели для числа наблюдений, попадающих внутрь или наружу интервала, без каких‑либо предположений о виде распределения.  Идея в том, чтобы подобрать размер выборки и ранги порядковых статистик так, чтобы вероятность накрытия не менее чем p доли  совокупности была не ниже\gamma. Условие записывается через биномиальное неравенство для числа точек в хвостах.  

В прикладных задачах это даёт простой критерий: при заданныхp,\gamma и объёме выборкиn можно либо подобрать рангиv,w (центральный интервал), либо ранг одной границы для одностороннего интервала, гарантируя заданный уровень покрытия при произвольном  законе распределения.  Такой подход хорошо подходит для задач, где важно гарантировать долю «хороших/плохих» объектов.

Классического Уилкса можно найти здесь, на странице 342. А по этой ссылке доступен пример современной модификации метода, основанной на интерполяции порядковых статистик, с целью повышения точности оценивания границ толерантных интервалов.

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

1. Воспроизведение в коде решений примеров из подраздела 5.7 ГОСТ

Подраздел содержит четыре примера задач и их решений.

Рисунок 5. Примеры задач с решениями из подраздела 5.7 ГОСТ Р 50779.29-2017.
Рисунок 5. Примеры задач с решениями из подраздела 5.7 ГОСТ Р 50779.29-2017.

Численная реализация биномиального условия  воспроизводит эти примеры.  Для этого используется функция min_n_np_ti(p, gamma, v, w), которая решает обратную задачу относительно выражения G.1: находит минимальныйn и фактические уровни доверия1-\alpha_{\text{ach}} при заданныхp,\gamma,v,w. Полученные результаты  совпадают решениями примеров из ГОСТ (ссылка на блокнот Colab с кодом в конце статьи).

  • Вход:  

    • p — требуемая доля покрытия

    •  gamma — требуемый уровень доверия

    • v — число наблюдений, отбрасываемых слева

    •  w — число наблюдений, отбрасываемых справа.

  • Внутренние обозначения:

    • q = 1-p — доля совокупности вне интервала

    • r=v+w — суммарное число допустимых точек в двух хвостах

    •  Y\sim\text{Bin}(n,q) — число выборочных точек, выходящих за интервал.

  • Алгоритм:

    • подбирается минимальноеn, при котором выполняется неравенство:

\Pr\{Y \le r-1\} = F_{\text{Bin}}(r-1;n,q) \le \alpha = 1-\gamma,

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

\Pr\{Y \ge r\}\ge\gamma.

  • Реализация:

    •  специальный случай: аналитическое решение для r=1

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

  • Функция возвращает минимальный объём выборки, такой, чтобы центральный интервал междуX_{(v+1)} иX_{(n-w)} обеспечивал накрытие доли не менееp с уровнем доверия не ниже\gamma при произвольном непрерывном распределении.

Рисунок 6. Собственные решения для примеров из подраздела 5.7.
Рисунок 6. Собственные решения для примеров из подраздела 5.7.

2. Гистограмма с границами Exact и GOST

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

  • Exact - точное вычисление  ранга верхней односторонней границы толерантного интервала по биномиальному распределению  точек внутри интервала.

  • GOST - реализация алгоритма, обратного выражению G.1 ГОСТ через суммарное число точек, отсекаемых снаружи интервала - в хвостах биномиального распределения.

По итогам  строится гистограмма Y_sorted с шагом 30 секунд и оценкой плотности (KDE), на которую наносятся две вертикали, соответствующие верхним односторонним границам, посчитанным для каждого из подходов:t^*_{\text{Exact}} = Y_{(k_{\text{Exact}})} и t^*_{\text{GOST}} = Y_{(k_{\text{GOST}})}.

  • В коде граница Exact вычисляется напрямую с применением функции binom.ppf(gamma, n, p).

  • Граница  GOST вычисляется через адаптированную функцию  min_n_np_ti(p, gamma, v, w) на основе той, что была использована для расчета примеров.

Рисунок 7. Два подхода к оцениванию верхней односторонней границы толерантных интервалов.
Рисунок 7. Два подхода к оцениванию верхней односторонней границы толерантных интервалов.

Оба подхода для (p,\gamma)=(0.9, 0.9) и n=80  дают одинаковый результат, совпадающий с табличным:

  • рангk=76 и

  • времяt^*=38{:}00.  

Это означает, что выбранное значение выдержки одновременно удовлетворяет обоим критериям, а код верно реализует методику ГОСТ (по крайней мере, на этих данных). Дополнительно рассчитываются достигнутые уровни(p,\gamma) для каждого из подходов.

3. Диаграмма «внутри — снаружи»

Рисунок 8.1. Поиск решения для границы двумя способами на ступенчатых CDF.
Рисунок 8.1. Поиск решения для границы двумя способами на ступенчатых CDF.

Диаграмма содержит две ступенчатые CDF:

  • F_p(t)=\Pr\{T\le t\} для T\sim\text{Bin}(n,p) (распределение точек «внутри» интервала) и

  • F_q(t)=\Pr\{Y\le t\} для Y\sim\text{Bin}(n,q),q=1-p (распределение точек «снаружи» относительно интервала).

На графике вертикальной линией отмечается соответствующий первому условию биномиальный квантиль: m_{\text{ex}} (и связанный с ним ранг k_{\text{Exact}}).
А также соответствующее второму условию количество точек снаружи: r^*=v+w (и связанный с ним ранг k_{\text{GOST}}).
При этом горизонтали на уровнях \gamma и \alpha=1-\gamma показывают, какие именно биномиальные неравенства реализуют оба подхода.  

  • Exact решает задачу «хотим, чтобы не болееk-1  точек было внутри интервала», а

  • GOST — «хотим, чтобы не менееr=v+w точек оказалось в хвостах».

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

В медицинской статистике и в задачах анализа надежности при решении биномиальных задач широко используется функция выживаемости (survival function). На рисунке 8.2 приведена диаграмма, демонстрирующая двоичный поиск решения на кривой хвостовой вероятности функции выживаемости.

Рисунок 8.2. Двоичный поиск верхней односторонней границы на кривой хвостовой вероятности.
Рисунок 8.2. Двоичный поиск верхней односторонней границы на кривой хвостовой вероятности.

Подписи у точек означают:

  • номер шага

  • анализируемый на текущем шаге ранг

  • значение хвостовой вероятности в этой точке

  • текущий диапазон анализа (в рангах).

    С помощью функции binom.sf() и двоичного поиска решение находится за шесть шагов.

4. Лестница рангов

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

Рисунок 9. Лестница рангов.
Рисунок 9. Лестница рангов.

На ней видно «насколько глубоко в хвост» уходит граница, какие и сколько единичных экстремальных наблюдений остается правее t^*, каково соотношение между объёмом выборки и количеством отсеченных хвостовых точекw.  Приk=76 правая часть вариационного ряда (доk=80) содержит всего четыре наблюдения, что согласуется:

  • с таблицей ГОСТ: w_{\text{one}} = v+w = 5, при v=0 и

  • мягкими требованиями:  (p,\gamma)=(0.9, 0.9).

5 Проверка устойчивости границ толерантного интервала

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

Существует, по крайней мере, два подхода для получения "возмущенных данных":

Первый подход — это параметрическое моделирование или параметрический бутстрэп (parametric model-based bootstrap):

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

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

  • далее для каждой псевдовыборки в соответствии с формулой G.1 ГОСТ рассчитывается искомая статистика - граница толерантного интервала.

    Таким образом строится модельно‑имитационное (не по исходным данным, а по подобранной и параметризованной на них модели) распределениеt^*.

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

Второй подход — это непараметрическое бутстрэп‑моделирование или непараметрический бутстрэп (nonparametric data-based bootstrap):

  • по исходной выборке Y_sorted многократно формируются псевдовыборки с возвращением непосредственно из исходных данных (сэмплинг с возвратом)

  • для каждой бутстрэп‑выборки определяется значение оценкиt^* тем же методическим правилом (например, формулой G.1 ГОСТ).

    Итоговое распределениеt^* строится исключительно на основе структуры внутренних взаимосвязей элементов исходной выборки, без каких‑либо предположений о виде распределения.

Достоинство этого подхода — независимость от модельных допущений и сохранение всей «экспериментальной» информации. Основной недостаток -  при малых объёмах выборки его возможности в описании  сложных внутренних зависимостей в данных ограничены.

Важно: оба метода формируют «эмпирическое» распределение итоговой статистики, однако основаны на принципиально разных предпосылках:

  • параметрический — на доверии к выбранной модели

  • непараметрический — исключительно на данных и их внутренней структуре.

В настоящей работе акцент вынуждено делается на непараметрическом бутстрэпе по причине особенностей эмпирического выборочного распределения времени выдержки (рисунок 2).

Рисунок 10. Результаты бутстрэп-моделирования оценки t*
Рисунок 10. Результаты бутстрэп-моделирования оценки t*

На рисунке 10 представлены результаты непараметрического бутстрэп-моделирования на 5 000 псевдовыборках.

Для примера Y_sorted, (p,\gamma)=(0.9 ,0.9) GOST‑процедура по формуле G.1 даёт оценку верхней односторонней границы  выдержкиt^*=38{:}00 мин. Бутстрэп‑моделирование показывает, что при повторении эксперимента 5 000 раз оценкаt^* распределяется по шести квантам времени:

С вероятностью ≈95 %  t^*попадает в диапазон 37:30–38:30 мин, медиана совпадает, а среднее значение близко к оцененной границеt^*=38{:}00 мин и равно38{:}02 мин.

При этом мода дискретного распределения приходится на 38:30 мин. Таким образом, при избранном “минимально жёстком” наборе требований (p,\gamma)=(0.9 ,0.9) оценка толерантной границы лежит между двумя соседними высокочастотными квантами (38:00 и 38:30), а значит является чувствительной к выборочным вариациям.

Вывод по результатам анализа устойчивости:

Необходимо перейти к более строгому набору требований (p,\gamma) и повторить процедуру оценивания по ГОСТ.

Для(p,\gamma)=(0.9 ,0.95) и кодовое решение и таблица E.1 ГОСТ дают:

k=77 =>t^*=38{:}30 мин,

которое совпадает с результатами бутстрэп-моделирования. И позволяет сделать новый вывод о том, что нормативное по ГОСТ решение t^*=38{:}30 мин при требованиях (p,\gamma)=(0.9 ,0.95) устойчиво к выборочным вариациям. Устойчивость оценена непараметрическим бутстрэп-моделированием на уровне 95% доверительного интервала.

Заключительные соображения

  • Формула Уилкса, лежащая в основе ГОСТ Р 50779.29‑2017, естественным образом выражается через биномиальное распределение и может быть реализована с помощью высокопроизводительной библиотеки scipy.stats, что позволяет обходиться без таблиц, выходить за рамки сетки стандарта по(p,\gamma,n)и расчитывать фактические достигнутые значения по уровням покрытия и доверия.
    Это не отменяет правила - пользоваться исключительно ГОСТ в нормативно чувствительных случаях.

  • Exact‑подход и ГОСТ‑подход используют одну и ту же биномиальную модель, но по‑разному формулируют критерий (через число точек внутри интервала или через суммарное число точек в хвостах), и в рассмотренном примере оба критерия привели к одинаковой оценке верхней границе 38:00 мин.

  • Анализ устойчивости практически важный этап. Его роль состоит не в том, чтобы ставить под сомнения результаты ГОСТа. С его помощью мы, в определенных пределах, можем правильнее сформулировать требования (как в статье), либо понять, что нам не хватает объема выборки.

  • Надеемся, что визуализации и код помогли  немного "заглянуть под капот" метода.

Ссылка на блокнот Colab

Открыть блокнот Colab

Сергей Сюр

Комментарии (0)