Обработка результатов статистических наблюдений с помощью толерантных интервалов
Толерантные интервалы
При определении понятия толерантного интервала, в качестве примера, будем рассматривать функцию плотности вероятности стандартного нормального распределения. На рисунках 1 и 2 изображены квантили распределения - значения, которые случайная величина не превышает с заданной вероятностью или
где
Иллюстрация к данному соотношению, где черный квадрат - соответствующий квантиль распределения:
Таким образом, можно сказать, что квантиль распределения представляет собой границу, ниже которой расположена доля
где
Можно переписать данное выражение в другой форме
где
Меньший и больший квантили назовем нижней и верхней границами соответствующего интервала, например
Перейдем непосредственно к определению понятия толерантного интервала. В Википедии, со ссылкой на соответствующую литературу, пишут:
Толерантный интервал является интервалом в выборочном пространстве наблюденных случайных величин. Он определяется достаточной статистикой на основе требования о том, чтобы при заданном доверительном уровне содержать вероятностную меру статистического распределения, не меньшую заданного уровня.[1]
Доверительный интервал определяется для некоторого параметра функции распределения и является интервалом в параметрическом пространстве. Он определяется достаточной статистикой на основе требования о том, чтобы вероятность того, что он содержит истинное значение неизвестного параметра была не меньше доверительного уровня.[1]
[1] Закс, 1975
Толерантные интервалы являются по своему смыслу доверительными интервалами, но не для параметров, явно участвующих в задании распределения случайных чисел, а для самих функций распределения. До настоящего момента мы еще не упоминали о доверительном уровне (или уровне доверия), который фигурирует в приведенных определениях. Поскольку речь идет о случайных величинах, то какие-либо утверждения необходимо делать в вероятностной форме. Итак, еще раз дадим определение толерантного интервала: толерантный интервал - это интервал, между границами которого с заданной вероятностью содержится доля распределения не меньшая заданной.
Стоит акцентированно отметить важную особенность толерантного интервала, а именно то, что толерантный интервал содержит в себе некоторую долю распределения, то есть долю всей генеральной совокупности, хотя строится он на основании выборочных значений.
Для дальнейшего использования введем обозначения. Долю совокупности (2), содержащуюся внутри интервала, переобозначим как:
и будем называть ее - доля накрытия. А вероятность, с которой данная доля (4) превышает некоторую заданную, будем называть уровнем доверия
где
Как уже говорилось, толерантный интервал строится на основании выборки некоторого размера
Из соотношения (6) видно, что соответствующие границы являются случайными числами. Поэтому, чтобы определить уровень доверия к доле накрытия
Извлечение выборки заданного размера из генеральной совокупности;
Нахождение оценок необходимых статистик (в общем случае - параметры сдвига и масштаба распределения);
Вычисление соответствующих границ толерантного интервала;
Вычисление доли генеральной совокупности
внутри интервала.
Находя долю интервалов, для которых выполняется условие
Уровень доверия - предел доли интервалов, накрывающих долю выбранной совокупности не менее заданной, при бесконечном повторении метода.[2]
[2] [ГОСТ Р 50779.10—2000, статья 2.61 ]
Другими словами:
где
Вообще говоря, толерантный интервал не обязательно должен быть центральным, то есть значения толерантных границ
Уровень доверия
Перейдем к вычислению уровня доверия. Разброс случайной величины удобно измерять в единицах параметра масштаба, в нашем случае - в единицах стандартного отклонения. Проиллюстрируем на примере:
В нашем случае
Реализуем процедуру расчета, а затем рассмотрим возможные варианты применения толерантных интервалов на практике. Для рабочих задач мне необходима большая статистика и малый шаг вычислительных сеток
Для начала определим некоторые функции, которые будут нам необходимы.
def cumulative_sum(x):
"""
Возвращает значение функции вероятности нормального распределения в точке x
"""
cumulative_sum = (1 + sc.special.erf(x / 2**0.5)) * 0.5
return cumulative_sum
В модуле stats
библиотеки scipy
также есть реализации функции вероятности для различных распределений.Теперь непосредственно реализуем процедуру вычисления толерантных границ. Вычисления будем производить для ряда значений коэффициентов
Зададим два одномерных массива, один из которых содержит значения толерантных коэффициентов
k1 = np.linspace(0.05, 10, 100, dtype=np.float64)
k2 = np.copy(k1)
Сгенерируем выборку случайных чисел с нормальным законом распределения размера sample_size
. Будем генерировать выборку не один раз, а некоторое число раз num_of_events
, равное числу статистических испытаний
# размер выборки n (или число наблюдений)
sample_size = 10
# число статистических испытаний N (сколько раз повторим процедуру)
num_of_events = 10000
Для каждого отдельного повтора по выборке найдем оценки математического ожидания и стандартного отклонения
sample = np.random.randn(sample_size, num_of_events)
sample_means = sample.mean(axis=0)
sample_stds = sample.std(axis=0)
sample_means
и sample_stds
имеют размерность равную num_of_events
, то есть содержат в себе значения среднего и выборочного стандартного отклонения для каждого из num_of_events
повторов.
В соответствии с выражением (6) вычислим толерантные границы (также для каждого повтора и каждой пары коэффициентов
mean_grid, k1_grid, k2_grid = np.meshgrid(sample_means, k1, k2)
std_grid, k1_grid, k2_grid = np.meshgrid(sample_stds, k1, k2)
l1 = mean_grid - k1_grid * std_grid
l2 = mean_grid + k2_grid * std_grid
Теперь вычислим значения доли накрытия
stat_coverages = cumulative_sum(l2) - cumulative_sum(l1)
Итоговая функция для расчета доли накрытия в скрытом блоке:
Функция для расчета доли накрытия
def calc_coverages(sample_size, num_of_events, tolerance_factor_lower, tolerance_factor_upper):
"""
Вычисляет долю накрытия генеральной совокупности с нормальным законом распределения
толерантными интервалами с заданными границами. Розыгрыш событий производится заданное
число раз для выборки заданного размера
Параметры:
1. sample_size - размер генерируемых выборок из нормального распределения
2. num_of_events - число статистических испытаний
3. tolerance_factor_lower - массив нижних толерантных коэффициентов
4. tolerance_factor_upper - массив верхних толерантных коэффициентов
Вычисление производится для каждой пары коэффициентов
Возвращаемые значения:
1. stat_coverages - статистически вычисленная доля накрытия генеральной совокупности (для каждого испытания)
толерантным интервалом с коэффициентами tolerance_factor_lower и tolerance_factor_upper
Размерность массива coverages [len(tolerance_factor_lower)][num_of_events][len(tolerance_factor_upper)]
"""
k1 = tolerance_factor_lower
k2 = tolerance_factor_upper
sample = np.random.randn(sample_size, num_of_events)
sample_means = sample.mean(axis=0)
sample_stds = sample.std(axis=0)
mean_grid, k1_grid, k2_grid = np.meshgrid(sample_means, k1, k2)
std_grid, k1_grid, k2_grid = np.meshgrid(sample_stds, k1, k2)
l1 = mean_grid - k1_grid * std_grid
l2 = mean_grid + k2_grid * std_grid
stat_coverages = cumulative_sum(l2) - cumulative_sum(l1)
return stat_coverages
Имея рассчитанные значения доли накрытия (4) мы найти число тех, которые не меньше заданного значения required_coverage
-
required_coverage = 0.95 # Любое значение в интервале [0, 1]
# ось вдоль которой расположены результаты каждого статистического испытания
event_ax = 1
successes = (stat_coverages > required_coverage).sum(axis=events_ax)
Последним шагом, в соответствии с выражением (7) вычислим уровень доверия как отношение числа успехов к общему числу статистических испытаний. Под успехом в нашем случая подразумевается событие, в котором рассчитанная доля накрытия интервалом оказалась не меньше заданной доли
confidence_level = successes / stat_coverages.shape[events_ax]
В результате мы получили зависимость уровня доверия к доле накрытия
И, например, для различных размеров выборки, на основании которой строится толерантный интервал:
Из полученных зависимостей следуют общие выводы:
Чем о большей части (доле
) генеральной совокупности необходимо сделать утверждение, тем шире должен быть соответствующий толерантный интервал (большие значения толерантных коэффициентов ) для обеспечения одного и того же уровня доверия . Еще раз отметим: уровень доверия - вероятность сделать ложный вывод о доле нашей совокупности;
Для примера - изменение положения точки, соответствующей уровню доверия
Чем меньше размер выборки, на основании которой мы делаем вывод относительной генеральной совокупности, тем шире границы соответствующего толерантного интервала (рисунок 10).
Обработка результатов наблюдений
Один из возможных вариантов применения - анализ надежности систем различного рода. Для примера, на рисунке 12 приведено изменение тока смещения операционного усилителя [3] от поглощенной дозы:
Оригинал:
На рисунке показано ухудшение тока смещения операционного усилителя PM 155 в зависимости от дозы. Было использовано 22 образца из трех различных производственных партий.Среднее значение и граница одностороннего толерантного интервала были рассчитаны для логнормального закона распределения,
и [3] [3] Christian Poivey. Radiation Hardness Assurance for Space Systems / NASA GSFC, 2002. «Обеспечение радиационной стойкости космических систем».
Смоделируем похожую ситуацию. Пусть необходимо сделать вывод о работоспособности партии (генеральной совокупности) изделий при воздействии определенного фактора, например, ионизирующего излучения космического пространства. Для испытаний на радиационную стойкость из всей партии извлекается выборка размера
где
Соответствующий код
loc = 0
scale = 0.01
voltage = lambda a, dose: a * dose
dose = np.linspace(0, 100, 100) # поглощенная доза
lcb = np.zeros_like(dose) - 1 # нижняя критическая граница ПКГ
ucb = np.zeros_like(dose) + 1 # верхняя критическая граница ПКГ
sample_size = 10
voltages = np.empty(shape=(sample_size, dose.shape[0]))
for i in range(sample_size):
a = sc.stats.norm.rvs(loc=loc, scale=scale)
voltages[i] = voltage(a, dose)
fig, ax = plt.subplots(figsize=(9, 7))
for v in voltages:
ax.plot(dose, v, lw=0.7)
average = voltages.mean(axis=0)
ax.plot(dose, lcb, lw=2, c='#00FF00', label='Критическое значение ПКГ')
ax.plot(dose, ucb, lw=2, c='#00FF00')
ax.plot(dose, average, lw=3, c='b', label='Среднее значение ПКГ')
Результат приведен на рисунке 13:
Зелеными линиями на рисунке обозначены верхняя и нижняя критическая граница параметра критерия годности, то есть значение напряжения, при достижении которого изделие считается вышедшим из строя (временная потеря работоспособности, необратимый отказ и так далее).
С помощью толерантных интервалов по результатам испытаний определим с
Выполним следующие действия:
Для каждого значения уровня воздействия
по результатам испытаний найдем толерантные коэффициенты
Для каждой пары коэффициентов
найдем значение (рисунок 8), для которого выполняется условие:
Итоговая функция для обработки результатов
def search_nearest(array, target):
return np.argmin(np.abs(array - target))
def get_coverages(confidence_arr, coverages_mesh, required_confidence, voltages, lcb, ucb, tolerance_mesh):
"""
Возвращает значение доли накрытия P, для которого уровень доверия максимально близок к заданному
Параметры:
1. confidence_arr - трехмерный массив значений уровня доверия, где первая ось соответствует
каждому значению доли накрытия; вторая и третья - толерантным коэффициентам k1 и k2
2. coverages_mesh - одномерный массив значений доли накрытия, для которых расчитан confidence_arr
3. required_confidence - требуемый уровень доверия к результату
4. voltages - массив, содержащий результаты испытаний (зависимости напряжения от дозы)
5. lcb - нижняя граница параметра критерия годности (массив)
6. ucb - верхняя граница параметра критерия годности (массив)
7. tolerance_mesh - одномерный массив значений толерантных коэффициентов (узлы расчетной сетки)
Возвращаемое значение:
1. coverages - массив значений доли накрытия (для каждого значения уровня воздействия),
для которых уровень доверия равен требуемому required_confidence
"""
mean = voltages.mean(axis=0)
std = voltages.std(axis=0)
k1_exp = (mean - lcb) / std
k2_exp = (ucb - mean) / std
k1_indices = np.array([search_nearest(tolerance_mesh, k) for k in k1_exp])
k2_indices = np.array([search_nearest(tolerance_mesh, k) for k in k2_exp])
coverages = np.empty_like(mean)
for i in range(k1_indices.shape[0]):
distance_init = 1e6
for j, arr in enumerate(confidence_arr):
confidence_level = arr[k1_indices[i]][k2_indices[i]]
distance = np.abs(confidence_level - required_confidence)
if distance <= distance_init:
distance_init = distance
coverages[i] = coverages_mesh[j]
return coverages
В результате получаем зависимость следующего вида:
Выбрав, например, максимально допустимую вероятность отказа
Заключение
В заключение хотелось бы сказать, что, конечно, область возможного применения толерантных интервалов не ограничивается одной лишь теорией надежности, и приведенный пример использования далеко не единственный, а просто один из многих. Как мне кажется, в принципе, данная идея может быть использована в любой области, затрагивающей вопросы статистической обработки данных. Например, возможно также осуществление обработки (корректировки поведения каких-либо систем на основе поступающих данных) не только в статике, но и в динамике и многое другого.
Мы рассмотрели только случай нормального распределения, но общая концепция может быть распространена и на другие законы распределения, в том числе и распределение Коши. В таком случае, необходимо будет выбрать соответствующие методы построения оценок параметров сдвига и масштаба. Также существуют методы построения так называемых «свободных от закона распределения» толерантных интервалов; подробнее можно почитать в книге:
[4] Кендалл М., Стюарт А. Статистические выводы и связи. М., Наука, Физматлит, Т. 2, 1973. — 899 с.
Вообще, все три тома хорошие, особенно могут быть полезны тем, кто занимается теоретической статистикой.
На этом, наверное, все. В данной статье я постарался поделиться с вами личным опытом знакомства с толерантными интервалами и понагляднее описать некоторые моменты, надеюсь, что кому-то приведенная информация сможет оказаться полезной. Буду рад вашему мнению, замечаниям или критике. Всем желаю здоровья и творческих успехов. Спасибо.