Анализ эксперимента по плану 2 в 4 степени

11.3. Анализ эксперимента по плану 2⁴

При небольшом числе факторов общее число опытов в эксперименте по плану 2^р довольно большое. Например, по плану 2⁵ в эксперименте будут 32 опыта, по плану 2⁶ эксперимент состоит из 64 опытов, а по плану 2¹⁰ эксперимент включает уже 1024 опыта. Выделяемые для проведения эксперимента средства обычно ограничены, так что часто выполнить все опыты можно только один раз, как в примере эксперимента с пневматическим прибором из предыдущего раздела. Выполненный так эксперимент называют экспериментом без повторения опытов. По результатам эксперимента стараются оценить не только воздействия факторов и их взаимодействий на отклик, но и ошибки результатов оценки. Однако без повторения опытов дисперсию ошибок эксперимента можно оценить только с использованием постулируемой модели. В разделе 7.3 показано, что несмещённую оценку дисперсии s² можно сделать, используя вектор несмещённой оценки параметров этой модели.

Кроме этого, по результатам такого эксперимента нельзя проверить адекватность постулируемой модели. Для проверки адекватности модели в разделе 6.4 некоторые опыты эксперимента повторялись. Поэтому в данном разделе рассмотрим эксперимент, в котором все опыты повторялись несколько раз. Этот эксперимент из области чистовой токарной обработки заготовок на станке с числовым программным управлением (ЧПУ).

Шероховатость поверхности является наиболее важной характеристикой качества многих изделий получаемых на токарных станках в машиностроении. При планировании процесса токарной обработки, вместе с заданием допуска на размер, требование к шероховатости является одним из важнейших ограничений. Проанализируем результаты исследований, проведённых с целью выяснения влияния на шероховатость поверхности при токарной обработке таких факторов, как твёрдость материала обрабатываемой заготовки, скорость продольной подачи, угол при вершине резца в плане, глубина резания и скорость вращения шпинделя станка [Feng (2001); Feng, Wang (2002)]. При этом радиус закругления вершины резца оставался постоянным и равным 0,819 мм.

Для выяснения степени влияния рассматриваемых факторов на переменные отклика, являющиеся средними арифметическими измеряемых шероховатостей поверхностей обработанных заготовок, был поставлен отсеивающий эксперимент. Эксперимент по дробному факторному плану 2^5–1 разрешающей способности V имел шестнадцать опытов в виде определённых операций токарной обработки [Feng (2001)]. Все опыты этого эксперимента выполнялись в случайной последовательности с повторением каждого три раза, и получаемые значения шероховатости записывались в микронах.

Эксперимент и его результаты

Регрессионным анализом данных показано, что твердость материала обрабатываемой заготовки, скорость продольной подачи, угол при вершине резца в плане и скорость резания оказывают существенное влияние на шероховатость поверхности [Feng (2001)]. Четыре взаимодействия между двумя разными факторами также значительно влияют на шероховатость. Однако фактор "глубина резания", принимавший в эксперименте значения 0.51mm и 1.02mm, не имел существенного влияния. Поэтому, исключая этот фактор, статистическую модель зависимости шероховатости от указанных выше факторов можно построить на основе результатов эксперимента выполняемого по полному факторному плану 2⁴ [Боков (2016)].

В опытах эксперимента требуемые по плану токарные работы выполнялись на цилиндрических заготовках длиной 76,2 мм и диаметром 31,75 мм. Твёрдость 86 HRB материала заготовок на верхнем уровне этого фактора обеспечивалась выбором в качестве материала никель-хром-молибденовой цементированной стали (AISI 8620), а твёрдость 52HRB материала заготовок на нижнем уровне этого фактора обеспечивалась выбором алюминиевого сплава (типа AD-33). В резцах использовались вставные твердосплавные (режущие) пластины с покрытием углеродно-нитридным титаном Ti (C, N), Al2O3 и TiN (карболой класса TP200). Угол при вершине резца в плане на верхнем уровне был 80º, а на нижнем уровне составлял 35º.

Опыты проводились на токарном станке с ЧПУ типа YAM CK-1. Скорость вращения шпинделя станка на верхнем уровне устанавливалась равной 1500 об/мин, а на нижнем уровне была равной 1000 об/мин. Скорость продольной подачи на верхнем уровне устанавливалась равной 0,127 мм/об, а на нижнем уровне была 0,051 мм/об.

Рекомендуемые материалы

Измерения шероховатости поверхности обработанных заготовок проводились с использованием профилометра Митутойо типа SJ-301. Три измерения выполнялись вдоль оси каждой обработанной заготовки по её обработанной поверхности, и каждое измерение отстояло приблизительно на 120° по спирали в последовательности один за другим. Усреднённое значение этих трёх измерений представляло наблюдаемое значение арифметического среднего шероховатости поверхности.

В рассматриваемом эксперименте по плану 2⁴ все опыты подлинно повторялись три раза. Под подлинными повторениями опыта подразумевается то, что вариация между выполняемыми при одинаковых условиях опытами эксперимента является отражением общей изменчивости от опыта к опыту. Рандомизация последовательности выполнения всех 48 опытов эксперимента гарантирует, что повторения являются подлинными. Но выполнение этого не всегда бывает легким. Для данного эксперимента опыт включал: (1) наладку станка, (2) загрузку программы обработки заготовки в соответствии с планом эксперимента, (3) установку на станке заготовки из требуемого по плану материала (4) выполнение токарной обработки при всех факторах, установленных на соответствующих уровнях для этого опыта и (5) измерение шероховатости поверхности обработанной заготовки. Подлинное повторение опытов означает выполнение всех этих операций для всех 48 опытов. Повторение только измерений шероховатости поверхности обработанной заготовки в каждом опыте дало бы только оценку дисперсии измерений. Обычно это только малая часть дисперсии эксперимента от опыта к опыту. Для выполнения оценки дисперсии ошибок эксперимента, результат которой используется для нахождения  стандартной ошибки результата оценки отдельного воздействия, каждый опыт эксперимента должен подлинно повторяться. Результаты выполненных таким образом в случайной последовательности 48 опытов представлены в таблице 11.3.1 элементами векторов у₁, у₂, у₃, а вектор  состоит из усреднённых значений результатов повторений всех опытов.

