Анализ факторных экспериментов

 

1. Полные факторные эксперименты и кодированные переменные. При проведении экспериментальных исследований часто применяются факторные эксперименты.

Пусть – функция отклика и задан план, матрица которого есть (i = 1, 2,..., k; u = 1, 2,..., N), где – значение переменной (фактора) в u ‑м опыте. Каждое из различных значений, принимаемых переменной в эксперименте, называют уровнем этой переменной. Обозначим через s_i число различных уровней фактора X_i.

 

Определение 1. Эксперимент, в котором уровни каждого фактора комбинируются со всеми уровнями других факторов, называется полным факторным экспериментом.

Полный факторный эксперимент состоит из различных экспериментов, поэтому его называют экспериментом типа .

 

Определение 2. План называется симметричным, если все факторы имеют одинаковое число уровней, т. е. В этом случае полный факторный эксперимент принято называть экспериментом типа s^k, где k – число факторов.

 

Предположим, что число различных значений, которые может принимать переменная , в каждом опыте равно 2, т. е. s = 2. Тогда говорят, что переменная в каждом опыте варьируется на двух уровнях. Обозначим эти уровни через и . Если , то называют верхним уровнем, а – нижним уровнем фактора .

Обозначим , и введем новые переменные

 

(i = 1, 2,..., k).

 

Переменные x_i называют кодированными переменными. Легко проверить, что они могут принимать лишь два значения 1 (верхний уровень) и –1 (нижний уровень). Заменяя переменные кодированными переменными , можно представить функцию отклика в виде

 

. (1)

 

2. Полные факторные эксперименты типа 2² и 2³. Если число независимых переменных равно двум, то равенство (1) приводится к виду . Результаты наблюдений, отвечающих всевозможным комбинациям уровней переменных x ₁ и x ₂, сводятся в таблицу.

 

Номер опыта				Наблюдения
	– 1	– 1		y 1
		– 1	– 1	y 2
	– 1		– 1	y 3
				y 4

 

Пусть функция отклика имеет вид

 

 

и в каждом варианте испытаний проводится по одному наблюдению. Тогда имеем полный факторный эксперимент типа 2². Матрица плана и матрица планирования:

 

, .

 

Столбцы матрицы X попарно ортогональны, следовательно, планирование является ортогональным. Отсюда, используя результаты п. 7 § 1, получаем, что МНК‑оценки параметров некоррелированны и имеют вид

 

,

 

где X _j – j -й столбец матрицы Х, при этом .

 

Пусть теперь число переменных равно трем и функция отклика имеет вид

 

.

 

Рассмотрим все возможные комбинации уровней кодированных переменных , и .

Номер опыта								Наблю- дения
	– 1	– 1	– 1				– 1	y 1
		– 1	– 1	– 1	– 1			y 2
	– 1		– 1	– 1		– 1		y 3
			– 1		– 1	– 1	– 1	y 4
	– 1	– 1			– 1	– 1		y 5
		– 1		– 1		– 1	– 1	y 6
	– 1			– 1	– 1		– 1	y 7
								y 8

 

Здесь, как легко убедиться, планирование является ортогональным. Поэтому МНК-оценки параметров b _j некоррелированны и имеют вид

 

,

 

при этом .

 

Пример 1. Полный двухфакторный эксперимент проводится на двух уровнях.

 

Переменные	X 1	X 2
Нижний уровень
Верхний уровень

 

Получены результаты исследований: , , , .

Для нахождения МНК-оценок параметров ортогонального плана пользуются расчетной таблицей.

 

i
	– 1	– 1			– 66	– 66
		– 1	– 1			– 68	– 68
	– 1		– 1		– 48		– 48

С помощью расчетной таблицы получаем требуемые оценки:

 

; ;

; .

 

Значит, оценка функции отклика есть

 

. (2)

 

Возвращаясь к исходным переменным, имеем:

 

, , , , значит, , . Подставляя найденные выражения в (2), получаем

 

= = .

 

Пример 2. Полный трехфакторный эксперимент проводится на двух уровнях.

 

Переменные	X 1	X 2	X 3
Нижний уровень
Верхний уровень

 

Результаты исследований: ; ; ; ; ; ; ; .

Составим расчетную таблицу для нахождения МНК‑оценок параметров ортогонального плана.

