Идентификация статистических моделей технологических процессов с заполнением пропусков в данных Текст научной статьи по специальности «Математика»

У

правление технологическими процессами

УДК 681.51.015

ИДЕНТИФИКАЦИЯ СТАТИСТИЧЕСКИХ МОДЕЛЕЙ ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ С ЗАПОЛНЕНИЕМ ОРООУСКОО В ДАННЫХ

Л. А. Кузнецов, А. М. Корнеев, М. Г. Журавлева

Липецкий государственный технический университет

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

ВВЕДЕНИЕ

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

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

точки зрения получения оптимальных оценок параметров итеративная стратегия EM-алгоритма (expectation-maximization algorithm), каждая итерация которого подразумевает выполнение двух шагов: шага E — вычисления математического ожидания и шага M — максимизации (метод Бака является частным случаем данного метода, реализующим его первые две итерации) [1]. Помимо оценивания параметров, с помощью EM-алгоритма для неполной многомерной нормальной выборки можно решать другие задачи статистического анализа, в том числе задачу построения множественной линейной регрессии. Рассмотрим подробнее методику EM-алгоритма в контексте идентификации моделей зависимости выходных свойств продукции, являющейся результатом сложного производственного процесса, от факторов технологии.

1. МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ ЕМ-АЛГОРИТМА

Пусть имеется к-мерная нормальная переменная х = (х1, х2, ..., хк), причем X = (Г1, Ґ2, ..., ґг) и (ї1, ї2, ...,

т

%_ г), где = У1р ..., ґп) — і-я технологическая величи-

т

на, і = 1, ..., г, Sj = (¿у, ..., ап) , — у-я выходная характеристика, у = 1, ..., к — г, полученные в результате пассивного эксперимента, п — число наблюдений. Она характеризуется параметрами 0 = (т, Е), где т = (ц1, ц2, ..., цк)т — вектор средних, Е = (аі7), і, І = 1, 2, ..., к — ковариационная матрица. Множество исходных данных удобно пред-

ставить в виде: х = (хнабл, х^), где хнабл - подмножество наблюдаемых значений факторов технологии и выходных свойств, хПр0П — подмножество пропущенных

значений. ТаКИм обрa3ом, хнабл = Кабл.Р ^абл^ ..., ^аб™^ где хнабл; — множество переменных (технологических величин и (или) выходных характеристик) с присутствующими значениями в наблюдении г, г = 1, 2, ..., и.

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

ной матрицы В = (Ь/ той же размерности, что исходная,

- - 1

х/ =

2 х1] , У = 1, 2, ..., к,

і і = 1

где п. — число наблюдаемых значений переменной у,

(х., если х. наблюдается, х>. = <1 >; >;

>; 1 0, если х. пропущено.

с элементами:

Ь = — 1/а

рр ' рр’

Ь. = Ь . = а. /а , р ф г,

1р рг гр рр’ ^ ’

Ь-- = а.. — а. ат./а„„, р ф г, р ф у.

г] у гр р/ рр’ ^ > г J

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

Обозначим исходную к-мерную переменную х, дополненную слева столбцом из единиц, через X. Приращенная ковариационная матрица строится следующим

_1 т

образом. Вычисляется произведение и XX, свертка которого по строке и столбцу у, j = 1, ..., к, позволяет получить регрессию всех остальных переменных на переменную с номером у. В общем случае свертка по нулевой строке и столбцу и первым г переменным позволяет получить многомерную регрессию переменных хг + 1, хг + 2, ..., хк на х1, х2, ..., хг. Результатом такой свертки будет блочная матрица:

swp[0, 1, ..., г]п 1ХтХ =

-Г G G Н

(2)

Если на т-й итерации рассчитаны параметры 0(т) = = (т(т), Е(т)), на шаге Е необходимо заменить пропущенные данные условными средними значениями по наблюдаемым в г-м опыте переменным хнабЛ

Е 2 х/хнабл, 0(т)| = 2 хі, У = 1, 2, ..., к,

Е 2

хііхіі\хнабл-