Таблица 11.3.1. Компоновка факторов и данные отклика эксперимента по плану 2⁴

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

Таблица 11.3.2. Матрица полной модели для эксперимента по плану 2⁴

Анализ эксперимента по усреднённым результатам опытов

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

=β₀+β₁x₁+β₂x₂+β₃x₃+β₄x₄+β₁₂x₁₂+β₁₃x₁₃+β₁₄x₁₄+β₂₃x₂₃+β₂₄x₂₄+β₃₄x₃₄

+β₁₂₃x₁₂₃+β₁₂₄x₁₂₄+β₁₃₄x₁₃₄+β₂₃₄x₂₃₄+β₁₂₃₄x₁₂₃₄+e,

где коэффициенты β с соответствующими подстрочными индексами - параметры модели, а нормированные факторы x₁ –твёрдости материала заготовок, x₂ –скорости продольной подачи, x₃ - угла при вершине резца, x₄ - скорости вращения шпинделя и e - ошибка опыта. При выбранных выше верхних и нижних уровнях включённых в модель факторов, их нормированные значения находятся по формуле (2.6.4), как и в (11.2.1). В матричном виде модель записывается так

=X_оβ_о+e,                                                    (11.3.1)

где матрица Х_о дана в таблице 11.3.2, вектор  приведён в таблице 11.3.1, β_о - вектор параметров модели и e - вектор ошибок опытов эксперимента.

В предыдущем разделе показано, что при двухуровневом плане эксперимента его линейная модель имеет параметры, являющиеся фактически воздействиями факторов и их взаимодействий на переменные отклика. Поэтому, считая верными допущения раздела 7.1, где Е(у)=Xq заменяется на Е(у)=X_оβ_о, методом наименьших квадратов можно найти вектор оценки параметров модели (11.3.1) по формуле

= (X_о^TX_о)^–1X_о^T.                                     (11.3.2)

А если воспользоваться полученными в результате выполнения эксперимента  числовыми значениями переменных отклика [Feng (2001)], то вектор оценки параметров получается

= [1,354  –0,201  0,084  –0,089  –0,229  0,038  –0,237  0,054  0,020  0,048  0,151

  0,031  0,050  0,057  –0,104  –0,022].         (11.3.3)

На основе результатов 16 опытов эксперимента по плану 2⁴ невозможно оценить дисперсию ошибок эксперимента, но обоснованным предположением может быть то, что некоторые результаты оценки воздействий от взаимодействий (т.е., параметров модели) получаются, в основном, из-за шума и поэтому их можно считать множеством результатов для сравнения с остальными результатами оценки параметров.

При анализе воздействий на основе данных факторных экспериментов без повторения опытов появляются две трудности: (1) иногда получаются большие результаты оценки воздействий от взаимодействий высокого порядка и (2) необходимо предусматривать отбор воздействий. Графическое изображение результатов оценки воздействий в виде рассмотренного в разделе 1.6 графика кумулятивных вероятностей их распределения может обеспечить эффективный путь преодоления обеих трудностей. Такие графики строятся обычно с использованием компьютерных программ. Например, в программе Mathcad эти графики строятся следующим образом.

На графике кумулятивных вероятностей распределения каждое значение данных строится в зависимости от такого его значения, которое имело бы нормальное распределение. Рассмотрим те элементы вектора (11.3.3), которые представляют результаты оценки воздействий только контролируемых факторов и их взаимодействий. Эти элементы обозначим вектором . Ожидаемые нормальные значения, называемые нормальными подсчётами (N_Score), могут быть оценены вычислением сначала ранговых подсчётов сортированных данных. В Mathcad это делается так

j:=0..last()

Effect:=sort()

r_j:=j+1

И, формируя пропорциональную оценку для каждого rank_j получаемую по следующей формуле, называемой иногда формулой Van der Waerden,

p_j:=,

нормальные подсчёты являются соответствующими p_j процентильными пунктами из стандартного нормального распределения

N_Score_j:=qnorm(p_j,0,1)

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

Положим, что данные анализируемого эксперимента по плану 2⁴ получились просто в результате случайной вариации около среднего и изменения уровней факторов не имели реального влияния на отклик. Тогда эти 15 результатов оценки воздействий, представляющих 15 контрастов между парами усреднённых переменных отклика, имели бы приблизительно нормальные распределения с центром в нуле и располагались бы на графике их кумулятивных вероятностей приблизительно по прямой линии. Однако на Рис.11.3.1 видно, что результаты оценки воздействий не все располагаются вдоль прямой. Довольно близко к прямой линии располагаются 8 результатов оценки, включая результат оценки воздействия от взаимодействия четырёх факторов. Это «линия ошибок», так как образующие её результаты оценки таковы, что могли бы быть результатом только шума. Однако остальные результаты оценки резко отклоняются от этой линии и не объясняются случайным происхождением. Поэтому может быть резонным допущение, что восемь результатов оценки по линии ошибок в основном из-за шума и поэтому составляют множество для сравнения с остальными результатами оценки.