i
x 1	–1		–1		–1		–1
x 2	–1	–1			–1	–1
x 3	–1	–1	–1	–1
x 1 x 2		–1	–1			–1	–1
x 1 x 3		–1		–1	–1		–1
x 2 x 3			–1	–1	–1	–1
x 1 x 2 x 3	–1			–1		–1	–1
yi	5,6	7,7	8,1	9,6	8,6	5,1	6,4	6,9
x 1 yi	–5,6	7,7	–8,1	9,6	–8,6	5,1	–6,4	6,9
x 2 yi	–5,6	–7,7	8,1	9,6	–8,6	–5,1	6,4	6,9
x 3 yi	–5,6	–7,7	–8,1	–9,6	8,6	5,1	6,4	6,9
x 1 x 2 yi	5,6	–7,7	–8,1	9,6	8,6	–5,1	–6,4	6,9
x 1 x 3 yi	5,6	–7,7	8,1	–9,6	–8,6	5,1	–6,4	6,9
x 2 x 3 yi	5,6	7,7	–8,1	–9,6	–8,6	–5,1	6,4	6,9
x 1 x 2 x 3 yi	–5,6	7,7	8,1	–9,6	8,6	–5,1	–6,4	6,9

 

Используя расчетную таблицу, так же, как и в предыдущем примере, получаем требуемые оценки параметров: ; ; ; ; ; ; ; . Следовательно,

 

. (3)

 

Возвращаемся к исходным переменным.

; ; ; ; ; , значит, , , . Подставляя эти выражения в (3), получаем окончательно

 

3. Факторные эксперименты с повторными наблюдениями. Пусть – функция отклика, ; m, m,..., m – спектр плана (число наблюдений в каждой точке x _l одно и то же и равно m), mn = N. Пусть x _l = (x _{1 l}, x _{2 l},..., x_kl) и – повторные наблюдения в точке x _l .

Матрицу плана можно представить в виде матрицы из m блоков следующим образом:

 

, где .

 

Таким образом, – матрица размеров с различными строками, – матрица размеров .

 

Определение. План называется полным факторным планом типа 2 ^k с повторными наблюдениями кратности m, если матрица является матрицей полного факторного плана типа 2 ^k.

 

Пример. Пусть – функция отклика и задан план: , , , , причем в каждой точке число наблюдений m = 2. Тогда матрица плана записывается в виде:

 

, где .

 

Так как – матрица плана полного факторного эксперимента типа 2², то данный план является полным факторным планом с повторными наблюдениями кратности 2.

Матрица планирования имеет вид

 

 

и является матрицей ортогонального планирования.

Если – матрица планирования, соответствующая функции отклика и матрице плана , то из ортогональности планирования для МНК-оценки вектора β имеем равенство , в котором , , а значит, (j = 0, 1,..., p). В этом случае оценки некоррелированны и имеют дисперсию .

 

4. Насыщенное и ненасыщенное планирование.

Определение 1. Пусть матрица планирования X имеет ранг, равный числу p неизвестных параметров в функции отклика. Факторный план называется насыщенным, если p = N, и ненасыщенным, если p < N.

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

Пример 1. Пусть и планируется полный факторный эксперимент типа 2². Тогда планирование будет насыщенным, т. к. ранг матрицы планирования X равен 4 и при этом число наблюдений N и число неизвестных параметров p также равны 4.

Пример 2. Если и имеется полный факторный эксперимент типа 2², то планирование будет ненасыщенным, т. к. ранг матрицы планирования

 

равен 3, а .

 

Определение 2. Пусть n – число точек спектра факторного плана с повторными наблюдениями кратности m (N = mn) и r – ранг матрицы планирования X. Тогда план называется насыщенным, если , и ненасыщенным, если .

 

5. Проверка гипотезы адекватности. Рассмотрим факторный эксперимент с кратными повторными наблюдениями y_ls (l = 1, 2,..., n; s = 1, 2,..., m). Проверим гипотезу H ₀ адекватности модели

 

, (4)

 

где – известные функции, задаваемые равенствами вида (1 m i ₁ < i ₂ <... < i_q m k), а – неизвестные параметры. Функции отклика (4) и матрице (i = 1, 2,..., k; u = 1, 2,..., N; N = mn) полного факторного плана соответствует матрица планирования (j = 1, 2,..., p ₀; u = 1, 2,..., N). Столбцы матрицы X ⁰ должны удовлетворять условиям:

 