, 0(

(т) | _

(т) (т)

2 ХіГ х

іі

+ С( т) + С11і ,

у, I = 1, 2, ..., к,

(1)

где

х( т) = х1

х>;, если х>; наблюдается,

>> >>

Е(х>.хнабл, 0(т)), если х. пропущено,

а с.т — добавочные ковариации:

0, если х>. или х>7 наблюдаются,

С( т) = С11і

СОУ(х.Хнабл, 0(т)), если х>. и х>7 пропущены.

„( т)

Условные средние и добавочные ковариации из

формулы (1) вычисляются по текущим оценкам параметров посредством применения оператора свертки к приращенной ковариационной матрице [1] таким образом, что наблюдаемые в г-м опыте переменные рассматриваются как независимые, а те, в которых значение г пропущено — как отклики. Применение оператора свертки к некоторой симметричной матрице А = (а/ по строке и столбцу р приводит к получению симметрич-

где матрица G, размерностью (г + 1) х (к — г), у-й столбец которой содержит вычисленные методом наименьших квадратов свободный член и коэффициенты наклона

регрессии х.

на х1, х2, ..., хг, у = 1, 2, ..., к — г; Н — размерностью (к — г) х (к — г) остаточная ковариационная матрица хг + 1, х г +2, ..., хк при фиксированных х1, х2, ..., хг; элементы матрицы Г, размерностью (г + 1) х (г + 1), после умножения на соответствующие остаточные дисперсии и ковариации из матрицы Н, позволяют получить ковариационную матрицу коэффициентов регрессии (последнее справедливо для полных данных, в работе применялся метод оценивания данной матрицы для неполных данных [1, с. 176]).

На каждой итерации пропущенные значения заполняются новыми, более точными, по которым на шаге М заново вычисляются параметры:

(т + 1) —1

Ц/ = п

хіт), у = 1, 2, ..., к,

т + 1

._-1і

Е 2

І'і/'хі1| хнабл

(т + 1) (т + 1)

Ці

і= 1

у, I = 1, 2, ..., к.

Процесс сходится к некоторой оптимальной совокупности оценок параметров, по которым строится окончательное уравнение регрессии. В частности, для множества (?1, <2, ..., <г) факторов технологии и одного отклика * после завершения работы ЕМ-алгоритма дополнительно проводится соответствующая свертка. При этом матрица G из матрицы (2) будет иметь размерность (г + 1) х 1, т. е. являться вектором, содержащим свобод-

П

П

ный член и оценки параметров регрессии ^, *2, ..., ^ на я, а матрица Н — состоять из единственного значения, остаточной дисперсии полученной модели регрессии.

2. ПРИМЕНЕНИЕ ПРОЦЕДУРЫ ЕМ-АЛГОРИТМА К АНАЛИЗУ РЕАЛЬНЫХ ПРОИЗВОДСТВЕННЫХ ДАННЫХ С ПРОПУСКАМИ

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

каждой из трех моделей, отклику комплектной модели (по критериям х , Колмогорова — Смирнова [3]). На основании полученной информации был выбран оптимальный способ заполнения пропущенных значений и построена результирующая модель. Далее подробно изложены основные моменты и результаты проведенного анализа массива реальных производственных данных.

Исследовался массив с пропусками, содержащий набор наиболее существенных переменных, характеризующих ход и результаты технологического процесса производства 309 рулонов стали марки 08Ю. Из исходного массива было выделено 254 комплектных наблюдения. Посредством разработанной программной реализации ЕМ-алгоритма с использованием оператора свертки строились регрессионные модели зависимости выходной характеристики «глубина сферической лунки» от набора параметров сквозной технологии производства автолиста. По комплектным данным методом пошагового регрессионного анализа была построена следующая контрольная модель:

у = 11,543 — 4,264х1 — 0,873х2 — 0,002х3 + 0,004х4 —

— 0,020х5 — 0,034х6 — 0,035х7 — 0,034х8 + 0,002х9, (3)

где у — выходная характеристика, мм; х1 и х2 — массовые доли элементов химического состава, %, [С] и [А1] соответственно; х3 — скорость полосы при горячей прокатке, м/мин; х4 — температура конца горячей прокатки, °С; х5 ... х8 — обжатия по клетям стана холодной прокатки: соответственно, в 1-, 2-, 4- и 5-й клетях, %; х9 — скорость полосы на выходе первой клети стана холодной прокатки, м/мин (свободный член измеряется в мм; коэффициенты при переменных имеют размерности, обратные размерностям соответствующих переменных,

Таблица 1

Фрагмент статистического исследования построенных моделей

m, % Вид модели Norma RSS R F 2 X p (x2) D p (D)

0 (3) - 61,200 0,820 49,897 -

1 А В С 2,936 0,205 0,202 76,700 61,200 61,200 0,768 0,820 0,820 588 722 888 4, 9, 9, 344 7,462 0,529 0,571 0,382 0,999 0,999 0,079 0,028 0,028 0,410 1,000 1,000

10 A B C 14,650 1,099 0,825 1170,000 62,400 62,000 0,816 0,817 48,381 48,909 149,051 4,867 1,351 0,000 0,676 0,987 0,724 0,043 0,039 0,000 0,971 0,989

25 A B C 20,342 3,539 1,916 3110,000 90.000 67.000 0,720 0,801 26,096 43,454 227,664 19,107 9,983 0,000 0,008 0,190 0,961 0,079 0,059 0,000 0,410 0,768

40 A B C 17,296 11,881 1,682 000 000 ,0 ,0 ,8 0, 0, 4, 846 54 61 0,806 44,841 242,561 172,694 5,368 0,000 0,000 0,615 0,984 0,838 0,047 0,000 0,000 0,938

Примечания. A — модель, рассчитанная по данным, заполненным безусловными средними; B — по данным, заполненным условными средними; C — в результате работы EM-алгоритма; m — процент искусственно вводимых пропусков; Norma — норма разности векторов оценок параметров; RSS — остаточный квадрат; R — коэффициент множественной корреляции; F — значение статистики критерия Фишера; % — значение статистики критерия %; D — значение статистики критерия Колмогорова — Смирнова. Доверительная вероятность обозначена символом "p". Число групп для расчета значения критерия %2 — восемь, получено по формуле Штургеса [4]. При проверке адекватности моделей регрессии использовалось критическое значение критерия Фишера FKp. 9, 244, p = 0 99 = 2,481.

20 188 16 8 14-

Рн П го 12

сЗ

« 10-

а 8-

О

К 642-

А А л

ДД А

Д^ л^А , А Ад

^ ДДА Д<|> д

Д Д дд АА ф

д л д ♦

ДДА д ♦♦ А ♦ *

.....*$«•!••••••••••* *• *Ф* *

0 5 10 15 20 25 30 35 40 45

Число пропусков, %

Рис. 1. Зависимость значений норм разностей векторов оценок параметров от числа пропусков

умноженные на мм). Общий квадрат для модели (3) составил 186,739.

В табл. 1 приведен фрагмент статистического исследования для различных моделей и нескольких случаев искусственно полученных данных с пропусками. В соответствии с результатами, модели вида С наиболее устойчивы к пропускам в данных и близки к контрольной модели. Рассчитанные по ним значения норм разностей векторов оценок параметров являются минимальными из соответствующих значений для моделей вида В и С; остаточные квадраты, коэффициенты корреляции оптимальны и менее всего отличаются от соответствующих значений для контрольной модели. Модели вида А наиболее неустойчивы. Уже при 1 % пропусков остаточный квадрат, рассчитанный по модели вида А, примерно в 1,25 раза превышает соответствующее значение для контрольной модели, а для большего количества пропусков модели вида А неадекватны. Модели вида В близки к моделям вида С при небольших процентах пропущенных значений, в рассматриваемом случае это справедливо для т т 10 %. Перечисленные выводы подтверждаются общими результатами исследования, которое состояло в искусственной генерации пропусков, числом от 1 до 40 %, с шагом 1 (при 40 % пропусков и выше скорость и результаты работы алгоритма начинают существенно зависеть от структуры пропусков). На каждом шаге набор из трех моделей вида А, В и С генерировался семь раз [4], структура пропусков изменялась случайным образом по равномерному закону, значения рассчитываемых характеристик (норм разностей векторов, остаточных квадратов и др.) усреднялись.

На рис. 1—2 для указанных в табл. 1 видов моделей приведены зависимости норм разностей векторов оценок параметров и остаточных квадратов от числа пропусков.

Как показывают рисунки, модели вида А характеризуются недопустимо большими значениями остаточных квадратов, а их векторы оценок параметров значительно отличаются от соответствующего вектора контрольной модели. Для числа пропусков от 30 % и выше модели вида В существенно неадекватны.

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

Рис. 2. Зависимость остаточных квадратов от числа пропусков

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

2 до 6 % пропусков и близки к нулю, начиная с 7 %; для моделей вида В в диапазоне, примерно, от 18 до 29 % пропусков наблюдается большой разброс, а затем — нулевые значения рассматриваемых коэффициентов; для моделей вида С неустойчивость начинает проявляться при значениях около 35 % пропусков и выше.

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

у = 10,697 - 3,983хх - 0,783х2 - 0,003х3 + 0,004х4 -

- 0,017х5 - 0,030х6 - 0,040х7 - 0,027х8 + 0,002х9, (4)

где использованы введенные для модели (3) обозначения. Некоторые характеристики контрольной (3) и результирующей (4) моделей приведены в табл. 2.

Рис. 3. Зависимость значений коэффициентов детерминации от числа пропусков

Таблица 2

Характеристики контрольной и результирующей моделей

Контрольная модель: n = 254; RSS = 61,158; SS = 186,739; R = 0,820; F = 49,897; p (F)9, 244 = 1,000

Значения t-статистик для оценок параметров, ?кр244 = 099 = 2,59 6

h t2 t3 t4 t5 t6 t7 t8 t9

-2,813 -1,241 -5,631 2,221 2,329 -4,399 3,673 -2,885 2,562

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

xi x2 x3 x4 x5 x6 x7 x8 x9

-7,235 - -0,003 - - -0,049 0,054 0,056 -

-1,293 -0,002 -0,019 0,016 -0,011

Результирующая модель: n = 309; m = 2,1 %; Norma = 0,420; RSS = 71,375; SS = 231,389; R = 0,832; F = 63,445; p (F)9, 299 = 1,000.

Значения t-статистик для оценок параметров, tKp299 = 0 99 = 2,592

h t2 t3 t4 t5 t6 t7 t8 t9

-2,789 -1,336 -6,386 2,856 2,128 -3,980 4,264 -2,357 3,200

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

xi x2 x3 x4 x5 x6 x7 x8 x9

-6,783 - -0,003 0,001 - -0,044 -0,058 - 0,001

-1,184 -0,002 0,007 -0,015 -0,021 0,004

Примечание. В таблице использованы введенные ранее обозначения; п — число наблюдений; ££ — общий квадрат.

Модели (3) и (4) адекватны и пригодны для прогнозирования и управления. Заполнение пропущенных значений привело к смещению границ доверительных интервалов для значимых оценок параметров, а также переходу отдельных параметров в разряд значимо отличающихся от нуля и, наоборот, с заданной доверительной вероятностью. Такие колебания обусловлены близостью к нулю некоторых оценок параметров (соответствующие значения ¿-статистик мало отличаются от критических) и оказывают незначительное влияние на конечный результат. В рассматриваемой модели на выходную характеристику в большей степени влияют скорость горячей прокатки и обжатия во второй и четвертой клетях стана холодной прокатки, об этом свидетельствуют значения соответствующих ¿-статистик для обеих моделей. В целом, модель (4) улучшает характеристики исходной модели, не изменяя ее структуру. На это указывает малое значение нормы разности векторов оценок параметров контрольной и результирующей моделей. Даже при небольшом общем проценте пропусков (2,1 %) заметно увеличилось значение статистики критерия Фишера (от 49,897 до 63,445) и незначительно вырос коэффициент множественной корреляции (от 0,820 до 0,832).

ЗАКЛЮЧЕНИЕ

Анализ данных с искусственно введенными пропусками показал максимальную устойчивость характеристик моделей, построенных путем реализации ЕМ-алго-

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

ЛИТЕРАТУРА

1. Литтл Р. Дж. А., Рубин Д. Б. Статистический анализ данных с пропусками. — М.: Финансы и статистика, 1990. — 336 с.

2. Бард Й. Нелинейное оценивание параметров. — М.: Статистика, 1979. — 349 с.

3. Гайдышев И. Анализ и обработка данных: специальный справочник. — СПб.: Питер, 2001. — 752 с.

4. Львовский Е. Н. Статистические методы построения эмпирических формул. — М.: Высшая школа, 1988. — 239 с.

8 (4742) 46-53-54;

е-таг'1: kuznetsov@stu.lipetsk.ru

Статья представлена к публикации членом редколлегии

А. А. Дорофеюком. □