Рис.11.3.1. График кумулятивных вероятностей распределения для результатов оценки воздействий

В это множество для сравнения входят результаты оценки воздействий от взаимодействий двух =0,038, =0,054, =0,020, =0,048, трёх =0,031, =0,050, =0,057 и четырёх =–0,022, факторов. Так как математические ожидания этих результатов оценки равны нулю, то их совместная оценка дисперсии выполняется по формуле

(СОО)²=(²+²+²+²+²+²+²+²)/8=1,781х10^–3

и является квадратом стандартной ошибки (СОО) результатов оценки воздействий остальных факторов и их взаимодействий. Поэтому СОО==0,042.

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

=β₀+β₁x₁+β₂x₂+β₃x₃+β₄x₄+β₁₃x₁₃+β₃₄x₃₄+β₂₃₄x₂₃₄+e

или в матричном виде

=Xβ+e,                                                      (11.3.4)

где X - матрица уровней факторов и отобранных их взаимодействий, а β - вектор параметров этой модели. Вектор оценки этих параметров методом наименьших квадратов находится по формуле = (X^TX)^–1X^T и получается

=[1,354  –0,201  0,084  –0,089  –0,229  –0,237  0,151  –0,104].           (11.3.5)

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

S_E=[I–X(X^TX)^–1X^T]=0,228

и из выражения (7.3.8) несмещённая оценка дисперсии получается s²=S_E/8=0,029. Обратная матрицы X^TX

(X^TX)^–1=.

Отсюда оценка дисперсии для результатов оценки параметров делается по формуле 0,063s²=1,796х10^–3, а стандартная ошибка результатов оценки параметров, равная корню квадратному из результата оценки их дисперсии, получается равной 0,042. Это значение совпадает с полученным методом множества для сравнения.

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

Таблица 11.3.3. Результаты оценки параметров со стандартными ошибками для данных эксперимента по плану 2⁴

Из таблицы 11.3.3 также видно, что все факторы эксперимента взаимодействуют друг с другом. Поэтому они должны рассматриваться вместе. Можно также графически представить взаимодействия факторов, как было сделано на Рис.11.2.2. При этом взаимодействие трёх факторов х₂, х₃ и х₄ можно рассматривать как взаимодействие фактора х₂ с взаимодействием двух факторов х₃ и х₄.

Построение и анализ графика кумулятивных вероятностей для распределения оценок воздействий является важным в отборе воздействий, не относящихся к шуму. Другой графический метод отбора воздействий предложен в [Lenth, (1989)]. Его описание можно найти также в [Box с соавт. (2005) стр. 203].

Для правильного использования статистических методов информация о результатах оценки воздействий факторов и их взаимодействий с возможной ошибкой должна анализироваться совместно с имеющимися знаниями предметной области. Графические методы дают ценные средства, позволяющие информации в данных и знаниям предметной области соответствующим образом взаимодействовать. Например, для используемых уровней факторов рассматриваемого эксперимента специалист по чистовой токарной обработке мог и не подозревать о таких больших воздействиях от взаимодействий двух факторов х₁ и х₃, х₃ и х₄, а также трёх факторов х₂, х₃ и х₄.

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

В результате анализа эксперимента с использованием его полной линейной модели (11.3.1) установлено, что из всех взаимодействий факторов только воздействия от взаимодействий двух факторов х₁ и х₃, х₃ и х₄, а также трёх факторов х₂, х₃ и х₄, имеют выделяющиеся из шума результаты оценки. Поэтому дальнейший анализ эксперимента с учётом повторных опытов имеет смысл проводить с использованием функции линейной модели (11.3.4). Все опыты эксперимента повторялись трижды, следовательно, с учётом повторных опытов эта модель в матричном виде записывается так

у=СXβ+e,                                                     (11.3.6)

где у= - вектор всех измеренных значений шероховатости, С= - матрица преобразования размеров тхп и ранга п=16, при т=3п=48, X - матрица модели, состоящая из уровней нормированных факторов и их отобранных взаимодействий с первым столбцом из единиц, как показано в таблице 11.3.2, β - вектор параметров модели и e - вектор ошибок эксперимента, которые, полагается, распределены нормально в виде N(0, Is²). В таблице 11.3.4 уровни факторов и опыты размещены в стандартном порядке. Матрица С состоит из трёх единичных матриц, так как все опыты повторялись трижды. Но если опыты повторяются не одинаковое число раз, то матрица С формируется в соответствии с числом повторений каждого опыта, например, как показано в разделе 6.4.

Таблица 11.3.4. Матрица модели для эксперимента по плану 2⁴

Модель (11.3.6) содержит линейную функцию вектора β неизвестных параметров, а элементами её матрицы X являются уровни нормированных факторов и их отобранных взаимодействий. В этом случае, как показано в предыдущем разделе, параметры модели выражают воздействия соответствующих им факторов на переменные отклика. По теореме Гаусса-Маркова оценка параметров получается наилучшей, если она выполняется методом наименьших квадратов. Однако для такой оценки вектора β необходимо соблюдение допущений раздела 7.1, где Е(у)=Xq нужно заменить на Е(у)=СXβ. Тогда формула оценки вектора β имеет вид

= (X^TС^TСX)^–1X^TС^TR

=(X^TX)^–1X^TС^TR/3.                                     (11.3.7)

В числовом виде, если воспользоваться полученными значениями переменных отклика [Feng (2001)], то вектор оценки параметров получается в виде