(j = 2, 3, …, p); (5)

(j = 1, 2, …, p); (6)

 

. (7)

 

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

Гипотеза состоит в том, что M y = X ⁰ β ⁰, где . Эта гипотеза проверяется при альтернативной гипотезе : M y ¹ X ⁰ β ⁰.

Для проверки гипотезы нужно вычислить отношение . Величина представляет собой несмещенную оценку для . Полагая в формуле (6) § 2 m ₁ = m ₂ =... = m_i = m, получаем

 

.

 

Таким образом,

 

. (8)

 

Оценка (5) § 2 дисперсии , связанная с неадекватностью модели, есть , где r – ранг матрицы X ⁰.

Пусть r = , тогда, учитывая (5) – (7), получаем

 

, (9)

 

где , j = 1, 2,..., p ₀.

Далее вычисляем и по таблице распределения Фишера по уровню значимости a и степеням свободы n – r и N – n находим .

Если , то гипотеза H ₀ принимается.

Если , то гипотеза H ₀ отклоняется.

 

Пример. Предположим, что зависимость прочности бетона R от двух факторов – расхода цемента X ₁ и расхода воды X ₂ – имеет вид: . Полный двухфакторный эксперимент проводится на двух уровнях:

 

Уровень	Расход цемента X 1	Расход воды X 2
Нижний
Верхний

 

Результаты исследований представлены таблицей.

 

i
	28,6	31,1	29,5	29,7
	44,3	47,8	46,2	46,1
	22,9	24,6	21,9	23,1
	38,7	38,7	35,3	36,7

МНК-оценки параметров ортогонального планирования находятся так же, как в примере 1 п. 3: ; ; , отсюда оценка функции отклика – , или, переходя к исходным переменным, .

Проверим, является ли полученная модель адекватной. Имеем N = 12, n = 4, m = 3, r = 3. По формулам (8) и (9)

 

,

 

.

Следовательно, . Далее по таблице критических точек распределения Фишера по уровню значимости a = 0,05 и степеням свободы n – r = 1 и N – n = 8 находим . Так как , то гипотеза, утверждающая, что модель адекватна, принимается.

 

 

§ 4. Оптимизация планов. Поиск экстремума функции отклика

 

1. Построение линейных оптимальных планов. Пусть функция отклика является полиномом первой степени от k переменных:

. (1)

 

Определение 1. План называется линейным, если ему соответствует функция отклика (1) и он позволяет получить несмещенные МНК-оценки параметров .

 

Если матрица линейного плана есть

, (2)

 

то матрица планирования для этого плана имеет вид

 

.

 

Определение 2. Линейный план называется линейным ортогональным планом, если столбцы матрицы X попарно ортогональны, т. е.

 

. (3)

 

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

Пусть матрица планирования X имеет ранг k + 1 и вектор наблюдений удовлетворяет условиям ; , тогда для ковариационной матрицы МНК-оценок вектора β имеем равенство

 

.

 

Будем предполагать, что

 

, , (4)

где величины c_i заданы; положим также .

Требуется выбрать матрицу планирования X или матрицу линейного плана , так чтобы дисперсии МНК-оценок параметров в классе линейных планов с ограничениями (4) были минимальными.

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

 

Теорема Бокса. Пусть функция отклика имеет вид (1), столбцы матрицы линейного плана удовлетворяют условиям (4) и ранг матрицы X равен k + 1. Тогда для МНК-оценок параметров выполняется неравенство

 

(), (5)

 

причем минимум дисперсий в классе линейных планов с ограничениями (4) достигается тогда и только тогда, когда столбцы матрицы X попарно ортогональны.

В частности, поскольку первый столбец матрицы X состоит из единиц (x _{0 u} = 1), из условия ортогональности (3) при j = 0 следует . Таким образом, получаем

Следствие. Необходимым условием минимума дисперсий МНК-оценок является выполнение при всех i = 1, 2,..., k равенства (условие симметрии плана).

 

2. Градиентный метод. Будем предполагать, что функция отклика непрерывна и имеет непрерывные частные производные 1-го порядка в ограниченной замкнутой области G Í Ñ ^k.

 

Определение 1. Функция f называется унимодальной в области G, если она в этой области имеет единственный экстремум.

 

Пусть изменение функции многих переменных вдоль некоторой траектории, проведенной из точки в точку и определенной в области G, задается уравнением , где l – параметр, пробегающий числовой отрезок [0; 1].

