CC BY
92
Аналитическая химия
УДК 543.25+378.02
МОДЕЛИРОВАНИЕ ДИСПЕРСИОННОГО АНАЛИЗА В ЭЛЕКТРОННЫХ ТАБЛИЦАХ И ОЦЕНКА МОЩНОСТИ КРИТЕРИЯ ПРИ РАЗЛОЖЕНИИ ОШИБОК НА СОСТАВЛЯЮЩИЕ
В.И. Голованов
Предложены принципы стохастического моделирования матриц плана для дисперсионного анализа. Обсуждаются перспективы использования стохастического моделирования при решении учебных и исследовательских задач химической метрологии.
Ключевые слова: дисперсионный анализ, метод Монте-Карло, статистический тест, мощность критерия, компьютерные технологии в аналитической химии, химическая метрология.
Сравнительно недавно опубликован национальный стандарт Российской Федерации [1], в котором впервые вводится понятие и способы применения мощности критерия (теста) при проверке статистических гипотез о средних и дисперсиях. Эти методики предназначены, в первую очередь, для статистической оценки качества. Показатели качества ряда продуктов во многом определяются показателями качества химического анализа [2, с. 107]. Однако в учебной и методической литературе по химической метрологии [2-4] раскрытию понятия и количественной оценке мощности критерия уделяется, с нашей точки зрения, недостаточно внимания. Так в [4, гл. 12.1.2] при раскрытии темы по статистическим тестам в качестве основных этапов теста называют: 1) выбор подходящего критерия, 2) выбор уровня значимости, 3) формулирование нуль-гипотезы (Н0) и альтернативной гипотезы (На), 4) вычисление тестовой статистики и ее сравнение с критическим значением. В этой последовательности не хватает заключительного этапа - не хватает оценки мощности теста, характеризующего надежность статистического вывода вероятностью непринятия Н0, когда верна На. Как считают авторы [5]: «Мощность критерия позволяет углубить выводы, полученные при испытании гипотез». Применительно к алгоритмам дисперсионного анализа, которые описаны в [6], аналогичные рекомендации высказывают авторы работ [7]. Критикуя ГОСТ Р ИСО 5725-2-2002, авторы [7] приходят в выводу: «Поэтому чтобы стандарты и информационные технологии «заработали» корректно, необходимо, прежде всего, преодолеть статистическую безграмотность специалистов, практиков».
Компьютерные технологии обучения при их реализации в электронных таблицах, содержащих богатый набор средств тестовой статистики и стохастического моделирования, должны эффективно использоваться при преподавании вузовских курсов химической метрологии (и не только...). Это, в свою очередь, ставит задачи по разработке нового инструментария для моделирования различных ситуаций (сценариев) при химическом анализе и обработке его результатов. Данная работа является продолжением серии публикаций [8, 9] по этой тематике.
Принципы и технологию построения стохастических моделей рассмотрим на примере построения матрицы плана однофакторного дисперсионного анализа. План включает измерения результирующего признака Ху (концентрации) в аналитических сериях у = 1, п, в каждой из которых выполняют / = 1, т параллельных измерений. Модель строим исходя из основных предположений дисперсионного анализа: выборки получены из нормальных совокупностей с математическими ожиданиями М(Х!).М(Хп) и одинаковой дисперсией а2. Эти положения записывают уравнением:
Ху = М( X) + еу + е1}. (1)
Случайную величину в, называют факторным эффектом, в, - остаточным эффектом. С учетом принципа, который называют методом обратных функций Монте-Карло [10, с. 321] уравнение (1) перепишем в виде:
X,, = М(Х) + N- (г) • а* + N- (г) • а или
1*1 (2)
X, = N (г, 0, а ) + ^\г, М(X), аД
где г - равномерно распределенные на интервале (0, 1) случайные числа. Уравнения (2) являются основой для генерирования выборочных матриц плана однофакторного дисперсионного анализа. Заметим, что метод обратной функции является универсальным методом генерации псевдослучайных чисел по любому закону распределения.
Для наглядности изложения нашего метода возьмем в качестве примера результаты, полученные в [2, с. 140]. Решается задача о вкладе погрешности пробоотбора антифрикционного сплава в общую погрешность определения массовой концентрации сурьмы. Для шести кернов сплава выполняют определение сурьмы (% 8Ь) в двойной повторности для каждого керна:
14,72
15,05
15,51
15,23
14,60
14,35
15,10
15,23
14,70
14,95
14,74
14,50
В результате дисперсионного анализа найдены характеристики разброса между сериями,
внутри серий и дисперсию пробоотбора: s 1 = s 2 + ns = 0,2219; s 2 = 0,0326; s = (s i - s 2)/n = = 0,0942 соответственно.
Воспользуемся результатом [2] для оценки параметров уравнений (2). Среднее значение результата анализа примем в качестве оценки M(A’) = 14,89. В качестве оценок а1 и а примем s1 и s . При вычислениях мощности критерия будем использовать параметр X2 = а21/а22 = 6,81.
На рис. 1 показана форма листа электронной таблицы, которая содержит массив ячеек, связанных формулами. Расчетные формулы в нотациях MS Excel приведены в табл. 1. Особенностью вычислений является то, что, во-первых, факторные эффекты генерируют независимые датчики вида НОРМОБР(СЛЧИС();0;$С$5а0.5), число которых равно числу уровней фактора, а остаточные эффекты генерируют датчики вида НОРМОБР(СЛЧИС();$С$3;$С$4А0.5)+В$8, число которых равно размерности матрицы плана. Во-вторых, факторные эффекты передаются по ссылке из массива погрешностей B8:G8 в ячейки матрицы плана B10:G11. В данном примере гипотеза об однородности дисперсий не проверяется, поскольку ее справедливость изначально заложена в модели. Правильность результатов дисперсионного анализа с использованием формы на рис. 1 проверяли расчетами с использованием утилиты Excel «Однофакторный дисперсионный анализ» из надстройки «Анализ данных». Результаты совершенно идентичны.
в
Стохастическая модель дисперсионного анализа Параметры модели
14 89 0 0326 00942
x(Sb). масс % о22(анализа)= о2’ (пробоотбора(=
Генерирование матрицы плана однофакторного дисперсионного анализа
1) погрешностей пробоотбора шести кернов баббита
0.0963 -0.1378 0.1827 0.5947 -0.2541 -0.2923
2) погрешностей анализа и их смешивания с погрешностями пробоотбора
14.83 15.07 14.82 15.39 14.62 14.13
14.87 14.70 15.02 15.40 14.56 14.59
Разложение ошибок (дисп анализ)
СрЗнач 14.85 14.89 14.92 15.40 14.59 14.36
Дисперсия 0.0006 0.0699 0.0204 0.0001 0.0019 0.1081
5^= 0.2434 F= 7,27
Мощность критерия
s22= 0.03 35
0.1049
>болыие> 4.39 = F(a.,f1,f2)
ЧислИсп= 498 ЧислУдач= 160
Мощность
68%
Рис. 1. Форма листа Excel для моделирования мощности критерия
Таблица 1
Формулы в ячейках на рис. 1
Ячейки Функции
B8:G8 НОРМОБР(СЛЧИС();0;$С$5А0.5)
B10:G10 НОРМОБР(СЛЧИС();$С$3;$С$4А0.5)+В$8
B11:G11 НОРМОБР(СЛЧИС();$С$3;$С$4А0.5)+В$8
B13:G13 СРЗНАЧ(В10:В11)
B14:G14 ДИСП(В10:В11)
B15 СЧёТ(В10:В11)*ДИШ(В13^13)
D15 СРЗНАЧ(В14^14)
G15 ЕСЛИ(В16^16;(В15^15)/СЧёТ(В13:В14);0)
B16 B15/D15
C16 ЕСЛИ(В16^16;"<меньше<";">больше>")
D16 FРАСПОБР(0.05; СЧЁТ(В 10:G10)-1; СЧЁТ(В10^10)*(СЧёТ(В10:В11)-1))
B18 (B18+1)
B19 ЕСЛИ(В16^16;В19;В19+1)
E19 1-B19/(B18-1)
Удобство использования нашего алгоритма состоит, во-первых, в возможности многократного прогона процедур анализа и, во-вторых, в автоматизации стадии статистического вывода. (Обратим внимание на блок ячеек A16:G16.) Предложенный нами алгоритм обладает достаточной гибкостью для того, чтобы перестраивать форму под задачи с матрицами плана большего или меньшего размера. В учебных целях можно легко изменять сценарий. Например, факторными эффектами могут быть погрешности при анализе в условиях внутрилабораторной прецизионности или межлабораторная ошибка.
Многие особенности дисперсионного анализа можно исследовать, не прибегая к эксперименту, например, исследовать его основные предположения. Апостериорное использование рассмотренного подхода может способствовать углублению обоснованности статистических выводов из лабораторных экспериментов.
Дополнительно к традиционному представлению результатов дисперсионного анализа разработанная нами вычислительная форма позволяет исследовать, с использованием техники Монте-Карло, эффективность разложения погрешностей. При многократных прогонах программы и подсчете процента случаев успешной дискриминации ошибки пробоотбора получаем оценку мощности (эффективности) теста Фишера. Мощность принято обозначать как дополнение к ошибке второго рода: 1 - р.
Подсчет числа случаев, когда принимается альтернатива На: а21 > а22 против Н0: а21 = а22, осуществляется с помощью сконструированного нами счетчика, который работает в итерационном режиме и изменяет показания при пересчете листа. Причем встроенный в Excel датчик случайных чисел Л(г|0,1) также обновляет значения при пересчете листа. Для того чтобы создать счетчик, пользователю следует установить на вкладке «Вычисления» переключатель «вручную», активировать флаг «итерации» и назначить предельное число итераций равное 1. Итерационный процесс осуществляется с использованием формул, записанных в ячейках В18; В19; Е19 (см. табл. 1). Для запуска счетчика следует нажимать/удерживать клавишу F9. Чтобы после завершения серии испытаний привести счетчик к исходному состоянию, следует установить фокус на ячейке В18 и двойным щелчком войти в формулу. Последующий ввод формулы клавишей «Enter» приведет к отображению числа 1. Аналогичная процедура установит в ячейке В19 число 0/1, в зависимости от сообщения «больше/меньше» в ячейке С16.
Вероятность отвергнуть Н0, когда верна альтернативная гипотеза На, вычисляется как относительная частота события «ЧислУдач» по формуле, записанной в ячейке Е19. На рис. 1 видим, что потребовалось около 500 испытаний модели, при которых показатель мощности в ячейке Е19 перестает изменяться. Для получения такого результата требуется не более 30 секунд.
Результаты испытания эффективности F-теста с данными [2] приведены на рис. 1 в ячейке Е19. Заметим, что специалисты фирмы StatSoft [11] считают приемлемой мощность теста не меньшую чем 80 %. Мощность теста следует усилить. Есть четыре основные стратегии повыше-
ния мощности: выбор более высокого уровня ошибки первого рода а, формулирование направленных гипотез, увеличение объема выборки и усиление эффекта. Для нашего случая наиболее важным представляется увеличение размерности матрицы плана, а также осуществление мероприятий по усилению изучаемого эффекта. Очевидно, что эффект будет усилен при повышении точности измерения результирующего признака, т. е. при уменьшении величины а2. В вычислительной среде нашего алгоритма доступным остается только увеличение размерности матрицы плана. Априори не ясно следует ли увеличить число параллельных опытов т или более эффективным будет увеличение числа проб п? (Хороший вопрос для проблемного обучения.) Однако имеется иной способ решения задачи, который следует из теории статистических тестов. Мощность критерия можно вычислить [1, 5], не используя технику Монте-Карло.
Исходя из принципов, которые обосновываются в руководствах по математической статистике, для нашей задачи запишем:
Ра (/1, /2)
^1-Р (fl, /2) =-
X2
(3)
То есть квантиль распределения Фишера - Снедекора для вероятности 1 - в можно вычислить с использованием таблиц Б-распределения. Поскольку уровень значимости в задачах химического анализа обычно принимают равным 5 %, то соответствующая ему квантиль известна. Известен также параметр X2 = а21/а22.Тогда искомое значение мощности можно найти по тем же таблицам, исходя из выражения:
(
1 -Р = P F(fi,/2)>
Fa( fl, /2 ) X2
Л
(4)
Встроенные в электронные таблицы библиотеки статистических функций, существенно облегчают задачу обращения к статистическим функциям. Поэтому уравнения (3) и (4) для работы в среде Excel можно записать одной формулой:
=РРАСП(БРАСПОБР(а/ /2)/X2; /1/),
(5)
которая сразу возвращает значение мощности критерия. В этой формуле в отличие от предыдущих вместо адресов ячеек целесообразнее использовать их имена.
С использованием нашей функции (5) весьма просто вычислить так называемые функции мощности, которые в виде графиков проводят в [1]. Пример организации вычислений таких функций показан на рис. 2.
При создании этой формы в В 6 записывают формулу РРАШ(РРАШОБР($Б$2;$В$3;$В$3)/А6;$В$3;$Б$3), которую затем копируют ниже по столбцу с помощью функции автозаполнения. Посредством формы на рис. 2 можно генерировать всевозможные графики для испытаний Фишера и не только при равенстве /1 = /2, как это делается в [1]. Более того, для отыскания параметров функции мощности практически отпадает необходимость в утомительной процедуре графической интерполяции.
Подставляя в ячейки на рис. 2 значения параметров для задачи с кернами, получаем 1 -р = 68 %, что не отличается от результата метода Монте-Карло на рис. 1. Согласие в результатах, с одной стороны, указывает на справедливость формулы (5), и, с другой стороны, говорит о надежности нашего алгоритма статистических испытаний.
Этот результат также косвенно указывает на качество использованного нами генератора случайных величин. Мы использовали самый простой генератор нормально распределенных величин без цифровой фильтрации «промахов» по правилу 3а, и округления данных для учета дискретности измерительных шкал.
Теперь решим задачу по оптимизации матрицы плана дисперсионного анализа в задаче с кернами. Для этого используем технику сценариев «Что если?» и вычисления по уравнению (5). Испытанные нами сценарии приведены в табл. 2. Видим, что увеличение числа проб по Сцена-
А В I С И
1 Вычисления функции мощности
2 Паоаметры функции а = 0.05
3 f1 = 5 f2 = 6
4 РеаМЬТаТ
5 X* 1 -р
6 1.0 0 050
7 1.5 0.112
8 2.0 0 183
9 2.5 0256
10 3.0 0 326
11 3.5 0 390
12 4.0 0 448
13 4.5 0 501
14 6.0 0 547
15 5.5 0 589
16 6.0 0 626
17 6.5 0 658
18 7.0 0 688
Рис. 2. Фрагмент формы для вычислений функции мощности
рию 2 на две единицы, против базового Сценария 1, позволяет удовлетворить требованиям эффективности теста: 1 - в > 80 %. Такого же результата достигаем в Сценарии 3. Здесь число проб меньше чем в Сценарии 2, но число опытов больше из-за увеличения числа параллелей до 3. Пожалуй, наиболее оптимален Сценарий 4, при котором число измерений не многим больше числа измерений по Сценарию 2, но трудозатраты на пробоотбор и пробоподготовку меньше (при очевидном выигрыше в надежности статистического вывода).
Таблица 2
Результаты испытаний «Что если?» (А,2 = 6,81, а = 0,05)
Сценарий т п тхп /1 /2 1 - в
1 2 6 12 5 6 0,68
2 2 8 16 7 8 0,80
3 3 6 18 5 12 0,80
4 3 7 21 6 14 0,85
5 3 8 24 7 16 0,89
Рис. 3 иллюстрирует возможности алгоритма генератора «Вычисления функции мощности» для визуализации результатов его работы. На рис. 3 показаны функции мощности для сценариев в табл. 2. Пунктирная линия соответствует X2 = 6,81.
Очевидно, что исследования мощности критерия методом Монте-Карло и методами тестовой статистики взаимно дополняют друг друга. Чем больше методов для решения одной и той аналитической задачи, тем надежнее не только обработка результатов, но и выводы из результатов химического анализа. Сочетание названных методов особенно полезно при преподавании химической метрологии и хемометрики, поскольку это отвечает требованиям современной дидактики, приоритетами которой являются компьютерные технологии, активные формы обучения, проблемное обучение, наглядные методы и повышение уровня мотивации к обучению.
Результаты работы могут быть востребованы при реализации концепции открытого образования и научных вычислительных веб-сервисов.
Литература
1. ГОСТ Р 50779.25-2005. Статистические методы. Статистическое представление данных. Мощность тестов для средних и дисперсий. - М., 2005. - 47 с.
2. Дерффель, К. Статистика в аналитической химии / К. Дерффель - М.: Мир, - 1984. - 268 с.
3. Дворкин, В.И. Метрология и обеспечение качества количественного химического анализа / В.И. Дворкин. - М.: Химия, 2001. - 263 с.
4. Аналитическая химия. Проблемы и подходы: пер. с англ. / под ред. Р. Кельнера,
Ж-М. Мерме, М. Отто, М. Видмера. - М.: Из-во АСТ, 2004. - Т. 2. - 728 с.
5. Головач, А.В. Критерии математической статистики в экономических исследованиях / А.В. Головач, А.М. Ерина, В.П. Трофимов. - М.: Статистика, 1973. - 135 с.
6. ГОСТ Р ИСО 5725-2-2002. Точность (правильность и прецизионность) методов и результатов измерений. - М., 2002. - Ч. 2. - 43 с.
7. Александровская, Л.Н. Использование ГОСТ требует бдительности / Л.Н. Александровская, О.М. Розенталь. - http://metrob.ru/HTML/Stati/gost-bdit.html (дата обращения 14.12.2012).
8. Голованов, В.И. Прогнозирование погрешностей фотометрии с использованием закона на-
0 5 10 15 а2 20
Рис. 3. Функции мощности для сценариев (номера на рисунке соответствуют номеру сценария в табл. 2)
копления ошибок и метода Монте-Карло / В.И. Голованов, Е.И. Данилина // Вестник ЮУрГУ. Серия «Химия». - 2010. - Вып. 3. - № 11 (187). - С. 20-26.
9. Голованов, В.И. Прогнозирование метрологических характеристик в титриметрии с использованием метода Монте-Карло / В.И. Голованов, Е.И. Данилина, Ю.С. Дворяшина // Вестник ЮУрГУ. Серия «Химия». - 2010. - Вып. 3. - № 11 (187). - С. 27-33.
10. Гмурман, В.Е. Теория вероятностей и математическая статистика / В.Е. Гмурман. - М.: Высшая школа, 1977. - 470 с.
11. Электронный учебник 81а!8ой. - http://www.statsoft.ru/home/textbook (дата обращения 14.12.2012).
Поступила в редакцию 12 июня 2012 г.
SIMULATION DISPERSION ANALYSIS IN SPREADSHEETS AND POWER TEST EVALUATION FOR COMPONENTS ERRORS DECOMPOSITION
The principles of stochastic simulation of the plan matrix for the dispersion analysis are proposed. The prospects for the use of stochastic modeling to solve problems of teaching and research in chemical metrology are discussed.
Key words: dispersion analysis, Monte Carlo method, statistical test, power test, computer technology in analytical chemistry, chemical metrology.
Golovanov Vladimir Ivanovich - Dr. Sc. (Chemistry), Professor, Head of Analytical Chemistry Subdepartment, South Ural State University. 76, Lenin avenue, Chelyabinsk, 454080.
Голованов Владимир Иванович - доктор химический наук, профессор, заведующий кафедрой аналитической химии, Южно-Уральский государственный университет. 454080, г. Челябинск, пр. им. Ленина, 76.
Е-mail: [email protected]