= [1,354  –0,201  0,084  –0,089  –0,229  –0,237  0.151  –0,104],

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

=1,354–0,201x₁+0,084x₂–0,089x₃–0,229x₄–0,237x₁₃+0,151x₃₄–0,104x₂₃₄.          (11.3.8)

Диагностическая проверка модели

В разделе 8.4 отмечено, что остатки играют важную роль в диагностической проверке модели. Вектор остатков для модели (11.3.6) находится по формуле =у–, где вектор оценки ожидаемых значений переменных отклика =СX. Построение графика кумулятивных вероятностей распределения остатков даёт диагностическую проверку любой модели эксперимента. Причём на этом же графике можно построить прямую линейной функции, аппроксимирующей зависимость нормальных подсчётов от элементов вектора остатков с использованием метода наименьших квадратов. Результаты этих построений показаны на Рис.11.3.2.

Рис.11.3.2. График кумулятивных вероятностей распределения остатков и аппроксимирующая прямая линия

Все значения остатков на Рис.11.3.2 располагаются вблизи или прямо по прямой линии, подтверждая предположение, что все не включённые в модель (11.3.6) воздействия от взаимодействий факторов легко объясняются случайным шумом. Эта проверка является полезной, так как подтверждает, что только три из всех воздействий от взаимодействий факторов действительно влияют на переменные отклика.

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

Рис.11.3.3. График зависимости остатков от результатов оценки ожидаемых значений переменных отклика.

Однако если величина дисперсии зависит от значений переменных отклика, то этот график покажет воронкообразный рисунок. Это наводит на мысль о необходимости преобразования переменных отклика [Myers с соавт. (2016) стр. 37]. Рекомендуется также построение графиков остатков в зависимости от порядка их получения и от уровней факторов. И если на этих графиках значения остатков выстраиваются не в случайном порядке, то это указывает на отклонение распределения ошибок от нормального закона и преобразование переменных отклика может способствовать достижению нормальности.

Проверка значимости регрессии

В разделе 9.1 показано, что проверка значимости регрессии заключается в проверке гипотезы H₀: β₁=0, где в данном случае β₁^Т= [β₁, β₂, β₃, β₄, β₁₃, β₃₄, β₂₃₄]. Для этой проверки используется статистика F_Rc, расчёт которой выполняется по приведённым в таблице 9.1.1 формулам. Однако в этих формулах не учитывается повторение опытов. Поэтому для расчёта статистики F_Rc необходимы формулы с учётом повторений опытов.

Для получения таких формул разделим матрицу Х в виде [1_п, X₁], а вектор β соответственно в виде . Тогда модель (11.3.6) можно записать так

у= [С1_п, СX₁]+e.                                  (11.3.9)

Аналогично выражению (7.5.28), нормальные уравнения для этой модели записываются в виде

[С1_п, СX₁]^Т[С1_п, СX₁]=[С1_п, СX₁]^Тy.                       (11.3.10)

Произведение [С1_п, СX₁]^Т[С1_п, СX₁] в левой части можно преобразовать так

[С1_п, СX₁]^Т[С1_п, СX₁] =[С1_п, СX₁]=,

так как С1_п=1_т. Матрица С^TС диагональная с элементами по диагонали равными числу повторений опытов. Все опыты повторялись трижды, поэтому диагональ состоит из чисел 3. В силу (П.2.21) и (П.2.22), можно поменять местами матрицы С^TС и X₁, а также С^TС и 1, и записать последнее выражение в виде

==,

так как 1_т^Т1_т=m. Взаимодействия факторов сами можно рассматривать в качестве новых нормированных факторов. Поэтому, как показано после (7.5.29), имеем 1^ТX₁=0^Т.

Правую часть выражения (11.3.10) можно записать в виде

[1_т, СX₁]^Тy=y==.

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

={[1_т, СX₁]^Т[1_т, СX₁]}^–1[1_т, СX₁]^Тy

==

=

или отдельными формулами

=                                                            (11.3.11)

и

= (X₁^TС^TСX₁)^–1X₁^TС^Tу

=(X^*TX^*)^–1X^*Tу,                                         (11.3.12)

где X^*=СX₁.

Как и в (9.1.1), скорректированная усреднённым значений переменных отклика сумма квадратов S_Tс= может быть разделена в виде

S_Tс=^ТX^*Tу+[–^ТX^*Tу]

=^ТX^*Tу+S_E             [в силу (7.5.35)]

=S_Rс+S_E.                                             (11.3.13)

где S_Rс=^ТX^*Tу - образующаяся в результате регрессии и скорректированная усреднённым сумма квадратов. Подставляя в выражение для S_Rс вместо  его выражение по формуле (11.3.12), получаем S_Rс в виде квадратичной формы

S_Rс=у^ТX^*(X^*TX^*)^–1X^*Tу

Для процедуры проверки гипотез на основе имеющей распределение F статистики представим и остальные суммы квадратов выражения (11.3.13) в виде квадратичных форм относительно вектора у. Это необходимо для того чтобы, используя теоремы главы 5, показать, что суммы S_Rс и S_E имеют распределения хи-квадрат и статистически независимы. Так как сумма =y^Т(I–Е/т)y и S_E=–^ТX^*Тy, то (11.3.13) можно записать так

y^Т(I–Е/т)y=S_Rс+S_E

=у^ТX^*(X^*TX^*)^–1X^*Tу+y^Т(I–Е/т)y–у^ТX^*(X^*TX^*)^–1X^*Tу

=y^ТH^*y+y^Т(I–Е/т–H^*)y,                               (11.3.14)