Траекторию назовем строго возрастающей (строго убывающей), если функция h(l) строго возрастает (строго убывает) на отрезке [0; 1].

 

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

 

Экстремум функции отыскивается методом спуска (или подъема) по исследуемой поверхности. Этот метод основан на построении последовательности точек , лежащих в области G и таких, что (метод спуска) или (метод подъема).

Наиболее часто применяемая в практике исследований разновидность метода спуска и подъема – так называемый градиентный метод, при котором последовательность точек определяется с помощью равенства

 

, (6)

 

где – градиент функции f в точке ; a – некоторое положительное число. В этом случае функция f предполагается строго унимодальной.

Различные варианты градиентного метода отличаются друг от друга способом выбора a. Одним из вариантов градиентного метода является метод наискорейшего спуска (или подъема). При оптимизации технологических процессов в теории планированияэксперимента используется его статистический аналог – метод Бокса и Уилсона. В этом методе используется не сам градиент, а его оценка.

 

3. Оценивание градиента. Пусть функция отклика

 

(7)

 

определена в области G Ì Ñ ^k.

Выберем произвольно точку . Используя эту точку, построим полный факторный эксперимент. Выберем для каждого i = 1, 2,..., k нижний уровень и верхний уровень так, чтобы было . Положим , введем кодированные переменные и выразим функцию отклика (7) через кодированные переменные:

 

. (8)

 

Под задачей оценивания градиента будем понимать определение оценки градиента функции отклика (8) в точке x ⁰ = (0, 0,..., 0). Предположим, что в окрестности точки x ⁰ функция (8) допускает разложение по формуле Маклорена:

 

(9)

где . Введем обозначения: , , , . Тогда равенство (9) перепишется в виде

 

.

 

Так как , задача оценивания градиента сводится к нахождению МНК-оценок параметров .

Пусть матрица полного факторного плана с центром в точке задана равенством (2).

Для простоты будем считать, что в замкнутой области функция отклика достаточно точно аппроксимируется линейной функцией, т. е.

 

. (10)

 

Тогда матрице плана и функции отклика (10) соответствует матрица планирования (j = 0, 1, 2,..., k; u = 1, 2,..., N; x _{0 u} = 1) и для МНК-оценки параметра имеем равенство

 

 

(j = 0, 1, 2,..., k),

 

где – наблюдения в точках плана.

Поскольку являются оценками компонент градиента , т. е. частных производных , то МНК-оценка градиента функции отклика в точке определяется равенством

 

.

 

Пример. Пусть h = f (x ₁, x ₂) – функция отклика (x ₁, x ₂ –кодированные переменные), матрица плана и результаты наблюдений . Используя аппроксимацию вида h» b₀ + b₁ x ₁ + b₂ x ₂ + b₃ x ₁ x ₂ (), найдем оценку градиента в центре плана (0, 0). Так как и , то планирование ортогонально, = (– y ₁ + y ₂ – y ₃ + y ₄): 4 = –3, = (– y ₁ – y ₂ + y ₃ + y ₄): 4 = –10 и, следовательно, получаем .

 

4. Метод Бокса и Уилсона. Будем предполагать, что функция отклика (7) в области G строго унимодальна. Поиск ее максимума может быть осуществлен методом, предложенным Боксом и Уилсоном.

Пусть – начальная точка поиска максимума. При переходе к кодированным переменным точке соответствует точка , а функция отклика (7) принимает вид (8).

Пусть – оценка градиента в точке , полученная с использованием факторного эксперимента, матрица плана которого задана равенством (2). Для поиска максимума сделаем некоторый шаг из точки в направлении оценки градиента . Положим

 

, (11)

 

где параметр , – МНК-оценка градиента, , (i = 1, 2,..., k).

Точке в системе координат кодированных переменных соответствует точка в системе координат исходных переменных, при этом имеем .

Выполним в точке наблюдения и найдем в ней оценку функции отклика

 

. (12)

 

Так как измерение функции отклика происходит с ошибкой, то в точке можно найти только ее оценку, а не точное значение. Предположим, что оценка значимо больше оценки функции отклика в точке , где согласно (10) .

В этом случае делаем второй шаг в направлении оценки градиента и т. д. В общем виде, учитывая (11), для l -го шага имеем:

 

, (i = 1, 2,..., k), (13)

 

где