где H^*=X^*(X^*TX^*)^–1X^*T.

По теореме 9.1.1 матрицы I–Е/т, H^*=X^*(X^*TX^*)^–1X^*T и I–Е/т–H^* обладают следующими свойствами: произведение H^*(I–Е/т)=H^*, матрица H^* идемпотентная и ранга р–1, матрица I–Е/т–H^* идемпотентная и ранга т–р, а произведение H^*[I–Е/т–H^*]=О.

По теореме 9.1.2, если вектор у имеет нормальное распределение N_m(СXβ, s²I), то S_Rс/s²=^ТX^*ТX^*/s² имеет нецентральное распределение c²(р–1, γ₁) с параметром не центральности γ₁=β₁^ТX^*ТX^*β₁/(2s²), а S_E/s²=[–^ТX^*ТX^*]/s² имеет центральное распределение c²(т–р). Кроме этого, по теореме 9.1.3 квадратичные формы S_Rc=y^ТH^*y и S_E=y^Т(I–Е/т–H^*)y статистически независимы.

Далее по теореме 9.1.4 статистика

F_Rc==                          (11.3.15)

приобретает следующие распределения:

Ø Если гипотеза H₀: β₁=0 ложна, то статистика F_Rc имеет нецентральное распределение F(р–1, т–р, γ₁) с параметр не центральности γ₁=β₁^ТX^*ТX^*β₁/(2s²).

Ø Если гипотеза H₀: β₁=0 верна, то параметр γ₁=0 и статистика F_Rc принимает центральное распределение F(р–1, т–р).

По этой статистике проверка гипотезы H₀: β₁=0 осуществляется так. Эта гипотеза ложна, если значение статистики F_Rc больше критического значения F_кp случайной переменной, имеющей центральное распределение F(p–1, т–p), при выбранной 1–α интегральной вероятности на интервале от 0 до F_кp.

По аналогии с таблицей 9.1.1, сведём расчёт статистики F_Rc по полученным здесь формулам в таблицу 11.3.5. В результате, проверка гипотезы H₀: β₁=0 по статистике F_Rс показывает, что полученное значение статистики F_Rс=26,41 больше критического значения F_кр=2,249 случайной переменной, имеющей центральное распределение F(7, 40) и равную 0,95 интегральную вероятность на интервале от 0 до 2,249. Поэтому гипотеза H₀: β₁=0 ложна, и надо сделать вывод, что, по крайней мере, один элемент вектора β₁ не равен нулю. Кроме этого, рекомендуется чтобы для адекватности оценки функции модели значение статистики F_Rс превосходило критическое значения F_кр в 10 раз [Box, Draper (2007) стр.281]. В данном случае F_Rс превосходит F_кр более чем в 10 раз, поэтому в результате оценки вектора параметров модели получена адекватная её функция.

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

Таблица 11.3.5. Расчёт статистике F_Rc для проверки гипотезы H₀: β₁=0

Проверка адекватности модели

Для того чтобы убедиться, что постулируемая функция модели обеспечивает приемлемую аппроксимацию истинной функции шероховатости поверхности, необходимо рассмотреть возможные доказательства адекватности модели. Проверка адекватности требует чтобы, по крайней мере, для одного набора уровней включённых в эксперимент факторов имелись реальные повторения опытов. В рассматриваемом эксперименте для каждого набора уровней факторов опыты выполнялись три раза с наблюдениями отклика, показанными в таблице 11.3.1. Повторные данные переменных отклика в виде элементов векторов у_1, у_2, у₃ получались в результате подлинных повторений опытов на разных заготовках и при переналадках станка, а не просто повторением измерений шероховатости поверхности обработанных заготовок. Имея результаты так повторённых опытов, можно проверить является ли модель (11.3.6) адекватной, используя имеющую распределение F статистику.

Для этого, как и в разделе 6.4, представим приведённые в таблице 11.3.1 данные шероховатости поверхности в виде дисперсионной модели классификации по одному признаку

у_ij=+e_ij,                                           (11.3.16)

где у_ij - значения шероховатости, полученные в опытах эксперимента (i=1, 2, ..., п и j=1, 2, 3),  - среднее переменной отклика трёх повторений i-го опыта и ошибки e_ij считается, что распределены нормально, независимо и одинаково в виде N(0, s²), где s² дисперсия ошибок. Таким образом, полагается, что значения переменных отклика в подгруппах повторных опытов являются выборками из имеющих распределения N(, s²) совокупностей наблюдений.

Рассмотренный в разделе 1.6 график кумулятивных вероятностей распределения может использоваться для проверки нормальности распределения большой выборки значений данных. В данном случае для выборки из 48 результатов опытов эксперимента построение этого графика позволяет выяснить имеют ли полученные данные эксперимента нормальное распределение. Эти данные на графике изображаются в зависимости от теоретического нормального распределения таким образом, чтобы точки данных формировали прямую линию. Отклонения от прямой линии указывают на отклонения от нормальности. На Рис.11.3.4 точки графика кумулятивных вероятностей распределения данных рассматриваемого эксперимента выстраиваются приблизительно в прямую линию. Поэтому можно считать, что эти данные имеют приблизительно нормальное распределение.

Кроме проверки, что данные эксперимента имеют нормальное распределение, перед использованием модели (11.3.16) необходимо также проверить является ли обоснованным сделанное допущение о равенстве дисперсий результатов в подгруппах повторных опытов. Для этого может использоваться, например, тест Левина [NIST/SEMATECH e-Handbook of Statistical Methods, раздел 1.3.5.10]. Этим тестом проверяется нулевая гипотеза H₀: s₁²=s₂²=…=s_п² против альтернативной гипотезы, что по крайней мере одна пара дисперсий не равны.

Рис. 11.3.4. График кумулятивных вероятностей распределения данных шероховатости поверхности (у), где Нормальные подсчёты являются шкалой стандартных отклонений от среднего.

Для выборки данных шероховатости поверхности анализируемого эксперимента объёмом 48, разделённой на 16 подгрупп с размером каждой подгруппы равным 3, тестовая статистика Левина находится по формуле

L=,                                       (11.3.17)

где δ_ij получаются из выражения δ_ij=|у_ij–|,  - усреднённое i-й подгруппы для δ_ij и - общее усреднённое для всех δ_ij.

По тесту Левина гипотеза равенства дисперсий ложна, если L > F_a(п–1, т–п,), где F_a(п–1, т–п,) - верхнее критическое значение распределения F(п–1, т–п) при уровне значимости α. Так как т=3п=48, то вычисление статистики этого теста дает L=1,611. А если взять значения 0,01, 0,05 и 0,10 уровня значимости α, то для этого теста критические значения распределения F со степенями свободы п–1=15 и т–п=32 получаются большими [L=1,611<F_0,10(15, 32) =1,707< F_0,05(15, 32) =1,992 < F_0,01(15, 32) =2,655]. Это означает, что гипотеза равенства дисперсий верна.

В матричном виде модель (11.3.16) можно записать так

y=C+e,                                                      (11.3.18)

где ^T= [_{, ,} ..., ]. С учётом графика кумулятивных вероятностей для распределения данных шероховатости поверхности и результата теста Левина можно считать, что вектор e ошибок имеет распределение N(0, Is²), где I - единичная матрица ранга т. Оценка вектора  методом наименьших квадратов выполняется по формуле

= (С^TС)^–1С^Tу=С^Tу/3.                                            (11.3.19)

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

S_PE= (у–С)^Т(у–С)

=у^T(I–H_С)у,                                                   (11.3.20)

где матрица H_С=С(С^TС)^–1С^T. В разделе 6.4 отмечено, что при проверке адекватности модели эту сумму квадратов называют суммой квадратов чистых ошибок.

Из выражения (11.3.20) видно, что S_РЕ является квадратичной формой относительно вектора у. Матрица I–H_С является идемпотентной и дисперсия D(у)=Is². Поэтому, так как S_PE/s²=у^T(I–H_С)у/s², то матрица (1/s²)(I–H_С)Is² также идемпотентная. Отсюда, так как вектор у ~ N(С, Is²), то по следствию 2 теоремы 5.5 квадратичная форма у^T(I–H_С)у имеет нецентральное распределение χ²[ранг(I–H_С)_, (С)^T(I–H_С)С/2] и

S_PE/s²~χ²[ранг(I–H_С)_, ^ТС^T(I–H_С)С/(2s²)].

В этом выражении, при ранг(I)=т и ранг(С)=п, параметр не центральности обращается в нуль и нецентральное распределение χ² сводится к центральному, так что

S_PE/s²~χ²(т–n).                                                         (11.3.21)

По теореме 5.2.1, так как Е(у) =C и D(у)=Is², а матрица (I–H_С) симметричная, то

Е(S_PE)=Е[у^T(I–H_С)у]=след[(I–H_С)Is²]+(C)^Т(I–H_С)C

=s²след(I–H_С)

=s²[след(I)–след(H_С)]           [в силу (П.11.1)]

=s²[т–ранг(H_С)]                   [по теореме П.13.4]

=s²[т–ранг(С)]                    [по пункту 5 теоремы П.8.3]

=s²(т–п).

Следовательно, MS_PE=S_PE/(т–п) даёт несмещённую оценку дисперсии s². Этот результат оценки дисперсии называют независящим от модели. Однако эта независимость распространяется только на те модели, которые можно построить на основе вектора у значений переменных отклика рассматриваемого эксперимента и для определённой матрицы С преобразования.

Для модели (11.3.6) сумма S_E квадратов остатков является суммой квадратов отклонений наблюдаемых значений отклика у от их оцениваемых ожидаемых значений =СX, то есть,

S_E=(у–СX)^T(у–СX),                                           (11.3.22)

где  находится из (11.3.7). Таким образом, подставляя (11.3.7) в (11.3.22), имеем

S_E=у^T(I–H_X)у,                                                           (11.3.23)

где Н_Х=СX(X^TС^TСX)^–1X^TС^T.

Произведение самого на себя вектора разности С–СX векторов оцениваемых ожидаемых значений шероховатости по моделям (11.3.18) и (11.3.6) известно как сумма S_LF квадратов неадекватности. Эта сумма представляется выражением

S_LF= (С–СX)^T(С–СX).                               (11.3.24)

И, подставляя в него вместо  и  их выражения из (11.3.19) и (11.3.7), получаем

S_LF =у^T(H_С–H_X)у.                                                      (11.3.25)

Здесь S_LF является квадратичной формой относительно у. Матрица H_С–Н_Х идемпотентная и D(у)=Is², поэтому матрица (1/s²)(H_С–Н_Х)Is² тоже идемпотентная. Таким образом, так как у ~ N(СXb, Is²), то по следствию 2 теоремы 5.5 квадратичная форма у^Т(Н_С–Н_X)у имеет нецентральное распределение χ²[ранг(Н_С–Н_Х), (СXb)^Т(Н_С–Н_Х)СXb/2] и, следовательно,

S_LF/s²~χ²[ранг(H_С–H_X), b^TX^TС^T(H_С–H_X)СXb/(2s²)]

~χ²[ранг(С)–ранг(X), (b^TX^TС^TСXb–b^TX^TС^TСXb)/(2s²)],                      (11.3.26)

что, при ранг(С)=n и ранг(X)=p, сводится к центральному распределению

S_LF/s²~χ²(n–p)            .                                                          (11.3.27)

По теореме 5.2.1, так как Е(у) =CXb и D(у)=Is², а матрица (H_С–Н_Х) симметричная, то

Е(S_LF)=Е[у^T(H_С–H_X)у]=след[(H_С–Н_Х)Is²]+(CXb)^Т(H_С–Н_Х)CXb

=s²след(H_С–H_X)

=s²[след(H_С)–след(H_X)]       [в силу (П.11.1)]

=s²[ранг(H_С)–ранг(H_X)]       [по теореме П.13.4]

=s²[ранг(С)–ранг(X)]             [по пункту 5 теоремы П.8.3 и пункту 1 теоремы П.4]

=s²(п–p).

Следовательно, MS_LF=S_LF/(п–p) также даёт несмещённую оценку дисперсии s². Но эта оценка зависит от конкретной рассматриваемой модели, так как выражение (11.3.25) содержит матрицу Х этой модели. Квадратичную форму у^T(H_С–H_X)у можно представить в виде у^T(H_С–H_X)у=у^T(H_С–I+I–H_X)у=у^T(I–H_X)у–у^T(I–H_С)у откуда

у^T(I–H_X)у=у^T(I–H_С)у+у^T(H_С–H_X)у.                         (11.3.28)

Следовательно, сумма S_E квадратов остатков равна сумме S_PE квадратов чистых ошибок плюс сумма S_LF квадратов неадекватности модели. А если взять математические ожидания левой и правой частей выражения (11.3.28), то имеем s²(т–р)=s²(т–п)+s²(п–p), что после сокращения s² даёт (т–р)=(т–n)+(n–p). Отсюда число степеней свободы распределения хи-квадрат квадратичной формы S_E равно сумме чисел степеней свободы распределений хи-квадрат квадратичных форм S_PE и S_LF.

В разделе 7.3 также установлено, что несмещенную оценку дисперсии s² можно сделать по формуле (7.3.8), которая в данном случае, в силу (11.3.23), принимает вид

s²=у^T(I–H_X)у/(т–р).                                                  (11.3.29)

Таким образом, все квадратичные формы выражения (11.3.28) могут быть использованы для несмещённой оценки дисперсии s², но результаты этих оценок получаются разными. Формула оценки s_pe²=S_PE/(т–п) не содержит матрицы модели, а формула оценки s_lf²=S_LF/(п–p) имеет меньше делитель, чем формула (11.3.29), и получается завышенной.

Поскольку произведения матриц (Н_С–Н_Х) и (I–Н_С) дают нулевую матрицу, то по следствию 1 теоремы 5.6.2 квадратичные формы S_LF=у^T(H_С–H_X)у и S_PE=у^T(I–H_С)у статистически независимы. А две имеющие центральные распределения χ² независимые квадратичные формы являются основой образования распределения F. То есть, если S_LF имеет распределение χ²(n–p) и S_PE имеет распределение χ²(т–n), то их отношение в виде [S_LF/(n–p)]/[S_PE/(т–n)] имеет распределение F со степенями свободы (n–p) и (т–n). Следовательно, для проверки адекватности модели имеющая это распределение статистика находится по формуле

F_LF= [S_LF/(n–p)]/[S_PE/(т–n)].                                     (11.3.30)

Расчёт статистики F_LF для проверки адекватности линейной регрессионной модели (11.3.6) для данных шероховатости поверхности [Feng (2001)] сведён в таблицу 11.3.4. И, поскольку уровень значимости α для этой проверки выбран равным 0,05, то статистика F_LF получена меньше её критического значения [F_LF=2,004<F_0,05(8, 32,) =2,244]. Это значит, что линейная регрессионная модель (11.3.6) является адекватной.

Таблица 11.3.4. Проверка адекватности линейной регрессионной модели

Анализ модели давать предсказания

В разделе 8.6 показано, что сумма квадратов остатков предсказания может быть использована для оценки способности модели давать предсказания. В силу (8.6.11), эта сумма квадратов рассчитывается по формуле

S_op==,                                                   (11.3.31)

где элементы вектора  находятся по формуле (8.6.9). Для анализируемой модели в этой формуле элементами h_ii являются элементы диагонали матрицы H_X. В данном случае все они равны 0,167. Поэтому расчёт вектора  выполняется по формуле =/(1–0,167). График зависимости обычных остатков  от остатков  показан на Рис.11.3.5. Точки на графике следуют вдоль прямой примерно с единичным наклоном. На нём нет удалённых от этой линии точек, указывающих на резко выделяющиеся значения.

Рис.11.3.5. График зависимости обычных остатков  от остатков .

На основании найденного вектора  получаем S_op=2,951. Сумма квадратов остатков предсказаний используется для вычисления статистики R_p² предсказания по формуле (8.6.12). В ней S_Tс находится из выражения (11.3.13) и статистика предсказания получается R_p²=1–S_op/S_Tс=0,744. На основании полученного результата можно ожидать, что модель (11.3.6) объяснит 74,4% изменений предсказываемых новых наблюдений шероховатости поверхности. Статистика предсказания получилась меньше статистики R²=S_Rc/S_Tс=0,822 из раздела 6.4, представляющей описываемую функцией модели часть вариации значений переменных отклика относительно их средних значений. И для рассматриваемой модели статистика R² показывает высокую меру достигнутой адекватности, так как максимально возможное значение статистики R² получается R_max²=(S_Tс–S_PE)/S_Tс=0,881.

Упражнения

11.1. В результате постановки эксперимента по плану 2⁴ по его результатам получены следующие результаты оценки воздействий факторов и их взаимодействий:

А=0,4; В= –7,6; С=14,1; D=66.7; АВ=16,7; АС=3,1; AD=5,2; BC=8,3; BD= –3,6; CD=14.3; АВС= –0,1; АBD= –4,7; ACD= 7,7; BCD= –2,3; АВСD=3,9.

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

11.2. (a) Если выполняется эксперимент по плану 2^p и каждый фактор устанавливается на уровнях –1 и +1, то какова взаимосвязь между результатами оценки воздействий факторов и коэффициентами регрессии модели первого порядка.

(б) Как в такой же ситуации взаимодействия двух фактора представляются в регрессионной модели?

11.3. Представьте, что отдельные опыты в составе факторного плана 2³ выполняются последовательно в стандартном порядке (в полном пренебрежении к рандомизации!) в то время как, одновременно, игнорируемая экспериментатором четвертая переменная принимает значения, указанные в следующей таблице.

Если для модели первого порядка у=b₀+b₁х₁+b₂х₂+b₃х₃+ε результаты оценки её параметров находятся методом наименьших квадратов, то определите смещение результатов оценки параметров модели, вызванное отсутствием члена b₄х₄ в модели.

11.4. Основная цель закалки стали это образование мартенсита. Переход из аустенита в мартенсит зависит от состава стали, который влияет на температуру образования мартенсита. Издавна считалось, что углерод и другие легирующие элементы понижают температурный диапазон, в котором образуется мартенсит. Предварительное исследование по факторному плану 2³ для определения влияний углерода, марганца и никеля на ‘температуру начала образования мартенсита’ дало данные в следующей таблице.

(а) Сделайте факторы нормированными.

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

(в) Найдите все результаты оценки воздействий факторов, двухфакторных и трехфакторных взаимодействий.

(г) Допустим, что стандартное отклонение наблюдений отклика составляет s=65. Которые результаты оценки в (в) большие по сравнению с их стандартными ошибками?

(д) Каковы практические выводы из вашего анализа?

11.5. В приведённой ниже таблице данные получены в результате исследования влияния трех факторов: ξ₁ - обороты в минуту (об/мин), ξ₂ - осевая нагрузка вниз (вес, фунт) и ξ₃ - скорость потока жидкости для сверления, на скорость (у) врезания сверла. Опыты представлены в случайном порядке, в котором они проводились. Пронормируйте факторы и определите тип используемого плана. Затем сделайте соответствующий анализ и перечислите полученные вами практические выводы.

11.6. Изготовитель бумажных пакетов желает проверить сопротивление (у) разрыву, для чего давно используется числовая шкала. Рассматриваются три фактора на каждом из двух нормированных уровнях, а именно: х₁ - тип бумаги, х₂ - влажность, х₃ - направление разрыва, и получаются три повторных наблюдения, как показано в приведённой ниже таблице. (Для удобства опыты переставлены в стандартный порядок из случайного, в котором они выполнялись.) Общая сумма квадратов 496,87 для 24 наблюдений. Выполните факторный анализ данных, оцените s² на основе результатов повторных опытов и сформулируйте выводы о факторах.

11.7. Посмотрим еще раз на данные в упражнении 11.5. Оцените параметры модели первого порядка у=b₀+b₁х₁+b₂х₂+b₃х₃+ε методом наименьших квадратов и выполните обычный анализ, включая анализ адекватности модели. Ваши выводы должны согласовываться с данными факторного анализа.

11.8. Посмотрим на данные упражнения 11.6. Мы могли бы формально использовать уравнение регрессии для анализа этих данных, но они гораздо более подходяще анализируются факторными методами. Почему?

11.9. Результаты оценки воздействий, полученные на основе данных выполненного по плану 2⁵ эксперимента, следующие:

Е=–224             BCD=–18       BE=29

С=–153             ABDE=–14    DЕ=30

ABCDE=–77    D=–9              ABCE=31

ACE=–58         B=–6               BCE=39

AD=–54                       CD=–4                       ACDE=47

BC=–53                        АВС=0           АС=53

ABD=–34         АE=2              АВСD=58

ACD=–33         BD=7              АВ=64

BDE=–28         CDE=12         СЕ=83

ABE=–22         BCDE=16      А=190

АDE=21

(а) Для этих результатов постройте график кумулятивных вероятностей их распределения.

(б) Какие воздействия вы думаете, отличаются от шума?

(в) Постройте график кумулятивных вероятностей нормального распределения оставшихся воздействий.

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

[Частичный ответ: (б) E, C, A. (г) грубо 40]

11.10. (Источник: Wasiloff Е., Hargitt С. Using DOE to determine AA battery life, Quality Progress, Март 1999 г., cтр.67-71.) Срок службы двух типов батареек для игрушечных машин проверялся путём постановки факторного эксперимента по плану 2³. Три фактора и их уровни были следующими:

1. Тип батареи: низкой стоимости (х₁=–1), высокой стоимости (х₁=+1).

2. Контакты разъема: стандартный Тамия (х₂=–1), Дин позолоченный (х₂=+1).

3. Начальная температура батареи: холодная (х₃=–1), окружающей среды (х₃=+1).

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

Проанализируйте данные и сделайте выводы о том, как срок службы батареи зависит от трех факторов и их взаимодействий.