Планирование эксперимента и обработка результатов

Федеральное агентство по образованию Государственное образовательное учреждение высшего профессионального образования «Национальный Исследовательский Томский политехнический университет» ______________________________________________________________ В.Н. Шкляр Планирование эксперимента и обработка результатов Конспект лекций для магистров по направлению 220200 «автоматизация и управление в технических (мехатронных) системах» Издательство Томского политехнического университета 2010 1 Лекция 1 Планирование эксперимента и обработка результатов 1. 2. 3. 4. 5. 6. Литература Планирование эксперимента в исследовании технологических процессов / Под ред. Э.К.Лецкого. - М.: Мир, 1977. Математическая статистика с элементами теории планирования эксперимента: учеб. пособие: Л.В. Горская, В.Н. Пиунова, В.С. Смирнова. Саратовский полит. институт, 1975, 103 с. Математическая теория планирования эксперимента./ Под ред. С.М. Ермакова.М.: Наука. Главная редакция физико-математической литературы, 1983. 392 с. Надёжность технических систем: Справочник/ Ю.К. Беляев, В.А. Богатырёв и др. Под ред. И.А. Ушакова М.: Радио и связь, 1985. 608 с. Дьяконов В.П. MATLAB 6/6.1/6.5 +Simulink 4/5 в математике и моделировании. Полное руководство пользователя.- М.: СОЛОН-Пресс, 2003. 576 с. ГОСТ 24026-80 « Исследовательские испытания. Планирование эксперимента. Термины и определения. М.,1980. Дополнительная 1. Исследования и изобретательство в машиностроении. Практикум/ Под общ. Ред. М.М. Кане. Мн.: УП «Технопринт», 2003. 237 с. 2. Адлер Ю.П. Введение в планирование эксперимента. М., 1969. Минько А.А. Статистический анализ в MS Excel. – М.: Изд. дом «Вильямс», 2004. – 448 с. 3. Тюрин Ю.Н. Анализ данных на компьютере / Ю.Н.Тюрин, А.А. Макаров. - М.: Финансы и статистика, 1995. Основные понятия и определения Эксперимент - система операций, воздействий и (или) наблюдений, направленных на получение информации об объекте при исследовательских испытаниях Опыт- воспроизведение исследуемого явления в определенных условиях проведения эксперимента при возможности регистрации его результатов План эксперимента - совокупность данных, определяющих число, условия и порядок реализации опытов Планирование эксперимента - выбор плана эксперимента, удовлетворяющего заданным требованиям. Фактор (ндп. Параметр) - переменная величина, по предположению влияющая на результаты эксперимента. Уровень фактора - фиксированное значение фактора относительно начала отсчета Основной уровень фактора - натуральное значение фактора, соответствующее нулю в безразмерной шкале. Нормализация факторов - преобразование натуральных значений факторов в безразмерные значения. 2 Априорное ранжирование факторов - метод выбора наиболее важных факторов, основанный на экспертной оценке Размах варьирования фактора - разность между максимальным н минимальным натуральными значениями фактора в данном плане Интервал варьирования фактора - половина размаха варьирования фактора Эффект взаимодействия факторов - показатель зависимости изменения эффекта одного фактора от уровней других факторов. Факторное пространство - пространство, координатные оси которого соответствуют значениям факторов. Область экспериментирования (область планирования) - область факторного пространства, где могут размещаться точки, отвечающие условиям проведения опытов. Активный эксперимент - эксперимент, в котором уровни факторов в каждом опыте задаются исследователем. Пассивный эксперимент - эксперимент, при котором уровни факторов в каждом опыте регистрируются исследователем, но не задаются. Последовательный эксперимент (ндп. шаговый эксперимент) - эксперимент, реализуемый в виде серий, в котором условия проведения каждой последующей серии определяются результатами предыдущих. Отклик (ндп. Реакция на Параметр ) - наблюдаемая случайная переменная, по предположению, зависящая от факторов. Функция отклика - зависимость математического ожидания отклика от факторов. Оценка функции отклика - зависимость, получаемая при подстановке в функцию отклика оценок значений ее параметров. Дисперсия оценки функции отклика - дисперсия оценки математического ожидания отклика в некоторой данной точке факторного пространства. Поверхность отклика (ндп. Поверхность регрессии) - геометрическое представление функции отклика. Поверхность уровня функции отклика - геометрическое место точек в факторном пространстве, которому соответствует некоторое фиксированное значение функции отклика. Область оптимума - область факторного пространства в окрестности точки, в которой функция отклика достигает экстремального значения. Рандомизация плана - один из приемов планирования эксперимента, имеющий целью свести эффект некоторого неслучайного фактора к случайной ошибке. Параллельные опыты - рандомизированные во времени опыты, в которых уровни всех факторов сохраняются неизменными. Временный дрейф - случайное или неслучайное изменение функции отклика во времени. Исследователь на этапе планирования эксперимента должен: · помнить, к какому классу относится моделируемая система (статическая или динамическая, детерминированная или стохастическая и т.д.); · определить, какой режим работы его интересует, стационарный (установившийся) или нестационарный; · знать, в течение какого промежутка времени следует наблюдать за поведением (функционированием ) системы; 3 знать, какой объём испытаний (то есть повторных экспериментов) сможет обеспечить требуемую точность оценок (в статистическом смысле) исследуемых характеристик системы. Таким образом, планирование модельных экспериментов преследует 2 основные цели: 1. Сокращение общего объёма испытаний при соблюдении требований к достоверности и точности их результатов; 2. повышение информативности каждого из экспериментов в отдельности. Поиск плана эксперимента проводится в факторном пространстве. Факторное пространство (ФП) – это множество внешних и внутренних параметров модели, значения которых исследователь может контролировать в ходе подготовки и проведения эксперимента. Значения факторов обычно называются уровнями. Если при проведении эксперимента исследователь может изменять уровни факторов, эксперимент называется активным, в противном случае – пассивным. Введём ещё несколько терминов, используемых в теории планирования эксперимента. Каждый из факторов имеет верхний и нижний уровни, расположенные симметрично относительно некоторого нулевого уровня. Точка в ФП, соответствующая нулевым уровням всех факторов, называется центром плана. Интервалом варьирования фактора называется некоторое число J, прибавление которого к нулевому уровню даёт верхний уровень, а вычитание – нижний. Как правило, план эксперимента строится относительно одного (основного) выходного скалярного параметра Y, который называется наблюдаемой переменной. Если моделирование используется как инструмент принятия решения, то в роли наблюдаемой переменной выступает показатель эффективности. При этом предполагается, что значение наблюдаемой переменной, полученное в ходе эксперимента, складывается из двух составляющих: Y = f ( x) + e( x) , · где f ( x ) - функция отклика (неслучайная функция факторов), e ( x ) - ошибка эксперимента (случайная величина). Дисперсия DY наблюдаемой переменной, которая характеризует точность измерений, равна дисперсии ошибки опыта, т.е. DY = De . DY называют дисперсией воспроизводимости эксперимента. Она характеризует качество эксперимента. Эксперимент называется идеальным при DY = 0 . Существует два основных варианта постановки задачи планирования имитационного эксперимента: · Из всех допустимых требуется выбрать такой план, который позволил бы получить наиболее достоверное значение функции отклика f ( x ) при фиксированном числе опытов; 4 Из всех допустимых требуется выбрать такой план, при котором статистическая оценка функции отклика может быть получена с заданной точностью при минимальном объёме испытаний. Решение задачи планирования в первой постановке называется стратегическим планированием, во второй – тактическим планированием. При проведении опытных исследований различают пассивный и активный эксперимент. Методология пассивного экспериментирования предполагает проведение большой серии опытных исследований с поочередным варьированием значений входных переменных x и анализом результатов измерений выходной переменной y (лабораторный эксперимент или эксперимент на пилотной установке). К пассивному эксперименту принято относить также и сбор опытных данных в режиме эксплуатации промышленной установки – т.н. промышленный эксперимент. Обработка результатов, и выбор вида эмпирической модели (уравнения регрессии), т.е. решение задачи структурной идентификации является достаточно сложной задачей. Это связано с тем, что вид уравнения регрессии необходимо определять по характеру изменения переменных на графике эмпирической линии регрессии, полученной по выборке экспериментальных данных. Для решения этой задачи для одной входной переменной x предложены эффективные методы, в которых предусматривается преобразование системы координат как для входной ( x ), так и для выходной переменной ( y ). При большем числе входных переменных ( x1,…xm ) надёжных методов определения вида уравнения регрессии (вида эмпирической модели) в настоящее время не существует. Активный эксперимент проводится по заранее составленному плану, в соответствии с которым ставится задача не только определения оптимальных условий проведения эксперимента, но и оптимизации процесса (оптимальное планирование эксперимента). При определении оптимальных условий проведения процесса с использованием эмпирических моделей (например, методом Бокса-Вильсона) выходная переменная ŷ является критерием оптимальности или целевой функцией. В теории активного экспериментирования выходную (зависимую) переменную принято называть функцией отклика, а входные (независимые) переменные – факторами. Соответственно - координатное пространство с координатами ( x1, x2, …xm ) - факторным пространством, а геометрическое изображение функции отклика в факторном пространстве – поверхностью отклика. Активный эксперимент планируется таким образом, чтобы упростить обработку его результатов методами регрессионного и корреляционного анализа. К достоинствам активного экспериментирования относятся: · возможность предсказания количества опытов, которые следуют провести; · определение точек факторного пространства, где следует проводить опыты; · отсутствие проблем, связанных с выбором вида уравнения регрессии; · 5 · возможность определения оптимальных параметров процесса экспериментально-статистическим методом; · сокращение объёма опытных исследований. Основы планирования многофакторного эксперимента В общем случае объект исследования можно представить в виде структурной схемы, показанной на рис.1. x1 x2 .. . xk wn . .. w1 w2 Объект y1 y2 .. . yp Рис.1 Представление объекта в виде такой схемы основано на принципе «черного ящика». Имеем следующие группы параметров: 1) управляющие (входные) xi которые называются факторами; 2) выходные параметры yi, которые называются параметрами состояния; 3) wi - возмущающие воздействия. Предполагается, что возмущающие воздействия не поддаются контролю и либо являются случайными, либо меняются во времени. Каждый фактор xi имеет область определения, которая должна быть установлена до проведения эксперимента. Комбинацию факторов можно представить как точку в многомерном пространстве, характеризующую состояние системы. На практике целью многофакторного эксперимента является установление зависимости y = f (x1 , x2 , ... xk ) , (1) описывающей поведение объекта. Чаще всего функция (1) строится в виде полинома y = a0 + a1 x1 + a2 x2 (2) или y = a 0 + a1 x1 + a 2 x 2 + a11 x12 + a 22 x 22 + a12 x1 x 2 . (3) Целью эксперимента может быть, например, построение зависимости (1) при минимальном количестве измерений значений управляющих параметров xi . На первом этапе планирования эксперимента необходимо выбрать область определения факторов xi . Выбор этой области производится исходя из априорной информации. Значения xi называются уровнями управляющего параметра. 6 Если выбрана линейная модель (2), то для построения аппроксимирующей функции достаточно выбрать основной уровень и интервал варьирования управляющего параметра xi . Для линейной модели интервал варьирования можно определить как x - xmin I = max 2 , а основной (нулевой) уровень - как среднее значение x0 = x max + x min 2 . Для упрощения планирования эксперимента принято вместо реальных натуральных) уровней xi использовать кодированные значения факторов. Для факторов с непрерывной областью определения это можно сделать при помощи следующего преобразования (нормирования факторов) xj = ~ x j - x j0 Ij , ~ x где x j - натуральное значение фактора; I j - интервал варьирования; j 0 - основной уровень; x j - кодированное значение. В результате x j принимает значения x = ±1 x =0 на границах j , на основном уровне j . Основная проблема состоит в выборе области варьирования, поскольку эта задача является неформализованной. Рассмотрим полный факторный эксперимент на примере линейной модели (2). Если число факторов k , то для проведения полного факторного эксперимента нужно N = 2 опытов, где 2 - число уровней, которого достаточно для построения линейной модели. Условие проведения этого эксперимента можно зафиксировать в матрице планирования (табл.1). k Таблица 1 Номер y x1 x2 опыта 1 -1 -1 y1 2 +1 -1 y2 3 -1 +1 y3 4 +1 +1 y4 Таким образом, для двух факторов построение матрицы планирования элементарно. Для большего числа факторов необходимо разработать правила построения таких матриц. Например, при появлении фактора x3 в табл.1 произойдут следую7 щие изменения (табл.2): при появлении нового столбца каждая комбинация уровней исходной таблицы проявится дважды. Таблица 2 Номер x3 y x1 x2 опыта 1 -1 -1 +1 y1 2 +1 -1 +1 y2 3 -1 +1 +1 y3 4 +1 +1 +1 y4 5 -1 -1 -1 y5 6 +1 -1 -1 y6 7 -1 +1 -1 y7 8 +1 +1 -1 y8 Это не единственный способ расширения матрицы планирования. Используют также перемножение столбцов, правило чередования знаков. Очень важны общие свойства матрицы планирования: · симметричность матрицы относительно центра эксперимента: xi = 0 . Тогда N åx i =1 · N 2 å xij = N - ji =0 условие нормировки, то есть сумма квадратов элементов каждого i =1 столбца равна числу опытов. Первые два свойства относятся к построению отдельных столбцов матрицы. · N å xij × xin = 0 - совокупность столбцов имеет следующее свойство, где j ¹ n . i =1 · Ротатабельность. Это означает, что точки (значения факторов) в матрице планирования подбираются так, что точность предсказания выходного параметра должна быть одинакова на равных расстояниях от центра эксперимента (нулевого уровня) и не зависеть от направления. Планирование эксперимента первого порядка для двух переменных. План эксперимента первого порядка для двух переменных показан на рис. 2. То есть искомая функция y = f ( x1 , x2 ) описывается модельно в виде плоскости y = a0 + a1 x1 + a2 x2 (4) или гиперболоида y = a 0 + a1 x1 + a 2 x 2 + a3 x1 x 2 . (5) 8 Расположение этой модели в пространстве показано на рис.2 поверхностью, проходящей через точки 1 – 2 – 3 – 4. Рис.2 Необходимые уровни для полного факторного эксперимента расположены в плоскости ( x1 , x2 ) . Для модели в виде гиперболоида этот план является предельно экономным. Для построения гиперболоида необходимо определить четыре коэффициента в модели (5). Это можно сделать, решая систему из четырех уравнений. Следовательно, необходимы все четыре опыта. В теории планирования эксперимента используется термин насыщенности. Если рассматривать модель (4) в виде плоскости, то план эксперимента является ненасыщенным (избыточным), так как необходимо определить только три коэффициента a0 , a1 и a 2 . В случае модели (5) (насыщенный эксперимент) решение системы единственно, и поверхность гиперболоида пройдет через все четыре экспериментальных значения yi . Следствием этого является то, что насыщенный эксперимент не позволяет усреднить случайные погрешности и не дает сведения об их размере. Для ненасыщенного плана (4) избыточное число опытов позволяет произвести усреднение и оценить размеры погрешности. Проведя плоскость через точки 1, 2 и 3, можно оценить погрешность, определив, на каком расстоянии от плоскости находится точка 4. Погрешность в других точках может быть оценена проведением плоскостей 1 – 3 – 4, 1 – 2 – 4 и 2 – 3 – 4. С другой стороны коэффициент a1 наклона поверхности к оси x1 может быть найден как из наклона прямой 1 – 2, так и из наклона прямой 3 – 4. Аналогично коэффициент a2 при x2 можно определить из наклона прямых 1 – 3 и 2 – 4. 9 Поскольку полученные таким образом значения a1 и a2 могут отличаться, ненасыщенный эксперимент позволяет провести их усреднение и оценить погрешность. Если уравнение плоскости представить в виде y = a0 + a1 ( x1 - x1 ) + a 2 ( x2 - x2 ) , (6) x1min + x1max x + x2 max x2 = 2 min 2 2 ; , то мы переносим начало коордигде нат в точку с координатами ( x1 , x2 ) . Тогда коэффициент a0 находится усреднениx1 = ем всех четырех значений yi как высота центра плоскости 1 – 2 – 3 – 4. Процесс переноса начало координат в центр пространства факторов с коорди- натами ( x1 , x2 ,K, xk ) очень важен при обработке данных любых экспериментов, описываемых моделью в виде гиперплоскости, так как позволяет получить более устойчивое усредненное значение для a0 . Важнейшим фактором является то, что в результате такого усреднения построенная плоскость удовлетворяет всем четырем значениям yi лишь в среднем. В любой точке может быть найдена погрешность отклонения экспериментальных данных относительно модели, и по этим четырем отклонениям можно вычислить СКО. Таким образом, один из четырех опытов является избыточным и может быть исключен. Но тогда план эксперимента становится неротатабельным, то есть неравноточным по всем направлениям. Если исключена точка 4 на рис.2, то в направлении 3 – 2 в плоскости факторов будет обеспечена большая точность, чем в направлении 1 - 0. В этом случае для восстановления ротатабельности точки 1, 2 и 3 в плоскости факторов должны быть равноудалены как друг от друга, так и от центра, то есть располагаться в вершинах равностороннего треугольника с центром в точке 0. В общем случае для линейной модели (4), эксперимент содержащий конечное число опытов позволяет получить только оценки для коэффициентов a0 , a1 и a2 . Подставив в уравнение модели (4) известные значения факторов xij и результаты опытов yi получим систему линейных алгебраических уравнений для определения ai . Если количество этих уравнений больше трех, то значения оценок a0 , a1 и a2 могут быть получены при помощи МНК: æ N ö a j = ç å xij yi ÷ N è i =1 ø , (7) где N - количество опытов. Здесь учтено, что xij принимают значения -1,+1. Для вычисления коэффициентов линейной модели по формуле (7) получим: 10 a1 = [(-1) y1 + (+1) y 2 + (-1) y 3 + (+1) y 4 ] 4 , (8) [(-1) y1 + (-1) y 2 + (+1) y3 + (+1) y 4 ] a2 = 4 Таким образом, для вычисления a1 и a2 можно использовать (8). Для определения a0 в формуле (6.2.4) найдем среднее значение 1 n y= å yi n =i 1 , равное 1 n 1 n x = x1i x 2 = å x2i y = a0 + a1 x1 + a2 x2 , где 1 n å n i =1 i =1 , . В случае симметричности матрицы планирования x1 = x2 = 0 , откуда y = a0 . Чтобы коэффициент модели вычислялся по единой формуле (7) в матрице планирования вводят фиктивную переменную x0 , которая принимает значение 1 во всех опытах и соответствует коэффициенту a0 . Коэффициент при независимых переменных xi указывает на силу влияния факторов: чем больше значение имеет коэффициент ai , тем большее влияние оказывает соответствующий фактор. В этом смысле результаты планирования эксперимента аналогичны факторному анализу. Для пассивных экспериментов факторный анализ может использоваться в качестве априорных данных при планировании. Планируя эксперимент, стремятся получить линейную модель, однако в выбранных интервалах варьирования априори не известно, что линейная модель адекватно описывает поведение системы. Нелинейность связана со смешанным взаимодействием. Формула (5) всегда может быть оценена по полному факторному эксперименту. Для полного факторного эксперимента N = 2 приведена в табл.3. 2 матрица планирования с учетом эффекта взаимодействия Таблица 3 Номер опыта 1 x0 x1 x2 x1 x2 y +1 -1 -1 +1 y1 2 +1 +1 -1 -1 y2 3 +1 -1 +1 -1 y3 4 +1 +1 +1 +1 y4 В этом случае коэффициент a12 также может быть вычислен по формуле (9): 11 a12 = [(+1) y1 + (-1) y2 + (-1) y3 + (+1) y4 ] 4 . (9) Столбцы x1 , x2 задают планирование эксперимента – по ним определяют результаты опыта; столбцы x0 , x1 x2 служат только для расчета. С ростом числа факторов число возможных взаимодействий возрастает. Напри3 мер, для факторного эксперимента N = 2 кроме x0 , x1 , x2 , x3 в матрице плани- рования появляются столбцы x1 x2 , x1 x3 , x2 x3 , x1 x2 x3 . Всего в матрице планирования оказывается восемь столбцов, следовательно, необходимо определять восемь коэффициентов. Все восемь коэффициентов необходимо определять в том случае, если учитывать смешанное взаимодействие. Если же модель задается в виде гиперплоскости (линейная модель), то достаточно определить четыре коэффициента: x0 , x1 , x2 , x3 . Полный факторный эксперимент оказывается избыточным и у экспериментатора возникает выбор: 1. Построить гиперплоскость по четырем экспериментам, а остальные четыре опыта использовать для оценки погрешности. 2. Провести эксперимент, состоящий из 4-х опытов, то есть реализовать экономный план эксперимента. Таким образом, в отличие от модели гиперболоида, которая требует определеk ние 2 неизвестных коэффициентов, модель гиперплоскости, содержит k + 1 коэффициент и требует соответствующего числа опытов, то есть полный факторный план (ПФП) для модели гиперплоскости сильно избыточен. Для построения гиперплоскости, следовательно, достаточно использовать лишь некоторую часть из ПФП. Эту часть в теории планирования эксперимента называют дробной репликой или дробным факторным планом (ДФП). Если дробление ПФП производится последовательным делением числа опытов на 2, то реплику называют регулярной. Число p последовательного деления называют дробностью реплики. Число опытов регулярного ДФП равняется n = 2 k- p . При p = 1 ДФП называ- ют полурепликой (или 1/2 реплика), при p = 2 – 1/4 реплика и т.д. Соответствующее число опытов и параметров планирования приведены в таблице 4. Таблица 4 Число Число коэфф. Число опыфакторов, модели, Вид плана тов ПФП k k +1 Число плана 2 3 4 5 6 3 4 5 6 7 4 4 8 8 8 4 8 16 32 64 ПФП Полуреплика Полуреплика Четвертьреплика 1/8 реплика опытов Избыточность 1 3 2 1 12 7 8 8 9 128 256 1/16 реплика 1/16 реплика 8 16 7 Для составления планов-таблиц регулярных дробных реплик часто используют так называемое правило двоичного кода. Оно гласит, что для модели в виде гиперболоида знаки “+” и “–“ в столбцах плана должны чередоваться по правилу чередования двоичных чисел в разряде двоичного кода, то есть в столбце x1 - через 1, в k -1 столбце x2 - через 2, в столбце x3 - через 4, в столбце xk - через 2 . Лекция 2 Построение эмпирических моделей по данным активного эксперимента При проведении опытных исследований различают пассивный и активный эксперимент. Методология пассивного экспериментирования предполагает проведение большой серии опытных исследований с поочередным варьированием значений входных переменных x и анализом результатов измерений выходной переменной y (лабораторный эксперимент). К пассивному эксперименту принято относить также и сбор опытных данных в режиме эксплуатации промышленной установки – т.н. промышленный эксперимент. Обработка результатов пассивного эксперимента проводится методами регрессионного и корреляционного анализа, и выбор вида эмпирической модели (уравнения регрессии), т.е. решение задачи структурной идентификации является достаточно сложной задачей. Это связано с тем, что вид уравнения регрессии необходимо определять по характеру изменения переменных на графике эмпирической линии регрессии, полученной по выборке экспериментальных данных. Для решения этой задачи для одной входной переменной x предложены эффективные методы, в которых предусматривается преобразование системы координат как для входной ( x ), так и для выходной переменной ( y ). При большем числе входных переменных ( x1,…xm ) надёжных методов определения вида уравнения регрессии (вида эмпирической модели) в настоящее время не существует. Активный эксперимент проводится по заранее составленному плану, в соответствии с которым ставится задача не только определения оптимальных условий проведения эксперимента, но и оптимизации процесса (оптимальное планирование эксперимента). При этом уравнения регрессии (эмпирические модели) описывают данные активного эксперимента, в основном, в двух ограниченных областях и имеют следующий вид: - вдали от экстремального значения выходной переменной y : 13 m m -1 m j =1 j =1 u = 2 yˆ I = a0 + å a j x j + åå a ju x j xu u> j - вблизи экстремального значения выходной переменной y («в почти стационарной области»): m m -1 m m yˆ II = a0 + å a j x j + åå a ju x j xu + å a j x 2j . j =1 j =1 u = 2 j =1 u> j Приведённые уравнения являются линейными относительно коэффициентов регрессии a и имеют достаточно простой вид. Они включают слагаемые с двойным взаимодействием входных переменных æ m -1 m ö ç åå a ju x j xu ÷ ÷ ç j =1 u = 2 è ø u> j и не учитывают взаимодействия более высоких порядков (тройные, четверные и т.д.), вероятность которых существенно меньше. Последнее уравнение включает слагаемые с квадратами входных переменных ö æ m ç å a j x 2j ÷ ÷ ç ø è j =1 и его коэффициенты получаются при обработке результатов активных эксперименII тов II-го порядка (верхний индекс II при yˆ : yˆ ) – например, ОЦКП – ортогонального центрального композиционного плана эксперимента. Предпоследнее уравнение не включает слагаемые с квадратами входных переменных и его коэффициенты получаются при обработке результатов активных экспериI ментов I-го порядка – верхний индекс I при yˆ : yˆ - например, ПФЭ – полный факторный эксперимент. При определении оптимальных условий проведения процесса с использованием эмпирических моделей (например, методом Бокса-Вильсона) выходная переменная ŷ является критерием оптимальности или целевой функцией. В теории активного экспериментирования выходную (зависимую) переменную принято называть функцией отклика, а входные (независимые) переменные – факторами. Соответственно - координатное пространство с координатами ( x1, x2, …xm ) - факторным пространством, а геометрическое изображение функции отклика в факторном пространстве – поверхностью отклика. Активный эксперимент планируется таким образом, чтобы упростить обработку его результатов методами регрессионного и корреляционного анализа. 14 Ортогональные планы экспериментов, используемые при активном экспериментировании, обеспечивают диагональный вид корреляционной матрицы C при регрессионном анализе и, соответственно, статистическую независимость коэффициентов регрессии. К другим достоинствам активного экспериментирования относятся: · возможность предсказания количества опытов, которые следуют провести; · определение точек факторного пространства, где следует проводить опыты; · отсутствие проблем, связанных с выбором вида уравнения регрессии; · возможность определения оптимальных параметров процесса экспериментально-статистическим методом; · сокращение объёма опытных исследований. 1. ПФЭ и обработка его результатов. Полный факторный эксперимент (ПФЭ) относится к экспериментам I-го порядI ка, т.к. описывающее его уравнение ŷ не включает факторы в квадрате. Для двух факторов ( x1 и x2 ) и без учёта взаимодействия факторов соответствующая эмпирическая модель может быть записана: yˆ I = a0 + a1 x1 + a2 x2 В соответствии с теорией ПФЭ при проведении опытных исследований каждый из факторов варьируется только на двух уровнях – минимальном (кодированное значение -1) и максимальном (кодированное значение +1). При этом реализуются возможные комбинации минимальных и максимальных значений факторов, в результате чего общее число опытов (n) в ПФЭ равно 2m и полный факторный эксперимент обычно называется ПФЭ типа 2m . Для определения числа опытов применяется формула: n = 2m В последнее уравнение включаются кодированные значения факторов zj вместо xj , значения которых получаются по следующей схеме кодирования: zj = x j - x (j0 ) Dx j j = 1,...m , где ( x (j0 ) = 0.5 x min + x max j j Dx j = x max - x min j j 2 , ) j = 1,...m В результате план проведения эксперимента с учётом вышесказанного и кодирования факторов имеет вид: (число факторов равно 2 - m = 2, число опытов n = 2m = 22 = 4) 15 p z0 z1 z2 y' 1 +1 -1 -1 y 1' 2 +1 +1 -1 y 2' 3 +1 -1 +1 y 3' 4 +1 +1 +1 y 4' n При этом уравнение регрессии, описывающее эти опытные данные, записывается с использованием кодированных факторов zj ( j = 0, 1, 2 ) и соответственно кодиро~ ~ ~ ванных коэффициентов регрессии a0 , a1 , a2 : yˆ = a~ + a~ z + a~ z . 1 1 2 2 В кодированном факторном пространстве в соответствии с указанным планом проведения эксперимента проведённые опыты представляются точками вершин квадрата: z2 (+ 1,+1) (- 1,+1) (0,0) (- 1,-1) z1 (+ 1,-1) Для параметрической идентификации кодированного уравнения регрессии используется метод регрессионного анализа, включающий три этапа: · определение кодированных коэффициентов уравнения регрессии a~ методом наименьших квадратов; · оценка значимости кодированных коэффициентов регрессии с использованием t – критерия Стьюдента; · проверка адекватности кодированного уравнения регрессии с использованием F – критерия Фишера. Реализация двух последних этапов возможна при выполнении свойства однородности дисперсий (одно из требований регрессионного анализа) и проведении параллельных опытов, например, в точке с координатами z1 = 0 и z2 = 0 (центр плана, на рисунке – тёмная точка). При проведении k параллельных опытов в центре плана ( y 0'S , s = 1,...k ) среднее ' значение yc определяется как среднее арифметическое результатов измерений во всех параллельных опытах: k yc' = åy S =1 k ' 0S . 16 2. Определение кодированных коэффициентов регрессии (ПФЭ) В этом случае используется применяемая при линейном регрессионном анализе матричная формула метода наименьших квадратов (МНК), которая с учётом кодирования факторов имеет вид: æ ~ T ~ ö a~ = çç F F ÷÷ ( m +1)´1 è ( m+1)´n n´( m +1) ø -1 ~ F T ' y , (m +1)´n n´1 где кодированная матрица, зависящая от независимых переменных для двух факторов включает только +1 и -1 и имеет вид: é z10 êz ~ F = z = ê 20 ( 4´3 ) ê z30 ê ë z40 z12 ù é + 1 z22 úú êê + 1 = z32 ú ê + 1 ú ê z42 û ë + 1 z11 z21 z31 z41 - 1 - 1ù + 1 - 1úú - 1 + 1ú ú + 1 + 1û Матрица z при активном экспериментировании называется матрицей планирования и обладает тремя оптимальными свойствами: · симметричности: сумма элементов всех столбцов матрицы, кроме первого (точнее, нулевого) равна нулю n åz i =1 ij = 0, j = 1,...m; · ортогональности: скалярное произведение двух любых столбцов матрицы равно нулю n z Tj zu = å zij ziu = 0 i =1 j , u = 0,1,...m u ¹ j ; · нормировки: скалярное произведение двух одинаковых столбцов матрицы равно n ( n = 2m в ПФЭ) n z z j = å zij2 = n T j i =1 j = 0,1,...m . Благодаря перечисленным оптимальным свойствам матрицы планирования формационная матрица в ПФЭ при m=2 равна z ин- 17 ~ ~T ~ I = F F = z (3´3 ) (3´4 ) (4´3 ) T (3´4 ) én 0 0ù z = ê 0 n 0 úú, ( 4´3 ) ê êë 0 0 n úû т.е. она является диагональной с одинаковыми элементами на главной диагонали, равными n=22=4. Соответственно, корреляционная матрица одинаковыми элементами главной диагонали: ~ æ~ ~ö C = çç F F ÷÷ (3´3) ø è T -1 C также будет диагональной и с é n -1 ê æ T ö = çz z÷ = ê 0 ø è ê0 ë -1 -1 n 0ù ú 0ú n -1 úû Результатом подстановки последних соотношений в матричную формулу для определения кодированных коэффициентов регрессии будет простая формула: n a~ j = åz ij yi' , j = 0,1,...m n При учёте взаимодействия двух факторов z1 и z2 кодированное уравнение регрессии принимает вид: i =1 yˆ = a~0 + a~1 z1 + a~2 z 2 + a~12 z1 z 2 и в матрицу планирования z включается ещё один дополнительный последний столбец, каждый элемент которого равен произведению элементов столбцов, соответствующих взаимодействующим факторам: é z10 êz ~ F = z = ê 20 ( 4´4 ) ê z30 ê ë z40 z11 z 21 z12 z22 z31 z32 z 41 z42 (z11 z12 )ù é+ 1 ( z21 z22 )úú êê+ 1 = ( z31 z32 )ú ê+ 1 ( z41 z42 )úû êë+ 1 - 1 - 1 + 1ù + 1 - 1 - 1úú - 1 + 1 - 1ú ú + 1 + 1 + 1û При этом матрица планирования сохраняет все три оптимальных свойства – симметричности, ортогональности и нормировки, а кодированный коэффициент уравнения регрессии при члене, характеризующем взаимодействие факторов, определяется по формуле: å (z z )y n a~ ju = ij iu ' i , j, u = 1,...m u > j n В теории ПФЭ доказывается, что при увеличении числа факторов ( m >2 ) матриi =1 z ца планирования (n´ p ) строится с использованием рассмотренной методики, в том числе и с учётом взаимодействия факторов (не только двойного, но и тройного, четверного и т.д.). 18 В этом случае число столбцов матрицы p зависит от числа учёта взаимодействий факторов n = 2m и матрица планирования сохраняет перечисленные оптимальные свойства. Поэтому для определения кодированных коэффициентов регрессии используются приведённые выше формулы. Для расчёта натуральных значений коэффициентов в кодированное уравнение регрессии вместо кодированных факторов zj ( j = 1, … m ) следует подставить выражения для последних через натуральные значения факторов x j ( j = 1, … m ) в соответствии с приведённой выше схемой кодирования. 3. Определение значимости кодированных коэффициентов регрессии (ПФЭ) Незначимость кодированных коэффициентов регрессии определяется с использотабл ванием квантиля t – распределения Стьюдента t b ( f ) при помощи неравенства: a~ j £ t bтабл ( fe ), S a~ j где β – доверительная вероятность (в инженерных расчётах равная 0,95); fe – число степеней свободы дисперсии воспроизводимости (при одной серии параллельных опытов равная k -1). Выборочное значение квадратного корня дисперсии кодированного коэффициента регрессии определяется по формуле: e ~ S a~ j = C jj S e , где Se - квадратный корень из дисперсии воспроизводимости, определяемой по k параллельным опытам в центре плана эксперимента: å (y k S e2 = j =1 2 ' 0S - yc' k -1 ) º SSe , fe где SSe - сумма квадратов дисперсии воспроизводимости; fe - число степеней свободы дисперсии воспроизводимости. Как было показано выше, диагональные элементы корреляционной матрицы в ПФЭ при кодировании факторов одинаковы и равны 1/n, вследствие чего S S a~ j = e . n В результате условие не значимости кодированных коэффициентов регрессии принимает вид: a~ j n £ t bтабл ( fe ) . Se 19 ~ Так как корреляционная матрица C в этом случае является диагональной, то кодированные коэффициенты регрессии статистически независимы и при одновременной не значимости нескольких кодированных коэффициентов регрессии они (в отличие от процедуры обработки пассивного эксперимента) могут быть сразу, все вместе, исключены из кодированного уравнения регрессии. 4. Проверка адекватности уравнения регрессии (ПФЭ) Проводится так же, как и при проведении пассивного эксперимента, с использованием табличного значения критерия Фишера, выбранного при доверительной вероятности β (чаще всего равной 0,95) и числе степеней свободы остаточной дисперсии fR и дисперсии воспроизводимости fe . Условие адекватности проверяется с использованием неравенства: S R2 £ Fbтабл ( f R , fe ), 2 Se где остаточная дисперсия, характеризующая точность уравнения, определяется по формуле: å (yˆ n 2 I i -y ' i ) SS R . n- p fR При этом fR = n - p , где n – число экспериментов при различных значениях факторов; p - число значимых коэффициентов регрессии. К недостаткам ПФЭ относится резкое увеличение числа опытов при возрастании количества факторов больше, чем 5 ( при m = 5 n = 25 = 32 ). Для проведения регрессионного анализа при пренебрежении целым рядом несущественных взаимодействий факторов достаточно проводить меньшее число опытов. В этом случае можно реализовать часть ПФЭ, т.н. дробный факторный эксперимент (ДФЭ), который здесь не рассматривается. S = 2 R i= 1 = 5. ОЦКП и обработка его результатов Ортогональный центральный композиционный эксперимент (ОЦКП) относится II к экспериментам II – го порядка, так как описывающее его уравнение ŷ включает факторы в квадрате и поэтому может описывать поверхности функций отклика в окрестности их экстремальных значений. Для двух факторов ( x1 и x2 ) с учётом только двойного взаимодействия факторов соответствующая эмпирическая модель может быть записана: yˆ II = a0 + a1 x1 + a2 x2 + a12 x1 x2 + a11 x12 + a22 x22 . В соответствии с методикой ортогонального центрального композиционного плана эксперимента (ОЦКП) здесь, также как и для ПФЭ, осуществляется кодирование факторов по приведённой выше схеме, и для обеспечения ортогонального 20 II свойства матрицы планирования эксперимента в уравнение регрессии ŷ включается некоторая постоянная S. В результате уравнение регрессии при m = 2 принимает вид: ( ) ( ) yˆ = a~0 + a~1 z1 + a~2 z 2 + a~12 z1 z 2 + a~11 z12 - S + a~22 z 22 - S . Для определения большего числа кодированных коэффициентов, чем при обработке ПФЭ, и описания поверхности функции отклика вблизи её экстремума («почти стационарной области»), количество опытов в этом случае увеличивается. При этом опыты, проводимые при ПФЭ n = 2m , дополняются опытами в «звёздных» точках факторного пространства nα = 2m и опытами в центре плана с координатами z1 = 0 и z2 = 0 ( nc ). «Звёздные» точки в факторном пространстве располагаются на осях координат на расстоянии + α и – α от центра плана эксперимента; причём величина α называется «звёздным» плечом и её значения, так же как величина S, определяются из условия ортогональности матрицы планирования z для ОЦКП. Общее число опытов N в ортогональном центральном композиционном эксперименте определяется по формуле: N = n + nα + nc , или с учётом приведённых выше равенств: N = 2m + 2m + nc . Для случая двух факторов (m = 2): N = 8 + nc . Расположение опытных точек в факторном пространстве для случая двух факторов в приведённой ранее кодированной системе координат может быть представлено: z2 (- 1,+1) (+ 1,+1) +a (0,0 ) -a -a +a (- 1,-1) z1 (+ 1,-1) План проведения экспериментов в этом случае может быть представлен: 21 p z0 z1 z2 z1 z2 z12 - S z 22 - S y' 1 +1 -1 -1 +1 1- S 1- S y1' 2 +1 +1 -1 -1 1- S 1- S 3 +1 -1 +1 -1 1- S 1- S y2' y3' 4 +1 +1 +1 +1 1- S 1- S y4' 5 +1 -a 6 +1 +a 0 a -S -S 0 a2 - S - S 7 +1 -a 8 +1 9 +1 M +1 N +1 n 2m 2m nc 2 y5' y6' 2 -S a -S y7' +a - S a2 - S -S -S y8' y9' -S -S -S -S M y N' Матрица планирования z представляет собой часть плана проведения эксперимента без горизонтальных и вертикальных заголовков таблицы и вектора наблюдения y ' (правого столбца). 6. Определение величины «звёздного плеча» α и S из условия ортогональности матрицы планирования z . z Матрица планирования ( N ´6 ) была бы ортогональной, если бы выполнялись следующие равенства: ìï z0T (z j - S ) = 0 í ïî j = 1,2 и (z )( T ) - S z22 - S = 0. Раскрывая первое равенство, можно получить: 2 1 z0 (z j - S ) = å z z - å zi 0 S = n + 2a 2 - NS = 0 N n 2 i 0 ij i =1 i =1 j = 1,2 . Откуда: n + 2a 2 S= N ( A). Раскрывая второе равенство, получаем: (z 2 1 -S ( ) (z T 2 2 n - n + 2a 2 Откуда: ) ( ) ( ) )S - S (n + 2a ) + NS = n - 2 NS T T - S = z12 z22 - z12 S - S T z22 + S T S = 2 2 2 + NS 2 = n - NS 2 = 0 22 S= n N (B ). Последнее выражение используется для определения S. Приравнивая правые части двух выражений для S, можно найти формулу для определения α : n + 2a 2 = N a2 = 1 2 ( ) Nn - n = n N ö næ N ç ÷. 1 ÷ 2 çè n ø В результате звёздное плечо α можно определить по формуле: a= ö n æç N ÷. 1 ÷ 2 çè n ø 4. Определение кодированных коэффициентов регрессии (ОЦКП) В соответствии с методом наименьших квадратов эти коэффициенты определяются по матричной формуле: ~ T a~ = C z y ' , где ~ æ T ö C =ç z z÷ è ø -1 Из-за свойства ортогональности матрицы планирования z необходимо определить только диагональные элементы информационной матрицы: T ~ I = z z, а затем диагональные элементы корреляционной матрицы: ~ ~ -1 C=I . 8. Определение диагональных элементов информационной и корреляционной матриц. Обобщая уравнение регрессии на случай m факторов и учитывая только все двойные взаимодействия факторов, число которых определяется по формуле: 23 Cm2 = m(m - 1) , 2! общее число коэффициентов уравнения регрессии для m факторов равно: m(m - 1) p = 1+ m + + m, 2! ~ диагональные элементы информационной матрицы I i00 = N - число таких элементов равно 1; определяются: ijj = n + 2α2 ( j = 1, …m ); m(m - 1) 2! . iju = n ( u > j ) - число таких элементов равно: Для определения при квадратах факторов ijj можно записать: ( ) i jj = n(1 - S ) + 2 a 2 - S + (N - n - 2 )S 2 = 2 2 = n - 2nS + nS 2 + 2a 4 - 4a 2 S + 2 S 2 + NS 2 - nS 2 - 2S 2 = ( ) 2 n4 NS = 2a 4 + n - 2 S n + 2a 2 S + NS 2 = 2a 4 + 1 -2 4 3 14243 = 2a 4 из равенства ( B )= 0 из равенства ( A )= NS Количество таких диагональных элементов – m. Диагональная матрица I имеет размер: (m + 1)(m + 2) , m(m - 1) p = 1+ m + +m= 2 2 что соответствует числу определяемых параметров p. ~ ~ В результате диагональная корреляционная матрица C = I -1 размером р х р для m факторов и с учётом их двойных взаимодействий имеет вид: é(N )-1 ê ê ê ê ê ê C =ê p´ p ê ê ê ê ê ê êë (n + 2a )-1 L (n + 2a )-1 n -1 L n -1 (2a ) 4 -1 L 1 ù ú ú m ú ú ú ú ú m(m - 1) ú 2 ú ú ú ú m ú 4 -1 ú 2a û ( ) Элементы корреляционной матрицы определяются по матричной формуле МНК: 24 ~ a~ = C z T y ' . Кодированные коэффициенты регрессии определяются: N a~0 = a~ ju = i =1 ij iu ' i a~ j = ; N å (z z )y N N å yi' n å (z n a~ jj = i =1 i =1 ij yi' n + 2a 2 ( j = 1,...m); m(m - 1) 2 ) (число коэффициентов u> j i =1 åz 2 ij ) - S yi' 2a 4 ( j = 1,...m). Для пересчёта этих коэффициентов регрессии в натуральные значения необходимо вместо кодированных факторов z подставить их натуральные величины xj в соответствии с приведённой схемой кодирования. 9. Определение значимости кодированных коэффициентов регрессии (ОЦКП) В отличие от ПФЭ значимость коэффициентов регрессии определяется по разным формулам для различных коэффициентов, так как диагональные элементы кор~ реляционной матрицы C отличаются друг от друга. С учётом общей формулы для определения не значимости коэффициентов регрессии a~ j £ t bтабл (Se ) , ~ C jj Se не значимость каждого вида коэффициентов регрессии определяется: a~0 Se a~ j Se a~ ju Se n £ t bтабл ( Se ) a~ jj Se N £ t bтабл (Se ) ; n + 2a 2 £ t bтабл ( S e ) ( j = 1,...m) m(m - 1) 2 (число коэффициентов ) 2a 4 £ t bтабл ( Se ) ( j = 1,...m). 25 10. Проверка адекватности уравнения регрессии (ОЦКП) Осуществляется с использованием критерия Фишера – так же, как и в случае с ПФЭ. 11. Определение экстремума функции отклика Уравнение регрессии с m факторами вида: m m -1 m m ( yˆ II = a~0 + å a~ j z j + åå a~ ju z j zu + å a~ jj z 2j - S j =1 j =1 u = 2 ) j =1 u> j может применяться для определения экстремума функции отклика с использованием необходимого условия экстремума функции многих переменных: ¶ yˆ II = 0 ¶ z1 L L L L ¶ yˆ II = 0. ¶zm Полученная система линейных уравнений (СЛАУ) позволяет расчётным путём оп- II ределить z j ( j = 1,...m) и после подстановки их величин в исходное уравнение ŷ получить максимальное или минимальное значение функции отклика. opt Лекция 3 Планы дробного факторного эксперимента (планы ДФЭ) При многофакторном эксперименте, особенно когда число факторов больше шести (n > 6), число опытов планов ПФЭ 2n (N = 2n) становится чрезмерным. Если нам не требуется определение всех коэффициентов неполного квадратичного полинома, то переходят к дробному факторному эксперименту (ДФЭ) – части полного факторного эксперимента. Так, например, если требуется определить лишь коэффициенты при самих факторах , n , в этом случае трето план ПФЭ 2 дает избыточную информацию. Так при буется определить коэффициентов, тогда как по плану ПФЭ необходимо 6 N провести = 2 =64 опыта. Хотя эта избыточная информация не является бесполезной, она позволяет более точно определить коэффициенты, но все же часто используют планы ДФЭ 2n-k , где k – показатель дробности плана ПФЭ. При k = 1 число опытов в плане ДФЭ в два раза меньше, чем в плане ПФЭ, поэтому такие планы называют полуреплика плана ПФЭ. Так при k=1 для плана ДФЭ 26-1 N =26-1 = 32, при k=2 для плана ДФЭ 26-2 N =26-2 = 16 и такой план называют четвертьрепликой, при k=3 для плана ДФЭ 26-3 N =26-3 = 8. 26 При выборе дробности плана k необходимо учитывать, что число опытов должно быть больше числа членов уравнения. В рассматриваемом случае величина k должна быть такой, что бы удовлетворялось условие . План ДФЭ строится, как и для плана ПФЭ, но с меньшим числом факторов. Оставшиеся факторы варьируются не произвольно, а так чтобы сохранялась ортогональность плана. Это обеспечивается, если оставшиеся факторы варьируются по выбранному генерирующему соотношению, например как произведение каких-либо факторов из первой группы. Но это приводит к тому, что в матрице Х будут существовать одинаковые столбцы. Следовательно, мы не сможем найти в чистом виде все коэффициенты неполного квадратичного полинома, а лишь определим совместную величину коэффициентов для одинаковых столбцов. Рассмотрим построение плана ДФЭ 23-1 . Здесь n = 3, к =1, N=23-1=4. Первые два фактора варьируем как и ранее для плана ПФЭ 22, а для третьего фактора выбираем генерирующее соотношение в виде . Для неполного квадратичного полинома количество столбцов плана составляет восемь. План является ортогональным, но в нем оказались четыре пары одинаковых столбцов. Поэтому можно определить только четыре коэффициента, отражающие совместные влияния двух одинаковых столбцов . Суммарные значения коэффициентов ; ; определяются аналогично. Это следствие того, что мы пытаемся определить полное количество коэффициентов – 8 по недостаточному числу опытов - 4. Однако, если заранее известно, что некоторые из членов уравнения равны нулю (пренебрежимо малы) или имеется априорная информация о величинах некоторых коэффициентов, то полученные коэффициенты могут быть вычленены. Так если , то . 27 Если можно допустить, что коэффициенты из их смешанной оценки сопоставимы, то для рассмотренного плана . Графическое изображение планов ПФЭ 23 и ДФЭ 23-1 в факторном пространстве (для трех факторов - трехмерное пространство) представлено на рис. 1. План ПФЭ 23 представлен кубом с восемью узлами (точками плана), а возможные планы ДФЭ 23-1 – проекциями этого куба на три плоскости. То есть из восьми узлов выбираются четыре (рис. 1, а). Из куба можно также выбрать четыре точки из восьми, не лежащие в одной плоскости, и сформировать план ДФЭ 23-1 (рис. 1, б). Насыщенные планы первого порядка. Насыщенным планом первого порядка – называется план, содержащий n+1 точку (опыт). Например, при n = 4, N=n + 1 = 5. То есть полином формируется в виде . Таким образом, насыщенный план – это предельно минимальный случай плана ДФЭ. Такие планы называются симплекс-планами. Для симплекс-плана при n = 1 N = 2 его геометрическое изображение представлено на рис. 2, а; при n=2, N=3 – на рис. 2, б; при n=3, N=4 – на рис. 2, в. Симплекс-планы обычно используются на стадии предварительного исследования. Рис. 2. Симплекс-план для n=1, N=2 (а); n=2, N=3 (б); n=3, N=4 (в) Симплекс-план не всегда является ортогональным. Симплекс-план называется правильным, если расстояние между двумя любыми точками плана одинаковое. Симплекс-план называется центрированным, если , для i=1, 2, …, n . Применимость планов ПФЭ и пути повышения точности полиномов По каким же признакам можно судить о допустимости использования неполного квадратичного полинома, построенного на основе планов ПФЭ 2n? 28 Такие полиномы дают поверхность отклика, которая проходит точно через все экспериментальные точки, по которым определяются коэффициенты. Так как точки планов ПФЭ располагаются на границах диапазонов варьирования факторов, то это означает, что поверхность отклика проходит через граничные точки. В любом сечении поверхности отклика, полученной по такому полиному, плоскостью при фиксированных всех факторах кроме одного и параллельной оси Y получается след в виде прямой линии. Возможны случаи, когда реальная поверхность отклика определяется полиномами второго и выше порядков В этом случае поверхность плана ПФЭ, совпадая с реальной поверхностью в граничных точках, может отличаться в других точках факторного пространства, например в центральной точке плана, т.е. . Поэтому одним из признаков неудовлетворительной аппроксимации полиномами по плану ПФЭ является расхождение результатов функции отклика с реальной функцией в центральной точке плана. Однако при многофакторном эксперименте возможны случаи, когда в реальности функция отклика зависит, в том числе, от квадратов факторов, у которых коэффициенты имеют разные знаки, например, для “седловидной” поверхности. При этом, несмотря на то, что эта поверхность явно нелинейная, результат опыта в центральной точке может оказаться достаточно близким к полученному результату по неполному квадратичному полиному плана ПФЭ. Однако расхождения будут возникать во всех других точках плана эксперимента. Поэтому нецелесообразность использования плана ПФЭ определяется нелинейностью каких-либо сечений поверхности отклика. Косвенным признаком может служить расхождение и в центральной точке плана. Если не удается получить полином по плану ПФЭ, хорошо аппроксимирующей реальную поверхность, то какие пути можно предложить для повышения точности полиномов? · Уменьшение диапазона варьирования факторов или его разбиение на поддиапазоны, для каждого из которых строится свой план ПФЭ и определяется свой полином. Путь достаточно трудоемок, но погрешность семейства планов ПФЭ снижается. · Выделение фактора, порождающий нелинейность, и построение для оставшихся n–1 факторов k планов ПФЭ, в каждом из которых выделенный фактор зафиксирован при некотором значении. На основе полученных k полиномов можно попытаться сформировать общий полином, коэффициенты которого являются функциями выделенного фактора. Этот путь также достаточно трудоемок. · Переход к плану ПФЭ с большим числом уровней варьирования факторов, например к планам с варьированием факторов на трех уровнях - планам ПФЭ 3n (рис. 3). В этом случае происходит резкое увеличение количества точек по сравнению с планом ПФЭ 2n . Так при n = 2 для ПФЭ 2n N=4, для ПФЭ 3n N=9; при n = 3 для ПФЭ 2n N=8, для ПФЭ 3n N=27; при n = 4 для ПФЭ 2n N=16, для ПФЭ 3n N=81 и т.д. 29 Рис. 3. Планы ПФЭ 32 (а) и ПФЭ 33 (б) Достраивание планов ПФЭ 2n до планов более высокого порядка (чаще всего второго) и построение полных квадратичных полиномов (с наличием квадратов факторов). · Преобразование метрики матричного пространства, то есть переход к новым факторам функционально связанным с прежними факторами, но не порождающими нелинейности. · Рис. 4. Графическое изображение планов ПФЭ 23 и ДФЭ 23-1 в факторном пространстве Планы ДФЭ, как и планы ПФЭ, являются рототабельными. Планы ДФЭ могут быть как насыщенными так и ненасыщенными. Достоинство планов ДФЭ заключается и в том, что если построенный на его основе неполный полином не удовлетворяет требованиям по точности, то план ДФЭ легко достраиваются до плана ПФЭ, без потери информации прежних опытах, с формированием более точного полинома. Пример построения плана ДФЭ. Построить план ДФЭ 24-1 и определить полином . Число факторов – 4. Нужно найти 8 коэффициентов полинома. Выбираем 8 из 16 опытов плана ПФЭ 24 таким образом, чтобы были определены независимые коэффициенты при самих факторах, смешанные коэффициенты при парных сочетаниях 30 факторов и в пренебрежении тройными и четверным сочетаниями факторов и при этом сохранялась ортогональность плана. Такой выбор позволяет сформировать план ДФЭ 24-1 как и план ПФЭ 23 , но с х4= х1х2х3 . План ДФЭ 24-1 представляется в виде Значения коэффициентов полинома составляют: . ; ; 31 ; Если принять, что , , , то полином имеет вид . Значения полинома в точках плана приведены в последнем столбце плана ДФЭ 2 . В нашем случае точность его достаточно высокая. 4-1 Лекция 4 Дробный факторный эксперимент При большом числе факторов (k>3) проведение полного факторного эксперимента связано с большим числом, экспериментов, значительно превосходящим число коэффициентов линейной модели. Если при получении модели можно ограничиться, линейным приближением, т. е. получить адекватную модель в виде полинома y=b0 + b1x1 + b2x2 +...+ bkxk, то число экспериментов можно резко сократить в результате использования дробного факторного эксперимента. Так, например, в полном факторном эксперименте типа 22 при линейном приближении коэффициент регрессии b12 можно принять равным нулю, а столбец x1x2 матрицы (табл. 4) использовать для третьего фактора x3. Таблица 4 Матрица планирования Номер экспер. 1 2 3 4 x0 x1 x2 x3 (x1x2) y + + + + + + - + + - + + y1 y2 y3 y4 В этом случае линейная модель будет определяться уравнением y=b0 + b1x1 + b2x2 + b3x3. Для определения коэффициентов этого уравнения достаточно провести четыре эксперимента вместо восьми в полном факторном эксперименте типа 32 23. План эксперимента, предусматривающий реализацию половины экспериментов полного факторного эксперимента, называют полурепликой. При увеличении числа факторов (k>3) возможно применение реплик большей дробности. Дробной репликой называют план эксперимента, являющийся частью плана полного факторного эксперимента. Дробные реплики обозначают зависимостю 2k-p, где p - число линейных эффектов, приравненных к эффектам взаимодействия. При p = 1 получают полуреплику; при p = 2 получают 1/4 - реплику; при p = 3 получают 1/8 - реплику и т. д. по степеням двойки. Так, например, если в полном факторном эксперименте 23 (табл. 5) один из эффектов взаимодействия (x1x2, x1x3, x2x3, x1x2x3) заменим четвертым фактором x4, то получим полуреплику 24-1 от полного факторного эксперимента 24. Если два эффекта взаимодействия заменить факторами x4 и x5, то получим ¼-реплику 25-2 от полного факторного эксперимента 25. Таблица 5 Матрица полного факторного эксперимента типа 2 3 Номер экспер. x0 x1 x2 x3 x1x2 x1x3 x2x3 x1x2x3 y 1 2 3 4 5 6 7 8 + + + + + + + + + + + + + + + + + + + + - + + + + + + + + - + + + + - + + + + - y1 y2 y3 y4 y5 y6 y7 y8 Можно получать 1/8 -реплику от полного факторного эксперимента 26, заменив три эффекта взаимодействия факторами x4, x5 и x6. Если заменить четыре эффекта взаимодействия факторами x4, x5, x6 и x7, то получим,1/16 -реплику 27-4 от полного факторного эксперимента 27. Реплики, которые используют для сокращения числа экспериментов в 2m раз, где m=1, 2, 3 ..., называют регулярными. В связи с тем, что в дробных репликах часть взаимодействий заменена новыми факторами, найденные коэффициенты уравнения регрессии будут являться совместными оценками линейных эффектов и эффектов взаимодействия. Так, например, если в матрице (табл. 4) вычислим элементы столбцов для произведений x1x3 и x2x3, то увидим, что элементы столбца x1x2 совпадают с элементами столбца x2, а элементы столбца x2x3 - с элементами столбца x1. Следовательно, коэффициенты b1, b2, b3 будут оценками совместных эффектов, а именно: b1 → β1 + β23; b2→ β2 + β13; b3→ β3 + β12. Коэффициент b1 является оценкой влияния фактора x1 и парного взаимодействия x2x3 на функцию отклика. Влияние фактора x1 в этом случае характеризуется величиной β1, а влияние взаимодействия - величиной β23. Оценки, в которых невозможно разделить линейный эффект и эффект взаимодействия, называют сме33 шанными. Линейные эффекты рекомендуется смешивать, прежде всего, с теми взаимодействиями, которые согласно априорной информации незначимы. Число несмешанных линейных эффектов в дробной реплике называют ее разрешающей способностью. Часто приходится решать задачи, в которых заранее можно полагать, что эффекты взаимодействия, хотя и малы по сравнению с линейными, но все же не равны нулю. В таких случаях необходимо заранее определить, какие коэффициенты являются смешанными оценками. Тогда в зависимости от условий поставленной задачи, подбирается такая дробная реплика, с помощью которой можно извлечь максимальную информацию из эксперимента. Прямая оценка разрешающей способности дробной реплики затруднена. Поэтому дробные реплики задают с помощью генерирующих соотношений. Генерирующим называют соотношение, которое показывает, какое из взаимодействий принято незначимым и заменено новым фактором. План типа 23-1 может быть представлен двумя полурепликами (табл. 6), которые задаются одним из следующих генерирующих соотношений: x3 = x1x2; x3 = - x1x2. Генерирующие соотношения умножим на новую независимую переменную x3: x32 = x1 x2 x3 : x32 = - x1 x2 x3 . Таблица 6 3-1 Две полуреплики плана типа 2 x3 = x1x2 Номер экспер. Номер экспер. 1 2 3 4 x1 x2 x3 + + + + - + + - 1 2 3 4 x3 = - x1x2 x1 x2 x3 + + + + - + + Поскольку всегда xi = 1 , получим следующие соотношения: 1 = x1x2x3; 1 = - x1x2x3. (3) В результате умножения генерирующего соотношения на новую переменную получают так называемый определяющий контраст. Для указанных выше полуреплик определяющими контрастами будут зависимости (3). Зная определяющий контраст, можно найти соотношения, задающие совместные оценки. Для этого необходимо умножить независимые переменные x1, x2 и x3 на определяющий контраст. Умножая определяющие контрасты (3) на x1, получим соотношения 2 x1 ×1 = x12 x2 x3 ; x1 ×1 = - x12 x2 x3 ; так как x1 = 1 , то 2 ния: x1 = x2x3; x1 = - x2x3. Умножая определяющие контрасты на x2 и x3, получаем следующие соотноше34 x2 = x1x3; x2 = - x1x3; x3 = x1x2; x3 = - x1x2. Это означает, что коэффициенты регрессии будут оценками b1 → β1 + β23; b1→ β1 – β23; b2 → β2 + β13; b2→ β2 – β13; b3 → β3 + β12; b3→ β3 – β12. 4-1 Полуреплика 2 может быть задана генерирующим соотношением x4= x1x2x3. Матрица планирования этой полуреплики представлена табл. 7. Определяющим контрастом полуреплики является соотношение 1= x1x2x3x4. Совместные оценки будут определяться следующим образом: x1= x2x3x4 b1 → β1 + β234; x2= x1x3x4 b2 → β2 + β134; x3= x1x2x4 b3 → β3 + β124; x4= x1x2x3 b4 → β4 + β123; x1x2 = x2x4 b12 → β12 + β34; x1x3 = x3x4 b13 → β13 + β24; x1x4 = x2x3 b14 → β14 + β23. 4-1 Полуреплика 2 может быть также задана генерирующим соотношением x4= x1x2. Матрица планирования этой полуреплики представлена табл. 8. Таблица 7 определяющим Полуреплика 24-1 с трастом 1= x1x2x3x4. Номер экспер. 1 2 3 4 5 6 7 8 x0 x1 x2 x3 x4 y + + + + + + + + + + + + + + + + + + + + - + + + + - y1 y2 y3 y4 y5 y6 y7 y8 Таблица.8 Полуреплика 24-1 с определяющим контрастом 1= x1x2x4 Номер экспер. 1 2 3 4 5 6 7 8 x0 x1 x2 x3 x4 y + + + + + + + + + + + + + + + + + + + + - + + + + y1 y2 y3 y4 y5 y6 y7 y8 кон- Определяющим контрастом полуреплики является соотношение 1= x1x2x4. Совместные оценки в этом случае будут определяться следующим образом: x1= x2x4 b1 → β1 + β24; x2= x1x4 b2 → β2 + β14; x3= x1x2 x3x4 b3 → β3 + β1234; x4= x1x2 b4 → β4 + β12; x1x3 = x2x3x4 b13 → β13 + β234; x2x3 = x1x3x4 b23 → β23 + β134; x3x4 = x1x2x3 b34 → β34 + β123. 35 В практических задачах тройные и более высокого порядка взаимодействия значительно чаще, чем двойные, бывают равны нулю и ими обычно можно пренебречь. Полуреплика 24-1, заданная генерирующим соотношением x4=x1x2x3, позволяет получить раздельные оценки четырех линейных эффектов и три совместные оценки парных взаимодействий. В этом случае раздельными оценками будут b1, b2, b3 и b4, так как тройными взаимодействиями β234, β134, β124 и β123 вследствие их незначимости можно пренебречь. В полуреплике, заданной генерирующим соотношением x4=x1x2, три линейных эффекта, а именно b1, b2, b4 - оказались смешанными с парными взаимодействиями. Разрешающая способность полуреплики, заданной генерирующим соотношением x4=x1x2x3, получилась значительно выше, чем у полуреплики, заданной генерирующим соотношением x4=x1x2. Следовательно, разрешающая способность полуреплики зависит от генерирующего соотношения, которым она задана. Для оценки разрешающей способности реплик (большой дробности (¼, ⅛ и т. д.) используют обобщающие определяющие контрасты. ¼-реплика 25-2 может быть задана следующими генерирующими соотношениями: x4=x1x2x3, x5=x2x3. Матрица планирования этой реплики представлена табл. 9. Таблица 9 Матрица планирования 25-2 Номер экспер. 1 2 3 4 5 6 7 8 x0 x1 x2 x3 x4 x5 y + + + + + + + + + + + + - y1 y2 y3 y4 y5 y6 y7 y8 + + + + - + + + + - + + + + + + + + Таблица 10 Матрица планирования 27-4 Номер экспер. 1 2 3 4 5 6 7 8 x0 x1 x2 x3 x4 x5 x6 x7 y + + + + + + + + + + + + - y1 y2 y3 y4 y5 y6 y7 y8 + + + + - + + + + + + + + - + + + + + + + + + + + + - Определяющими контрастами реплики являются соотношения 1= x1x2x3x4; 1= x2x3x5. Перемножив определяющие контрасты, получим третье соотношение 1= x1x4x5. Полная характеристика разрешающей способности рассматриваемой реплики будет определяться обобщающим определяющим контрастом, имеющим вид 1= x1x2x3x4= x2x3x5= x1x4x5. Схему смешивания оценок находим последовательным умножением обобщающего определяющего контраста на x1, x2, x3 и т. д. x1= x2x3x4= x1x2x3x5= x4x5 b1 → β1 + β234 + β1235 + β45; x2= x1x3x4= x3x5=x1x2x4x5 b2 → β2 + β134 + β35 + β1245; x3= x1x2x4= x2x5=x1x3x4x5 b3 → β3 + β124 + β25 + β1345; 36 x4= x1x2x3= x2x3x4x5= x1x5 b4 → β4 + β123 + β2345 + β15; x5= x1x2x3x4x5=x2x3= x1x4 b5 → β5 + β12345 + β23 + β14; x1x2= x3x4= x1x3x5=x2x4x5 b12 → β12 + β34 + β135 + β245; x1x3= x2x4= x1x2x5=x3x4x5 b13 → β13 + β24 + β125 + β345. Для 1/16 реплики генерирующими соотношениями x4= x1x2x3; x5= x1x2; x6= x1x3; x7= x2x3 матрица планирования представлена табл. 10. Определяющими контрастами этой реплики будут соотношения 1) 1= x1x2x3x4; 2) 1= x1x2x5; 3) 1= x1x3x6; 4) 1= x2x3x7. Если попарно перемножить определяющие контрасты 1х2; 1х3; 1х4; 2х3; 2х4; 3х4, то получим 1= x3x4x5; 1= x2x4x6; 1= x1x4x7; 1= x2x3x5x6; 1= x1x3x5x7; 1= x1x2x6x7. Произведения определяющих контрастов по три: 1x2x3; 1х2х4; 2х3х4; 1х3х4 будут равны 1= x1x4x5x6; 1= x2x4x5x7; 1= x5x6x7; 1= x3x4x6x7. Умножая определяющие контрасты по четыре, получим 1= x1x2x3x4x5x6x7. Чтобы полностью характеризовать разрешающую способность данной реплики, запишем обобщающий определяющий контраст 1= x1x2x3x4= x1x2x5=x1x3x6=x2x3x7=x3x4x5=x2x4x6= x1x4x7= x2x3x5x6= x1x3x5x7= x1x2x6x7= x1x4x5x6= x2x4x5x7= x5x6x7=x3x4x6x7= x1x2x3x4x5x6x7. Если эффектами взаимодействия, начиная с тройных, можно пренебречь, то коэффициенты будут оценками: b1 → β1 + β25 + β36 + β47; b2 → β2 + β15 + β37 + β46; b3 → β3 + β16 + β27 + β45; b4 → β4 + β35 + β26 + β17; b5 → β5 + β12 + β34 + β67; b6 → β6 + β13 + β24 + β57; b7 → β7 + β23 + β14 + β56. Таким образом, получаем весьма сложную систему смешивания. Все линейные эффекты оказались смешанными с несколькими парными взаимодействиями, поэтому разрешающая способность этой дробной реплики очень низкая. Пользоваться такой репликой можно лишь в том случае, если все парные взаимодействия близки к нулю. Выбор дробной реплики зависит от конкретной задачи. Для получения линейной модели рекомендуют выбирать дробные реплики с возможно большей разрешающей способностью, т. е. реплики, у которых линейные эффекты смешаны с эффектами взаимодействия близкими к нулю. При выборе дробной реплики важно учитывать насыщенность плана, т. е. соотношение между числом опытов и числом коэффициентов, определяемых по результатам этих экспериментов. Дробная реплика, полученная заменой всех эффектов взаимодействия новыми факторами, называется насыщенной. Применение насыщенных планов требует минимального числа экспериментов. Число экспериментов в матрице насыщенной дробной репли37 ки равно числу коэффициентов линейной модели. Гипотезу адекватности модели в этом случае проверить невозможно, так как число степеней свободы равно нулю. Например, 1/16-реплика от полного факторного эксперимента 27 (табл. 10) является насыщенной, так как линейная модель не содержит коэффициентов, которые необходимо определить по результатам восьми экспериментов. При этом не остается степеней свободы дли проверки адекватности модели. Дробные реплики широко применяют при получении линейных моделей. Эффективность применения дробных реплик зависит от удачного выбора системы смешивания линейных эффектов с эффектами взаимодействия. При построении дробных реплик используют следующее правило: новый фактор, введенный в планирование, нужно поместить в столбец матрицы, принадлежащий взаимодействию, которым можно пренебречь. Лекция 5 Проведение эксперимента и обработка его результатов После выбора плана эксперимента, основных уровней и интервалов варьирования факторов переходят к эксперименту. Каждая строка матрицы – это условия эксперимента. Для исключения систематических ошибок рекомендуется эксперименты, предусмотренные матрицей, проводить в случайной последовательности. Порядок проведения следует выбирать по таблице случайных чисел (табл. 11). Например, если требуется провести восемь экспериментов, то из случайного места таблицы последовательно выписывают числа, лежащие в интервале от 1 до 8, при этом не учитываются уже выписанные и числа больше восьми. Так, например, начиная с числа 87 (1-я строка табл. 11), получаем следующую последовательность реализации экспериментов: Номер опыта в матрице 1 Порядок реализации экспериментов 7 2 2 3 8 4 3 5 6 7 8 1 4 5 6 Таблица 11 87 07 03 13 00 24 61 27 96 Фрагмент таблицы случайных чисел 63 88 23 62 51 07 69 59 02 76 85 37 84 37 47 32 25 21 33 48 84 37 37 29 38 37 89 01 59 47 64 04 99 59 96 20 83 48 94 44 08 67 79 41 61 07 78 61 89 42 58 88 22 16 12 90 62 41 11 59 85 18 42 84 05 99 85 75 67 80 05 57 53 99 25 13 63 89 15 76 30 41 13 61 05 49 08 25 87 15 24 29 71 14 82 09 31 60 40 88 70 98 34 69 33 11 09 76 21 53 57 44 69 88 00 04 31 41 57 61 45 83 65 21 99 92 35 88 58 24 46 80 99 36 22 23 48 82 38 78 06 Для компенсации влияния случайных погрешностей каждый эксперимент рекомендуется повторить n раз. Эксперименты, повторенные несколько раз при од38 них и тех же значениях факторов, называют параллельными. Под дублированием понимают постановку параллельных экспериментов. Обычно число n параллельных экспериментов принимают равным 2-3, иногда – 4-5. При проведении исследований приходится иметь дело с тремя вариантами дублирования экспериментов: 1. с равномерным дублированием экспериментов; 2. с неравномерным дублированием экспериментов; 3. без дублирования экспериментов. При равномерном дублировании все строки матрицы планирования имеют одинаковые числа параллельных экспериментов. В случае неравномерного дублирования числа параллельных экспериментов неодинаковы. При отсутствии дублирования параллельные эксперименты не проводятся. Наиболее предпочтительным из трех вариантов дублирования является первый. При этом варианте эксперимент отличается повышенной точностью, а математическая обработка экспериментальных данных - простотой. Характер дублирования влияет на содержание математической обработки результатов наблюдений. Рассмотрим методику обработки результатов эксперимента для каждого из трех вариантов дублирования. Обработка результатов эксперимента при равномерном дублировании. Для каждой строки матрицы планирования по результатам n параллельных экспериментов находят y j среднее арифметическое значение параметра оптимизации: 1 n y j = å y ju , n u =1 где u - номер параллельного эксперимента; yju - значение параметра оптимизации в u-м параллельном эксперименте j-й строки матрицы. С целью оценки отклонений параметра оптимизации от его среднего значения s2 для каждой строки матрицы планирования вычисляют дисперсию j эксперимента по данным n параллельных экспериментов. Статистической дисперсией называют среднее значение квадрата отклонений случайной величины от ее среднего значения: 1 n ( y ju - y j ) 2 . s = å n - 1 u= 1 2 j (4) Ошибка sj эксперимента определяется как корень квадратный из дисперсии sj = + 1 n ( y ju - y j ) 2 . å n - 1 u =1 В этом случае ошибка при большом рассеянии будет значительной. Рассеяние результатов эксперимента определяется влиянием неуправляемых факторов, погрешностями измерений и другими причинами. Большое рассеяние изучаемой величины может произойти из-за наличия в эксперименте сомнительных результатов. Для проверки сомнительных, т. е. резко выделяющихся результатов, используют специальные критерии; одним из таких критериев является отношение U (ГОСТ 11.002-73). Чтобы оценить принадлежность резко выделяющихся результатов yj max 39 или yj min к данной нормальной совокупности и принять решение об исключении или оставлении их в составе выборки, находят отношение y - yj y - y j min U max = j max , U min = j sj s j или где yj max - наибольшее значение параметра оптимизации среди его значений, полученных в n параллельных экспериментах j-й строки матрицы планирования; yj min наименьшее значение параметра оптимизации среди его значений, полученных в n параллельных экспериментах j-й строки матрицы планирования. Результат сравнивают с величиной b, взятой из ГОСТ 11.002-73 (табл. 1) для числа n параллельных экспериментов и принятого уровня значимости a. Число n параллельных экспериментов и объем выборки n в рассматриваемом случае понятия равноценные. Если Umax³b, то сомнительный результат может быть исключен, в противном случае его считают нормальным и не исключают. Аналогично производится оценка результата yj min: если Umin³b, то сомнительный результат признают анормальным; при Umin s22 . Если наблюдаемое значение Fp-критерия меньше табличного Fт (табл. 12) для соответствующих чисел степеней свободы и принятого уровня значимости, то дисперсии однородны. Однородность ряда дисперсий проверяют по критерию Кохрена или по критерию Бартлета. При равномерном дублировании экспериментов однородность ряда дисперсий проверяют с помощью G-критерия Кохрена, 40 представляющего собой отношение максимальной дисперсии к сумме всех дисперсий: 2 N smax 2 Gp = 2 2 = s s2. / max å j 2 s1 + s2 + ... + s N j=1 Дисперсии однородны, если расчетное значение Gp-критерия не превышает табличного значения Gт-критерия. В табл. 13 - N показывает число сравниваемых дисперсий, а n - число параллельных опытов. Если Gp> Gт, то дисперсии неоднородны, а это указывает на то, что исследуемая величина y не подчиняется нормальному закону. В этом случае нужно попытаться заменить y случайной величиной q = f(y), 2 достаточно близко следующей нормальному закону. Если дисперсии s j эксперимен2 тов однородны, то дисперсию s y воспроизводимости вычисляют по зависимости 1 s = N 2 y N ås , 2 j (5) где N - число экспериментов или число строк матрицы планирования. j =1 Таблица 12 Значения F-критерия Фишера при 5 % - ном уровне значимости Число степеней свободы для меньшей дисперсии 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 Значения критерия при числе степеней свободы для большей дисперсии 2 3 4 5 6 12 24 » 164,4 199,5 215,7 224,6 230,2 234,0 224,9 249,0 18,5 19,2 19,2 19,3 19,3 19,3 19,4 19,4 10,1 9,6 9,3 9,1 9,0 8,9 8,7 8,6 7,7 6,9 6,6 6,4 6,3 6,2 5,9 5,8 6,6 5,8 5,4 5,2 5,1 5,0 4,7 4,5 6,0 5,1 4,8 4,5 4,4 4,3 4,0 3,8 5,5 4,7 4,4 4,1 4,0 3,9 3,6 3,4 5,3 4,5 4,1 3,8 3,7 3,6 3,3 3,1 5,1 4,3 3,9 3,6 3,5 3,4 3,1 2,9 5,0 4,1 3,7 3,5 3,3 3,2 2,9 2,7 4,8 4,0 3,6 3,4 3,2 3,1 2,8 2,6 4,8 3,9 3,5 3,3 3,1 3,0 2,7 2,5 4,7 3,8 3,4 3,2 3,0 2,9 2,6 2,4 4,6 3,7 3,3 3,1 3,0 2,9 2,5 2,3 4,5 3,7 3,3 3,1 2,9 2,8 2,5 2,3 4,5 3,6 3,2 3,0 2,9 2,7 2,4 2,2 4,5 3,6 3,2 3,0 2,8 2,7 2,4 2,2 4,4 3,6 3,2 2,9 2,8 2,7 2,3 2,1 4,4 3,5 3,1 2,9 2,7 2,6 2,3 2,1 4,4 3,5 3,1 2,9 2,7 2,6 2,3 2,1 254,3 19,5 8,5 5,6 4,4 3,7 3,2 2,9 2,7 2,5 2,4 2,3 2,2 2,1 2,1 2,0 2,0 1,9 1,9 1,8 41 22 24 26 28 30 40 60 120 ¥ 4,3 4,3 4,2 4,2 4,2 4,1 4,0 3,9 3,8 3,4 3,4 3,4 3,3 3,3 3,2 3,2 3,1 3,0 3,1 3,0 3,0 3,0 2,9 2,9 2,8 2,7 2,6 2,8 2,8 2,7 2,7 2,7 2,6 2,5 2,5 2,4 2,7 2,6 2,6 2,6 2,5 2,5 2,4 2,3 2,2 2,6 2,5 2,5 2,4 2,4 2,3 2,3 2,2 2,1 2,2 2,2 2,2 2,1 2,1 2,0 1,9 1,8 1,8 2,0 2,0 2,0 1,9 1,9 1,8 1,7 1,6 1,5 1,8 1,7 1,7 1,7 1,6 1,5 1,4 1,3 1,0 Таблица 13 Значения G-критерия при 5 % - ном уровне значимости n-1 N 1 2 3 4 5 6 7 8 9 4 0,9065 0,7679 0,6841 0,6287 0,5895 0,5598 0,5365 0,5175 0,5017 6 0,7808 0,6161 0,5321 0,4803 0,4447 0,4184 0,3980 0,3817 0,3682 8 0,6798 0,5157 0,4377 0,3910 0,3595 0,3362 0,3185 0,3043 0,2926 10 0,6020 0,4450 0,3733 0,3311 0,3029 0,2823 0,2666 0,2541 0,2439 12 0,5410 0,3924 0,3624 0,2880 0,2624 0,2439 0,2299 0,2187 0,2098 15 0,4709 0,3346 0,2758 0,2419 0,2195 0,2034 0,1911 0,1815 0,1736 20 0,3894 0,2705 0,2205 0,1921 0,1735 0,1602 0,1501 0,1422 0,1357 По результатам эксперимента вычисляют коэффициенты модели. Свободный член b0 определяют по формуле 1 b0 = N N åy . j (6) Коэффициенты регрессии, характеризующие линейные эффекты, вычисляют по зависимости 1 bi = N j =1 N åx j= 1 ij yj. (7) Коэффициенты регрессии, характеризующие эффекты взаимодействия, определяют по формуле 42 1 bil = N N åx x ij lj yj, (8) j =1 где i, l - номера факторов; xij, xlj - кодированные значения факторов i и l в j-м эксперименте. Формулы (6), (7), (8) получены в результате использования метода наименьших квадратов. Коэффициенты b0, bi, bij - это оценки теоретических коэффициентов b0, bi, bil регрессии. Оценки, найденные с помощью метода наименьших квадратов, являются наилучшими в том смысле, что они распределены нормально со средними значениями, равными теоретическим коэффициентам, и с наименьшими возможными дисперсиями. Вычислив коэффициенты модели, проверяют их значимость. Проверку значимости коэффициентов можно производить двумя способами: 1. сравнением абсолютной величины коэффициента с доверительным интервалом; 2. с помощью t-критерия Стьюдента. При проверке значимости коэффициентов первым способом для определения доверительного интервала вычисляют дисперсии коэффициентов регрессии. Дисперсию s2{bi} i-го коэффициента определяют по зависимости s 2 {bi } = 1 2 sy . nN (9) Доверительный интервал Dbi - находят по формуле Dbi = ±t τ s{bi }, (10) где tt - табличное значение критерия при принятом уровне значимости и числе сте2 пеней свободы f, с которым определялась дисперсия s y ; при равномерном дублировании экспериментов число степеней свободы находится по зависимости f = (n - 1) N, где N - число экспериментов в матрице планирования, а n-число параллельных экспериментов; s{bi} - ошибка в определении i-го коэффициента регрессии, вычисляе2 s b = + s { } {bi }. Значения t приведены в табл. 14. i мая по формуле Таблица 14 Значения t - критерия при 5% - ном уровне значимости Число степеней свободы Значения t Число степеней свободы Значения t Число степеней свободы Значения t Число степеней свободы Значения t 1 12,71 9 2,26 17 2,11 25 2,06 2 4,30 10 2,23 18 2,10 26 2,06 3 3,18 11 2,20 19 2,09 27 2,05 4 2,78 12 2,18 20 2,09 28 2,05 5 2,57 13 2,16 21 2,08 29 2,05 6 2,45 14 2,14 22 2,07 30 2,04 7 2,37 15 2,13 23 2,07 40 2,02 8 2,30 16 2,12 24 2,06 60 2,00 43 Коэффициент значим, если его абсолютная величина больше доверительного интервала. При проверке значимости коэффициентов вторым способом вычисляют tр - критерий по зависимости tp = bi s{bi } и сравнивают его с табличным tт. Коэффициент значим, если tp>tт для принятого уровня значимости и числа степеней свободы, с которым определялась дисперсия s 2y . Критерий Стьюдента t вычисляют для каждого коэффициента регрессии. Статистически незначимые коэффициенты могут быть исключены из уравнения. После расчета коэффициентов модели и проверки их значимости определяют дисперсию s ад адекватности. Остаточная дисперсия, или дисперсия адекватности, ) характеризует рассеяние эмпирических значений y относительно расчетных y , определенных по найденному уравнению регрессии. Дисперсию адекватности определяют по формуле 2 s ад2 = N ) nå ( y j - y j ) 2 j =1 f = N ) nå ( y j - y j ) 2 j =1 N - (k + 1) , (16.11) где y j - среднее арифметическое значение параметра оптимизации в j-м экспери) менте; y j - значение параметра оптимизации, вычисленное по модели для условий jго опыта; f - число степеней свободы, равное N - (k + 1); k - число факторов. Проверку гипотезы адекватности найденной модели производят по Fкритерию Фишера: sад2 Fp = 2 . sy (12) Если значение FPFт гипотеза адекватности отвергается. Таким образом, обработка результатов эксперимента при равномерном дублировании экспериментов может быть представлена следующей схемой: 1) для каждой строки матрицы планирования по формуле 1 n у = å y jи n и =1 вычисляют среднее арифметическое значение y j параметра оптимизации; 2 2) по формуле (4) определяют дисперсию s j каждого опыта матрицы планирования; 3) используя критерий Кохрена, проверяют гипотезу однородности дисперсий s 2j опытов; 44 4) если дисперсии опытов однородны, то по формуле (5) вычисляют диспер2 s y сию воспроизводимости эксперимента; 5) по формулам (6), (7), (8) определяют коэффициенты уравнения регрессии; 6) по зависимости (9) находят дисперсии s2{bi} коэффициентов регрессии; 7) по формуле (10) устанавливают величину доверительного интервала Dbi; 8) проверяют статистическую значимость коэффициентов регрессии; 2 9) по зависимости (11) определяют дисперсию s ад адекватности; 10) с помощью F-критерия проверяют гипотезу адекватности модели. В заключение необходимо отметить, что использование критериев Кохрена, Стьюдента и Фишера предполагает нормальное распределение результатов эксперимента. Обработка результатов эксперимента при неравномерном дублировании. Результаты отдельных экспериментов иногда получаются ошибочными, и их приходится исключать. Вследствие этого числа параллельных экспериментов оказываются неодинаковыми. Бывают и другие случаи, когда по тем или иным причинам не удается провести одинаковое число параллельных опытов в каждом из основных экспериментов. При неодинаковых числах параллельных экспериментов нарушается ортогональность матрицы планирования и, как следствие, изменяются формулы для определения коэффициентов регрессии и их ошибок. Расчет коэффициентов регрессии и их ошибок при неодинаковых числах параллельных опытов усложняется. Обработка результатов эксперимента при неравномерном дублировании производится по следующей схеме: 1. Для каждой строки матрицы планирования находят y j - среднее арифметическое значение параметра оптимизации 1 yf = nj nj åy ju , u =1 где nj - число параллельных экспериментов в j-й строке матрицы. 2 s 2. Для каждой строки матрицы вычисляют дисперсию j эксперимента; 1 n ( y ju - y f ) 2 . s = å n j - 1 u =1 3. Проверяют с помощью критерия Бартлета гипотезу однородности дисперсий. 2 Для этого подсчитывают дисперсию s y воспроизводимости эксперимента по форj 2 j муле 45 s y2 = 1 æN 2 ö ç åsj f j ÷ N ø å f j è j =1 j =1 2 где fj - число степеней свободы, с которым определялась дисперсия s j i-го эксперимента. После этого определяют величину где N ö æ Q = c -1 ç f lg s y2 - å f j lg s 2j ÷ , j =1 ø è é N 1 æç n 1 1 ö÷ù c = 0 ,4343ê1 + ; f fj . = å å ú ç ÷ j =1 êë 3( N - 1) è j =1 f j f øúû Бартлет показал, что величина Q приближенно подчиняется c2 - распределению с (N - 1) степенями свободы, где N - число сравниваемых дисперсий. 2 Если Q меньше c т (табл. 15) для данного числа (N - 1) степеней свободы и принятого уравнения значимости, то дисперсии однородны, и наоборот. Критерий Бартлета основан на нормальном распределении. Если распределение случайной величины не подчиняется нормальному закону, то проверка однородности дисперсий может привести к ошибочным результатам. Таблица 15 Значения c2 при 5 % - ном уровне значимости Число степеней свободы 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Значения c2 Число степеней свободы 3,84 5,99 7,82 9,49 11,07 12,59 14,07 15,51 16,92 18,31 19,68 21,0 22,4 23,7 25,0 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Значения c2 26,3 27,6 28,9 30,1 31,4 32,7 33,9 35,2 36,4 37,7 38,9 40,1 41,3 42,6 43,8 Рассмотрим применение критерия Бартлета для проверки однородности дисперсий. Матрица планирования предусматривала выполнение четырех экспериментов. Первый был повторен пять раз, второй - шесть, третий и четвертый - по четыре раза. При этом дисперсия первого эксперимента равна 3,5; второго – 4,22; третьего 46 — 5,88; четвертого — 11,36. Необходимо проверить, верна ли гипотеза об однородности дисперсий. 2 Дисперсия s y параметра оптимизации s 2y = 1 æ N 2 çç å s j f N å fj è j =1 ö 3,5 × 4 + 4,22 × 5 + 5,88 × 3 + 11,36 × 3 j ÷÷ = = 5,79. 15 ø j =1 Вычисляем величину с: é 1 æ 1 1 1 1 1 öù c = 0,4343ê1 + ç + + + - ÷ú = 0,485. 3 ( 4 1 ) è 4 5 3 3 15 øû ë Определяем Q: N 1æ ö Q = ç f lg s y2 - å f j lg s 2j ÷ = j= 1 cè ø 1 [15 lg 5,79 - 4 lg 3,5 - 5 lg 4,22 - 3 lg 5,88 - 3 lg 11,36] = 1,37 = 0 ,485 Табличное значение c т для трех степеней свободы (N –1 = 3) и 5 % уровня 2 значимости равно 7,82. Так как Q< c т , то гипотеза однородности дисперсий принимается. 4. Вычисляют коэффициенты bi уравнения регрессии, дисперсии s2{bi} коэффициентов регрессии и ошибки s{bi} в определении коэффициентов. 5. Для каждого коэффициента регрессии находят расчетное значение tкритерия 2 tp = {b } . s{b } i i Сравнивают расчетное значение tp с табличным значением tт критерия. Табличное значение критерия находят для принятого уровня значимости и числа степеней свободы f, которое в рассматриваемом случае определяют по зависимости N N j =1 j =1 f = å fj = å ( n j - 1). Коэффициент значим при tр>tт и незначим при tpFт гипотеза адекватности отвергается. Обработка результатов эксперимента при отсутствии дублирования Обработку результатов эксперимента в этом случае производят по следующей схеме. 2 s 1. Для вычисления дисперсии y воспроизводимости эксперимента выполня- ют несколько параллельных опытов в нулевой точке (в центре плана). При постановке опытов в нулевой точке все факторы находятся на нулевых уровнях. По ре2 s зультатам исследований в центре плана вычисляют дисперсию y воспроизводимо- сти эксперимента 1 én 2ù ( ) s = y y å u ú, n0 - 1 êë u =1 û 2 y где n0 - число параллельных экспериментов в нулевой точке; yu - значение параметра оптимизации в u-м опыте; y - среднее арифметическое значение параметра оптимизации в n0 параллельных экспериментах. 2. Закончив эксперимент, вычисляют коэффициенты модели. Свободный член b0 определяют по формуле 1 b0 = N N åy . j j =1 Коэффициенты регрессии, характеризующие линейные эффекты, вычисляют по зависимости 1 bi = N N åx y . ij j j =1 Коэффициенты регрессии, характеризующие эффекты взаимодействия, определяют по формуле bil = 1 N N åx x ij lj yj, j =1 где i, l - номера факторов; j - номер строки или опыта в матрице планирования; yjзначение параметра оптимизации в j-м опыте; xij, xlj - кодированные значения (±1) факторов i и l в j-м опыте. 3. Проверяют статистическую значимость коэффициентов уравнения регрессии. Проверку значимости коэффициентов можно производить двумя способами: 1) сравнением абсолютной величины коэффициента с доверительным интервалом; 2) с помощью t-критерия. 48 При проверке значимости коэффициентов первым способом для определения доверительного интервала вычисляют дисперсии коэффициентов регрессии по зависимости s 2 {bi } = 1 2 s y, N (13) где s {bi} - дисперсия i-го коэффициента регрессии; N - число строк или опытов в матрице планирования. Из формулы (13) следует, что дисперсии всех коэффициентов равны. Доверительный интервал Dbi определяют по формуле (10). Значение t-критерия, входящего в эту формулу, находят по таблице для принятого уровня значимости и числа степеней свободы f, которое определяют по зависимости f = n0 - 1. Коэффициент регрессии значим, если его абсолютная величина больше доверительного интервала. При проверке значимости коэффициентов вторым способом вычисляют критерий tp b tp = i s{bi } и сравнивают его с табличным tт. Коэффициент значим, если tp>tт для принятого уровня значимости и числа степеней свободы, определенного по формуле f = n0 -1. Критерий Стьюдента t вычисляют для каждого коэффициента регрессии. Статистически незначимые коэффициенты регрессии могут быть исключены из уравнения. 2 2 4. Определяют дисперсию s а , адекватности по формуле N N ) 2 ) (yj - yj ) ( y j - y j )2 å å . sад2 = j =1 = j =1 f N - (k + 1) ) где yj - значение параметра оптимизации в j-м эксперименте; y j - значение параметра оптимизации, вычисленное по модели для условий j-го эксперимента; f - число степеней свободы, которое для линейной модели определяется по зависимости f = N (k + 1), где k - число факторов. 5. Проверяют гипотезу адекватности модели но F-критерию, используя для определения Fр- критерия формулу (16.12). Если FpFт гипотеза адекватности отвергается. В этом случае для получения адекватной модели принимают одно из следующих решений: 1. переходят к планированию второго или более высокого порядка; 2. уменьшают интервалы варьирования и ставят новый эксперимент, повторяя эти действия до получения адекватной линейной модели. Если линейная модель адекватна, то переходят к методу крутого восхождения. Необходимо заметить, что крутое восхождение эффективно тогда, когда все коэффициенты при факторах значимы. Незначимость некоторых коэффициентов может получиться вследствие: 1. неудачно выбранных интервалов варьирования; 49 2. включения факторов, не влияющих на параметр оптимизации; 3. большой ошибки эксперимента. Принятие решения в данной ситуации зависит от того, какая из трех указанных причин имеет место. Если имеет место первая причина, то изменяют интервалы варьирования по незначимым факторам и ставят новую серию экспериментов. Если принята вторая, то не влияющие факторы стабилизируют и исключают из экспериментов. Если принята третья причина, то увеличивают число параллельных экспериментов. Увеличение их числа приводит к уменьшению дисперсии коэффициентов и величины доверительного интервала, в результате все или часть коэффициентов могут оказаться значимыми. Возможен случай, когда все коэффициенты, кроме b0, незначимы, а модель адекватна. Такая ситуация чаще всего возникает из-за слишком узких интервалов варьирования или вследствие большой ошибки эксперимента. В этой ситуации возможны два решения: 1) расширение интервалов варьирования или 2) повышение точности эксперимента путем улучшения методики проведения и увеличения числа параллельных экспериментов. Лекция 6 БАЗОВЫЕ ПОНЯТИЯ И ОПЕРАЦИИ ОБРАБОТКИ ДАННЫХ (ЭД) ЭКСПЕРИМЕНТАЛЬНЫХ 1. Эмпирическая функция распределения Методы обработки ЭД опираются на базовые понятия теории вероятностей и математической статистики. К их числу относятся понятия генеральной совокупности, выборки, эмпирической функции распределения. Под генеральной совокупностью понимают все возможные значения параметра, которые могут быть зарегистрированы в ходе неограниченного по времени наблюдения за объектом. Такая совокупность состоит из бесконечного множества элементов. В результате наблюдения за объектом формируется ограниченная по объему совокупность значений параметра x1, x2, …, xn. С формальной точки зрения такие данные представляют собой выборку из генеральной совокупности. Будем считать, что выборка содержит полные наработки до системных событий. Наблюдаемые значения xi называют вариантами, а их количество – объемом выборки n. Для того чтобы по результатам наблюдения можно было делать какиелибо выводы, выборка должна быть репрезентативной (представительной), т. е. правильно представлять пропорции генеральной совокупности. Это требование выполняется, если объем выборки достаточно велик, а каждый элемент генеральной совокупности имеет одинаковую вероятность попасть в выборку. Пусть в полученной выборке значение x1 параметра наблюдалось n1 раз, значение x2 – n2 раз, значение xk – nk раз, n1+n2+ … +nk=n. Совокупность значений, записанных в порядке их возрастания, называют вариационным рядом, величины ni – частотами, а их отношения к объему выборки wi =ni/n – относительными частотами (частостями). Очевидно, что сумма относительных частот равна единице. 50 Под распределением понимают соответствие между наблюдаемыми вариантами и их частотами или частостями. Пусть nx – количество наблюдений, при которых случайные значения параметра Х меньше x. Частость события X40) в целях удобства обработки и хранения сведений прибегают к группированию ЭД в интервалы. Количество интервалов следует выбрать так, чтобы в необходимой мере отразилось разнообразие значений параметра в совокупности и в то же время закономерность распределения не искажалась случайными колебаниями частот по отдельным разрядам. Существуют нестрогие рекомендации по выбору количества у и размера h таких интервалов, в частности: в каждом интервале должно находиться не менее 5 – 7 элементов. В крайних разрядах допустимо всего два элемента; количество интервалов не должно быть очень большим или очень маленьким. Минимальное значение интервалов должно быть не менее 6 – 7. При объеме выборки, не превышающем несколько сотен элементов, количество интервалов задают в пределах от 10 до 20. Для очень большого объема выборки (n>1000) количество интервалов может превышать указанные значения. Некоторые исследователи рекомендуют пользоваться соотношением y=1,441ln(n)+1; при относительно небольшой неравномерности длины интервалов удобно выбирать одинаковыми и равными величине h=(xmax – xmin)/у, где xmax – максимальное и xmin – минимальное значение параметра. При существенной неравномерности закона распределения длины интервалов можно задавать меньшего размера в области быстрого изменения плотности распределения; при значительной неравномерности лучше в каждый разряд назначать примерно одинаковое количество элементов выборки. Тогда длина конкретного интервала будет определяться крайними значениями элементов выборки, сгруппированными в этот интервал, т.е. будет различна для разных интервалов (в этом случае при построении гистограммы нормировка по длине интервала обязательна - в противном случае высота каждого элемента гистограммы будет одинакова). Группирование результатов наблюдений по интервалам предусматривает: определение размаха изменений параметра х; выбор количества интервалов и их величины; подсчет для каждого i-го интервала [xi–xi+1] частоты ni или относительной частоты (частости ni) попадания варианты в интервал. В результате формируется представление ЭД в виде интервального или статистического ряда. 52 Графически статистический ряд отображают в виде гистограммы, полигона и ступенчатой линии. Часто гистограмму представляют как фигуру, состоящую из прямоугольников, основаниями которых служат интервалы длиною h, а высоты равны соответствующей частости. Однако такой подход неточен. Высоту i-го прямоугольника zi следует выбрать равной ni/(nh). Такую гистограмму можно интерпретировать как графическое представление эмпирической функции плотности распределения fn(x), в ней суммарная площадь всех прямоугольников составит единицу. Гистограмма помогает подобрать вид теоретической функции распределения для аппроксимации ЭД. Полигоном называют ломаную линию, отрезки которой соединяют точки с координатами по оси абсцисс, равными серединам интервалов, а по оси ординат – соответствующим частостям. Эмпирическая функция распределения отображается ступенчатой ломаной линией: над каждым интервалом проводится отрезок горизонтальной линии на высоте, пропорциональной накопленной частости в текущем интервале. Накопленная частость равна сумме всех частостей, начиная с первого и до данного интервала включительно. Пример 2. Имеются результаты регистрации значений затухания сигнала xi на частоте 1000 Гц коммутируемого канала телефонной сети. Эти значения, измеренные в дБ, в виде вариационного ряда представлены в табл. 3. Необходимо построить статистический ряд. Таблица 3 i xi i xi i xi i xi 1 25,79 12 26,85 23 27,40 34 28,37 2 25,98 13 26,90 24 27,49 35 28,38 3 25,98 14 26,91 25 27,64 36 28,50 4 26,12 15 26,96 26 27,66 37 28,63 5 26,13 16 27,02 27 27,71 38 28,67 6 26,49 17 27,11 28 27,78 39 28,90 7 26,52 18 27,19 29 27,89 40 28,99 8 26,60 19 27,21 30 27,89 41 28,99 9 26,66 20 27,28 31 28,01 42 29,03 10 26,69 21 27,30 32 28,10 43 29,12 11 26,74 22 27,38 33 28,11 44 29,28 Решение. Количество разрядов статистического ряда следует выбрать минимальным, чтобы обеспечить достаточное количество попаданий в каждый из них, возьмем y = 6. Определим размер разряда h = (xmax – xmin)/у =(29,28 – 25,79)/6 = 0,58. Сгруппируем наблюдения по разрядам, табл. 4. Таблица 4 i xi ni Fi =ni/n fi = ni/h 1 25,79 5 0,11 4 0,196 2 26,37 9 0,205 0,353 3 26,95 10 0,227 0,392 4 27,5 3 9 0,205 0,353 5 28,12 5 0,11 4 0,196 6 28,70 6 0,136 0,235 На основе статистического ряда построим гистограмму, рис. 2, и график эмпирической функции распределения, рис. 3. График эмпирической функции распределения, рис. 3, отличается от графика, представленного на рис. 1 равенством шага изменения варианты и величиной шага 53 приращения функции (при построении по вариационному ряду шаг приращения кратен 1/ n , а по статистическому ряду – зависит от частости в конкретном разряде). Рис.2. Гистограмма распределения Рис. 3. Эмпирическая функция распределения Рассмотренные представления ЭД являются исходными для последующей обработки и вычисления различных параметров. 2. Оценки параметров распределения и их свойства Значение параметра, вычисленное по ограниченному объему ЭД, является случайной величиной, т. е. значение такой величины от выборки к выборке может меняться заранее не предвиденным образом. Следовательно, в результате обработки ЭД определяется не значение параметра q , а только лишь его приближенное значе~ ние q – статистическая оценка параметра. Получить статистическую оценку параметра теоретического распределения означает найти функцию от имеющихся результатов наблюдения, которая и даст приближенное значение искомого параметра. Различают два вида оценок – точечные и интервальные. Точечными называют такие оценки, которые характеризуются одним числом. При малых объемах выборки точечные оценки могут значительно отличаться от истинных значений параметров, поэтому их применяют при большом объеме выборки. Интервальные оценки задаются двумя числами, определяющими вероятный диапазон возможного значения параметра. Эти оценки применяются для малых и для больших выборок. Рассмотрим вначале точечные оценки. Применительно к каждому оцениваемому параметру закона распределения генеральной совокупности существует множество функций, позволяющих вычислить искомые значения. Например, оценку математического ожидания можно вычислить, взяв среднее арифметическое выборочных значений, половину суммы крайних членов вариационного ряда, средний член выборки и т.д. Указанные функции отличаются качеством оценок и трудоемкостью реализации. Качество оценок характеризуется такими свойствами, как состоятельность, несмещенность, эффективность и достаточность. ~ Состоятельность характеризует сходимость по вероятности оценки q к истинному значению параметра q при неограниченном увеличении объема выборки n.: 54 · Смещение оценки должно стремиться к нулю при объёме выборки, стре~ мящейся к бесконечности, т.е. M q i - qi = ji ® 0 ; ~ D q i ® 0 при n ® ¥ . Оценка является состоятельной, если она удовлетворяет закону больших чисел, т.е. с увеличением объёма выборки значение оценки стремится к значению параметра генеральной совокупности по вероятности: ~ P q i - q i < 0 ® 1 , при n ® ¥ . ( ) ( ) ( ) · Для определения состоятельности оценки достаточно выполнения двух условий Свойство состоятельности проявляется при неограниченном увеличении n, а при небольших объемах ЭД наличие этого свойства еще недостаточно для применения оценки. Несмещенность характеризует отсутствие систематических (в среднем) отклонений оценки от параметра при любом конечном, в том числе и малом, объеме вы~ борки, т. е. M(q ) =q . Использование статистической оценки, математическое ожидание которой не равно оцениваемому параметру, приводит к систематическим ошибкам. Не всегда наличие смещения плохо. Оно может быть существенно меньше погрешности регистрации значений параметра или давать дополнительную гарантию выполнения тре~ бований к значению параметра (если даже при положительном смещении оценка q меньше предельно допустимого значения, то несмещенное значение тем более будет отвечать этому условию). В таких ситуациях допустимо применение смещенных оценок, если они вычисляются проще, чем несмещенные. Но даже несмещенная оценка может быть удалена от истинного значения. Эффективность характеризует разброс случайных значений оценки около истинного значения параметра. Среди всех оценок следует выбрать ту, значения которой теснее сконцентрированы около оцениваемого параметра. Для многих применяемых способов оценивания выборочные распределения параметров асимптотически нормальны, поэтому часто мерой эффективности служит дисперсия оценки. ~ Оценка q является эффективной для параметра q , если она имеет наименьшую из возможных дисперсию при заданном объёме выборки п. Достаточность характеризует полноту использования информации, содержа~ щейся в выборке. Другими словами, оценка q будет достаточной, если все другие независимые оценки на основе данной выборки не дают дополнительной информации об оцениваемом параметре. Эффективная оценка обязательно является и достаточной. Рассмотренные свойства применимы также и к ЭД, которые характеризуются многомерными распределениями вероятностей. Подходы к формированию оценок разработаны в теории несмещенных оценок, предложенной А. Н. Колмогоровым и С. Рао. В данной теории предполагается известным с точностью до параметра T вид функции плотности распределения наблюдаемой величины f(x, Т). Вид распределения устанавливается исходя из априор55 ных соображений, например, на основе общепринятых суждений о характере безотказной работы технических средств. Тогда задача сводится к нахождению такой функции от результатов наблюдений, которая дает несмещенную и эффективную оценку. 23. Оценки моментов и квантилей распределения Для характеристики эмпирического распределения можно использовать оценки центральных и начальных моментов. Применение находят моменты до четвертого порядка включительно, так как точность выборочных моментов резко падает с увеличением их порядка, в частности, дисперсия начальных моментов порядка р зависит от моментов порядка 2р. Она становится значительной для моментов высокого порядка даже при больших объемах выборки. Выборочные значения моментов определяют непосредственно по выборке или по сгруппированным данным. Пусть Х – изучаемый количественный признак, x1 , x2 ,......, xk имеющаяся выборка объёма n = n1 + n2 + .... + nk . Для произвольного действительного числа х обозначим пх – количество вариант выборки, меньших числа х. n Относительная частота появления события (Х<х) равна x . n n Функция F * ( x ) = x , определяющая для каждого x Î n x относительную частоту n события (Х<х), называется эмпирической функцией распределения или выборочной функцией распределения. Функция распределения F ( x ) генеральной совокупности называется теоретической функцией распределения. Теоретическая функция распределения характеризует вероятность события (Х<х), а эмпирическая функция характеризует относительную частоту данного события. Эмпирическая функция распределения обладает всеми свойствами теоретической функции распределения. Рассмотрим выборочные числовые характеристики. Выборочным средним называется число k n x в = å xi i . i =1 n Выборочной дисперсией называется число k 2 n i . Dв = å xi - xв i =1 n Начальный выборочный момент р-го порядка - это число k p n mkB = å ( xi ) × i . i =1 n Центральный выборочный момент р-го порядка – это число k p n B mkC = å xi - x в × i . i =1 n ( ) ( ) 56 Важно заметить, что для интервального вариационного ряда в приведённых выше формулах вместо значений вариант xi берётся середина соответствующего интервала ~ xi , а число ni в этом случае есть количество наблюдённых значений, попавших в i-ый интервал. . Главное назначение средних величин (оценок начальных моментов и в первую очередь первого момента распределения) состоит в их обобщающей функции. Это обобщение позволяет заменить множество различных индивидуальных значений показателя средней величиной, характеризующей всю однородную совокупность. Иначе говоря, средняя величина является типической характеристикой варианты в конкретной выборке. Иногда средняя величина обобщает и неоднородные совокупности данных. Например, может применяться такой показатель как среднее количество обработанных запросов на сервере в течение суток, хотя очевидно, что дневная загрузка сервера сильно отличается от загрузки в ночное время. Указанный показатель имеет смысл для оценки ресурса накопителей на жестких дисках. Каждый элемент ЭД формируется под влиянием как общих закономерностей, так и особых условий и случайных событий. Следовательно, в обработке ЭД большой интерес представляют вопросы оценки величин, характеризующих вариацию значений параметра у разных объектов или у одного и того же объекта в разные моменты времени. Вариацией какого-либо параметра (показателя) в совокупности наблюдений называется различие его значений у разных элементов этой совокупности. Вариация является характерным свойством большинства информационных параметров АСОИУ. Именно это свойство является объектом исследования большинства методов обработки ЭД. Для характеристики вариации нет единого показателя, в этих целях применяются моменты распределения выше первого, производные от них величины, размах выборки, квантили и др. Лекция 7 ПРОВЕРКА СТАТИСТИЧЕСКИХ ГИПОТЕЗ 1. Сущность задачи проверки статистических гипотез Статистическая гипотеза представляет собой некоторое предположение о законе распределения случайной величины или о параметрах этого закона, формулируемое на основе выборки. Примерами статистических гипотез являются предположения: генеральная совокупность распределена по экспоненциальному закону; математические ожидания двух экспоненциально распределенных выборок равны друг другу. В первой из них высказано предположение о виде закона распределения, а во второй – о параметрах двух распределений. Гипотезы, в основе которых нет никаких допущений о конкретном виде закона распределения, называют непараметрическими, в противном случае – параметрическими. Гипотезу, утверждающую, что различие между сравниваемыми характеристиками отсутствует, а наблюдаемые отклонения объясняются лишь случайными ко57 лебаниями в выборках, на основании которых производится сравнение, называют нулевой (основной) гипотезой и обозначают Н0. Наряду с основной гипотезой рассматривают и альтернативную (конкурирующую, противоречащую) ей гипотезу Н1. И если нулевая гипотеза будет отвергнута, то будет иметь место альтернативная гипотеза. Различают простые и сложные гипотезы. Гипотезу называют простой, если она однозначно характеризует параметр распределения случайной величины. Например, если l является параметром экспоненциального распределения, то гипотеза Н0 о равенстве l=10 – простая гипотеза. Сложной называют гипотезу, которая состоит из конечного или бесконечного множества простых гипотез. Сложная гипотеза Н0 о неравенстве l>10 состоит из бесконечного множества простых гипотез Н0 о равенстве l=bi, где bi – любое число, большее 10. Гипотеза Н0 о том, что математическое ожидание нормального распределения равно двум при неизвестной дисперсии, тоже является сложной. Сложной гипотезой будет предположение о распределении случайной величины Х по нормальному закону, если не фиксируются конкретные значения математического ожидания и дисперсии. Проверка гипотезы основывается на вычислении некоторой случайной величины – критерия, точное или приближенное распределение которого известно. Обозначим эту величину через z, ее значение является функцией от элементов выборки z=z(x1, x2, …, xn). Процедура проверки гипотезы предписывает каждому значению критерия одно из двух решений – принять или отвергнуть гипотезу. Тем самым все выборочное пространство и соответственно множество значений критерия делятся на два непересекающихся подмножества S0 и S1. Если значение критерия z попадает в область S0, то гипотеза принимается, а если в область S1, то гипотеза отклоняется. Множество S0 называется областью принятия гипотезы или областью допустимых значений, а множество S1 – областью отклонения гипотезы или критической областью. Выбор одной области однозначно определяет и другую область. Принятие или отклонение гипотезы Н0 по случайной выборке соответствует истине с некоторой вероятностью и, соответственно, возможны два рода ошибок. Ошибка первого рода возникает с вероятностью a тогда, когда отвергается верная гипотеза Н0 и принимается конкурирующая гипотеза Н1. Ошибка второго рода возникает с вероятностью b в том случае, когда принимается неверная гипотеза Н0, в то время как справедлива конкурирующая гипотеза Н1. Доверительная вероятность – это вероятность не совершить ошибку первого рода и принять верную гипотезу Н0. Вероятность отвергнуть ложную гипотезу Н0 называется мощностью критерия. Следовательно, при проверке гипотезы возможны четыре варианта исходов, табл. 1. Таблица 1. Гипотеза Н0 Верна Решение Вероятность Примечание Принимается 1- a Доверительная вероятность Отвергается a Вероятность ошибки первого рода 58 Неверна Принимается b Отвергается 1-b Вероятность ошибки второго рода Мощность критерия Например, рассмотрим случай, когда некоторая несмещенная оценка параметра q вычислена по выборке объема n, и эта оценка имеет плотность распределения f( q ), рис. 1. Рис. 1. Области принятия и отклонения гипотезы Предположим, что истинное значение оцениваемого параметра равно Т. Если рассматривать гипотезу Н0 о равенстве q =Т, то насколько велико должно быть различие между q и Т, чтобы эту гипотезу отвергнуть. Ответить на данный вопрос можно в статистическом смысле, рассматривая вероятность достижения некоторой заданной разности между q и Т на основе выборочного распределения параметра q . Целесообразно полагать одинаковыми значения вероятности выхода параметра q за нижний и верхний пределы интервала. Такое допущение во многих случаях позволяет минимизировать доверительный интервал, т.е. повысить мощность критерия проверки. Суммарная вероятность выхода параметра q за пределы интервала с границами q 1–a/2 и q a/2, составляет величину a . Эту величину следует выбирать такой, чтобы выход за пределы интервала был маловероятен. Если оценка параметра попала в заданный интервал, то нет основания подвергать сомнению проверяемую гипотезу. Гипотезу равенства q =Т можно принять. Но если после получения выборки окажется, что оценка выходит за установленные пределы, то в этом случае есть серьезные основания отвергнуть гипотезу Н0. Отсюда следует, что вероятность допустить ошибку первого рода равна a (равна уровню значимости критерия). Если предположить, например, что истинное значение параметра в действительности равно Т+d , то согласно гипотезе Н0 о равенстве q =Т – вероятность того, что оценка параметра q попадет в область принятия гипотезы, составит b , рис. 1. При заданном объеме выборки вероятность совершения ошибки первого рода можно уменьшить, снижая уровень значимости a . Однако при этом увеличивается вероятность ошибки второго рода b (снижается мощность критерия). Аналогичные 59 рассуждения можно провести для случая, когда истинное значение параметра равно Т–d. Единственный способ уменьшить обе вероятности состоит в увеличении объема выборки (плотность распределения оценки параметра при этом становится более "узкой"). При выборе критической области руководствуются правилом Неймана – Пирсона: следует так выбирать критическую область, чтобы вероятность a была мала, если гипотеза верна, и велика в противном случае. Однако выбор конкретного значения a относительно произволен. Употребительные значения лежат в пределах от 0,001 до 0,2. В целях упрощения ручных расчетов составлены таблицы интервалов с границами q 1–a/2 и q a/2 для типовых значений a и различных способов построения критерия. При выборе уровня значимости необходимо учитывать мощность критерия при альтернативной гипотезе. Иногда большая мощность критерия оказывается существеннее малого уровня значимости, и его значение выбирают относительно большим, например 0,2. Такой выбор оправдан, если последствия ошибок второго рода более существенны, чем ошибок первого рода. Например, если отвергнуто правильное решение "продолжить работу пользователей с текущими паролями", то ошибка первого рода приведет к некоторой задержке в нормальном функционировании системы, связанной со сменой паролей. Если же принято решения не менять пароли, несмотря на опасность несанкционированного доступа посторонних лиц к информации, то эта ошибка повлечет более серьезные последствия. В зависимости от сущности проверяемой гипотезы и используемых мер расхождения оценки характеристики от ее теоретического значения применяют различные критерии. К числу наиболее часто применяемых критериев для проверки гипотез о законах распределения относят критерии хи-квадрат Пирсона, Колмогорова, Мизеса, Вилкоксона, о значениях параметров – критерии Фишера, Стьюдента. Лекция 8 1. Типовые распределения При проверке гипотез широкое применение находит ряд теоретических законов распределения. Наиболее важным из них является нормальное распределение. С ним связаны распределения хи-квадрат, Стьюдента, Фишера, а также интеграл вероятностей. Для указанных законов функции распределения аналитически не представимы. Значения функций определяются по таблицам или с использованием стандартных процедур пакетов прикладных программ. Указанные таблицы обычно построены в целях удобства проверки статистических гипотез в ущерб теории распределений – они содержат не значения функций распределения, а критические значения аргумента z( a ). Для односторонней критической области z( a )=z1–a, т.е. критическое значение аргумента z( a ) соответствует квантили z1–a уровня 1– a , рис 3, так как ¥ z (a ) z -¥ . ò f ( z )dz = a = 1 - ò f ( z )dz 60 Рис. 3. Односторонняя критическая область Для двусторонней критической области, с уровнем значимости a , размер левой области a 2, правой a 1 ( a 1+ a 2= a ), рис.4. Значения z( a 2) и z( a 1) связаны с квантилями распределения соотношениями z( a 1)= z1–a1, z( a 2)= za2, так как z (a1 ) ò f ( z )dz = 1 - a1 , -¥ z (a 2 ) ò f ( z )dz = 1 - a 2 -¥ Для симметричной функции плотности распределения f(z) критическую область выбирают из условия a 1= a 2= a /2 (обеспечивается наибольшая мощность критерия). В таком случае левая и правая границы будут равны |z( a /2)|. Рис. 4. Двусторонняя критическая область Нормальное распределение Этот вид распределения является наиболее важным в связи с центральной предельной теоремой теории вероятностей: распределение суммы независимых случайных величин стремится к нормальному с увеличением их количества при произвольном законе распределения отдельных слагаемых, если слагаемые обладают конечной дисперсией. Кроме того, А.М. Ляпунов доказал, что распределение параметра стремится к нормальному, если на параметр оказывает влияние большое количество факторов и ни один из них не является превалирующим. Функция плотности нормального распределения é ( x - m1 )2 ù 1 f (x ) = exp êú 2m2 û 2pm2 ë 61 – унимодальная, симметричная, аргумент х может принимать любые действительные значения, рис. 3.5. Рис. 5. Плотность нормального распределения Функция плотности нормального распределения стандартизованной величины u имеет вид é u2 ù 1 f (u ) = exp ê- ú 2p ë 2û Вычисление значений функции распределения Ф(u) для стандартизованного неотрицательного аргумента u ( u ³ 0 ) можно произвести с помощью полинома наилучшего приближения Ф(u)= 1– 0,5(1 + 0,196854u + 0,115194u2 + 0,000344u3 + 0,019527u4) – 4 Такая аппроксимация обеспечивает абсолютную ошибку не более 0,00025. Для вычисления Ф(u) в области отрицательных значений стандартизованного аргумента u (u<0) следует воспользоваться свойством симметрии нормального распределения Ф(u) = 1 – Ф(– u). Интеграл вероятностей связан с функцией нормального распределения стандартизованной величины u соотношением Ф(u) = 0,5 + F(u). Распределение хи-квадрат Распределению хи-квадрат ( c 2 -распределению) с k степенями свободы соответствует распределение суммы n c 2 = å ui2 i =1 квадратов n стандартизованных случайных величин ui , каждая из которых распределена по нормальному закону, причем k из них независимы, n ³ k . Функция плотности распределения хи-квадрат с k степенями свободы é k æ k öù f ( x ) = ê 2 2 Г ç ÷ú è 2 øû ë -1 k (x )2 -1 e - k 2 , x ³ 0, (4) где Г(k/2) – гамма-функция. 62 Функция плотности при k, равном одному или двум, – монотонная, а при k>2 – унимодальная, несимметричная, рис. 6. Рис. 6. Плотность распределения хи-квадрат Математическое ожидание и дисперсия величины c 2 равны соответственно k и 2k . Распределение хи-квадрат является частным случаем более общего гаммараспределения, а величина, равная корню квадратному из хи-квадрат с двумя степенями свободы, подчиняется распределению Рэлея. С увеличением числа степеней свободы (k>30) распределение хи-квадрат приближается к нормальному распределению с математическим ожиданием k и дисперсией 2k. Распределение Стьюдента Распределение Стьюдента (t-распределение, предложено в 1908 г. английским статистиком В. Госсетом, публиковавшим научные труды под псевдонимом Student) характеризует распределение случайной величины t= (u 2 1 u0 + u 22 + ... + u k2 ) k , где u0, u1, …, uk взаимно независимые нормально распределенные случайные величины с нулевым средним и конечной дисперсией. Аргумент t не зависит от дисперсии слагаемых. Функция плотности распределения Стьюдента f (t ) = Г [(k + 1) 2] é t ù 1+ pk Г (k 2 ) êë k úû 2 - ( k +1) 2 (5) Величина k характеризует количество степеней свободы. Плотность распределения – унимодальная и симметричная функция, похожая на нормальное распределение, рис. 7. 63 Рис. 7. Плотность распределения Стьюдента Область изменения аргумента t от минус до плюс бесконечности. Математическое ожидание и дисперсия равны 0 и k/(k–2) соответственно, при k>2. По сравнению с нормальным распределение Стьюдента более пологое, оно имеет меньшую дисперсию. Это отличие заметно при небольших значениях k, что следует учитывать при проверке статистических гипотез (критические значения аргумента распределения Стьюдента превышают аналогичные показатели нормального распределения). Таблицы распределения содержат значения для односторонней (пределы интегрирования от r(k; a ) до бесконечности) ¥ ò f (t )dt = a r ( k ,a ) или двусторонней (пределы интегрирования от – r(k; a ) до r(k; a ) r ( k ,a ) ò f (t )dt = a -r ( k ,a ) критической области. Распределение Стьюдента применяется для описания ошибок выборки при k<30. При k, превышающем 100, данное распределение практически соответствует нормальному, для значений k из диапазона от 30 до 100 различия между распределением Стьюдента и нормальным распределением составляют несколько процентов. Поэтому относительно оценки ошибок малыми считаются выборки объемом не более 30 единиц, большими – объемом более 100 единиц. При аппроксимации распределения Стьюдента нормальным распределением для односторонней критической области вероятность a Р{t>t(k; a )} = u1–a(0, k /(k–2)), где u1–a(0, k/(k–2)) – квантиль нормального распределения. Аналогичное соотношение можно составить и для двусторонней критической области. Распределение Фишера 64 Распределению Р.А. Фишера (F-распределению Фишера – Снедекора) подчиняется случайная величина Х = [(y1/k1)/(y2/k2)], равная отношению двух случайных величин у1 и у2, имеющих хи-квадрат распределение с k1 и k2 степенями свободы. Область изменения аргумента х от 0 до бесконечности. Плотность распределения ék ù f (x ) = ê 1 ú ë k2 û k1 2 Г [(k1 + k 2 ) 2] ( k1 -2 ) 2 é k1 ù xú x ê1 + Г (k1 2 )Г (k 2 2 ) ë k2 û - ( k1 + k2 ) 2 (6) В этом выражении k1 обозначает число степеней свободы величины y1 с большей дисперсией, k2 – число степеней свободы величины y2 с меньшей дисперсией. Плотность распределения – унимодальная, несимметричная, рис. 8. Рис. 8. Плотность распределения Фишера Математическое ожидание случайной величины х m1 = k2 /(k2–2) при k2>2, дисперсия т2 = [2k22(k1 + k2 –2)]/[k1(k2 –2)2(k2–4)] при k 2 > 4. При k1>30 и k2>30 величина х распределена приближенно нормально: с центром распределения (k1–k2)/(2k1k2) и дисперсией (k1+k2)/(2k1k2). 3. Проверка гипотез о законе распределения Обычно сущность проверки гипотезы о законе распределения ЭД заключается в следующем: Имеется выборка ЭД фиксированного объема, выбран или известен вид закона распределения генеральной совокупности. Необходимо оценить по этой выборке параметры закона, определить степень согласованности ЭД и выбранного закона распределения, в котором параметры заменены их оценками. Пока не будем касаться способов нахождения оценок параметров распределения, а рассмотрим только вопрос проверки согласованности распределений с использованием наиболее употребительных критериев. 3.1. Критерий хи-квадрат К. Пирсона Использование этого критерия основано на применении такой меры (статистики) расхождения между теоретическим F(x) и эмпирическим распределением Fп(x), ко65 торая приближенно подчиняется закону распределения c 2 . Гипотеза Н0 о согласованности распределений проверяется путем анализа распределения этой статистики. Применение критерия требует построения статистического ряда. Итак, пусть выборка представлена статистическим рядом с количеством разрядов y. Наблюдаемая частота попаданий в i-й разряд ni. В соответствии с теоретическим законом распределения ожидаемая частота попаданий в i-й разряд составляет Fi. Разность между наблюдаемой и ожидаемой частотой составит величину (ni–Fi). Для нахождения общей степени расхождения между F(x) и Fп(x) необходимо подсчитать взвешенную сумму квадратов разностей по всем разрядам статистического ряда 2 y (n - F ) 2 i i . c =å (7) i =1 Fi Величина c 2 при неограниченном увеличении n имеет распределение хи-квадрат (асимптотически распределена как хи-квадрат). Это распределение зависит от числа степеней свободы k, т.е. количества независимых значений слагаемых в выражении (7). Число степеней свободы равно числу y минус число линейных связей, наложенных на выборку. Одна связь существует в силу того, что любая частота может быть вычислена по совокупности частот в оставшихся y–1 разрядах. Кроме того, если параметры распределения неизвестны заранее, то имеется еще одно ограничение, обусловленное подгонкой распределения к выборке. Если по выборке определяются f параметров распределения, то число степеней свободы составит k = y–f–1. Очевидно, что чем меньше расхождение между теоретическими и эмпирическими частотами, тем меньше величина критерия. Область принятия гипотезы Н0 определяется условием c 2 < c 2 (k ,a ), где c 2 (k ,a ) – критическая точка распределения хиквадрат с уровнем значимости a . Вероятность ошибки первого рода равна a , вероятность ошибки второго рода четко определить нельзя, потому что существует бесконечно большое множество различных способов несовпадения распределений. Мощность критерия зависит от количества разрядов и объема выборки. Критерий рекомендуется применять при n>200, допускается применение при n>40, именно при таких условиях критерий состоятелен (как правило, отвергает неверную нулевую гипотезу). Теоретическое значение вероятности Fi попадания случайной величины в i-й интервал определяется по формуле xi æ ( x - m1 )2 ö 1 ÷dx . Fi = ò expç 2s 2 ÷ø s 2p xi -1 çè (8) Для нормального закона возможные значения случайной величины лежат в диапазоне от минус до плюс бесконечности, поэтому при расчетах оценок вероятностей крайний левый и крайний правый интервалы расширяются до минус и плюс бесконечности соответственно. Вычислить значения функции нормального распределения можно, воспользовавшись стандартными функциями табличного процессора или полиномом наилучшего приближения. 66 3.2. Критерий А.Н. Колмогорова Для применения критерия А.Н. Колмогорова ЭД требуется представить в виде вариационного ряда (ЭД недопустимо объединять в разряды). В качестве меры расхождения между теоретической F(x) и эмпирической Fn(x) функциями распределения непрерывной случайной величины Х используется модуль максимальной разности d n = max F ( x ) - Fn ( x ) (9) А.Н. Колмогоров доказал, что какова бы ни была функция распределения F(x) величины Х при неограниченном увеличении количества наблюдений n функция распределения случайной величины d n n асимптотически приближается к функции распределения ( ) K (l ) = P d n < l = å (- 1) exp(- 2k 2 l2 ) . ¥ k k = -¥ (10) Иначе говоря, критерий А.Н. Колмогорова характеризует вероятность того, что величина d n n не будет превосходить параметр l для любой теоретической функции распределения. Уровень значимости a выбирается из условия P d n n > l = a = 1 - K (l ), (11) в силу предположения, что почти невозможно получить это равенство, когда существует соответствие между функциями F(x) и Fп(x). Критерий А.Н. Колмогорова позволяет проверить согласованность распределений по малым выборкам, он проще критерия хи-квадрат, поэтому его часто применяют на практике. Но требуется учитывать два обстоятельства: 1. В соответствии с условиями его применения необходимо пользоваться следующим соотношением d n = max(d n+ , d n- ), ( где d n+ = max 1£ j £n ) j j -1 - F ( x ) ; d n- = max F ( x ) . 1£ j £ n n n (12) 2. Условия применения критерия предусматривают, что теоретическая функция распределения известна полностью – известны вид функции и значения ее параметров. На практике параметры обычно неизвестны и оцениваются по ЭД. Но критерий не учитывает уменьшение числа степеней свободы при оценке параметров распределения по исходной выборке. Это приводит к завышению значения вероятности соблюдения нулевой гипотезы, т.е. повышается риск принять в качестве правдоподобной гипотезу, которая плохо согласуется с ЭД (повышается вероятность совершить ошибку второго рода). В качестве меры противодействия такому выводу следует увеличить уровень значимости a , приняв его равным 0,1 – 0,2, что приведет к уменьшению зоны допустимых отклонений. 67 Пример. Проверить с помощью критерия А.Н. Колмогорова гипотезу о том, что ЭД, представленные в табл. 3, подчиняются нормальному распределению при уровне значимости a=0,1. Решение. Исходные данные и результаты вычислений сведены в табл. 3. Необходимые вычисления можно провести с использованием табличного процессора: значение эмпирической функции распределения Fn ( xi ) ; i = 1,44 ; значения теоретической функции F(xi) – это значение функции нормального распределения в точке xi. Таблица 3 i xi Fn(xi ) F(xi ) dn+ dni xi Fn(xi ) F(xi ) dn+ dni xi Fn(xi ) F(xi ) dn+ dni xi Fn(xi ) F(xi ) dn+ dn- 1 25,79 0,023 0,036 0,014 0,036 12 26,85 0,273 0,246 0,027 0,004 23 27,40 0,523 0,456 0,067 0,044 34 28,37 0,773 0,817 0,044 0,067 2 25,98 0,046 0,055 0,009 0,032 13 26,90 0,296 0,263 0,032 0,010 24 27,49 0,546 0,492 0,053 0,031 35 28,38 0,795 0,819 0,024 0,046 3 25,98 0,068 0,055 0,013 0,010 14 26,91 0,318 0,267 0,051 0,028 25 27,64 0,568 0,555 0,013 0,010 36 28,50 0,818 0,851 0,032 0,055 4 26,12 0,091 0,073 0,018 0,005 15 26,96 0,341 0,284 0,057 0,034 26 27,66 0,591 0,561 0,030 0,007 37 28,63 0,841 0,879 0,038 0,061 5 26,13 0,114 0,075 0,038 0,016 16 27,02 0,364 0,305 0,059 0,036 27 27,71 0,614 0,583 0,031 0,008 38 28,67 0,864 0,888 0,024 0,047 6 26,49 0,136 0,144 0,008 0,031 17 27,11 0,386 0,337 0,050 0,027 28 27,78 0,636 0,610 0,026 0,003 39 28,90 0,886 0,928 0,042 0,064 7 26,52 0,159 0,151 0,008 0,014 18 27,19 0,409 0,371 0,038 0,015 29 27,89 0,659 0,656 0,003 0,019 40 28,99 0,909 0,939 0,030 0,053 8 26,60 0,182 0,170 0,012 0,011 19 27,21 0,432 0,378 0,054 0,031 30 27,89 0,682 0,656 0,026 0,003 41 28,99 0,932 0,940 0,008 0,031 9 26,66 0,204 0,188 0,016 0,006 20 27,28 0,455 0,406 0,049 0,026 31 28,01 0,705 0,701 0,003 0,020 42 29,03 0,955 0,944 0,010 0,013 10 26,69 0,227 0,196 0,032 0,009 21 27,30 0,477 0,412 0,065 0,042 32 28,10 0,727 0,731 0,004 0,027 43 29,12 0,977 0,954 0,024 0,001 11 26,74 0,250 0,211 0,039 0,016 22 27,38 0,500 0,447 0,053 0,031 33 28,11 0,750 0,735 0,015 0,008 44 29,28 1,000 0,968 0,032 0,009 В данном примере максимальные значения dn+ и dn- одинаковы и равны 0,067. Из табл. распределения Колмогорова для a = 0,1 найдем значение la = 1,22 . Для n=44 критическое значение отклонения определится из соотношения d m (a = 0 ,1) n = la ; d m (0 ,1) = 1,22 = 0 ,184. 44 Поскольку величина max dn=0,067 меньше критического значения, гипотеза о принадлежности выборки нормальному закону не отвергается. Критерий Р. Мизеса В качестве меры различия теоретической функции распределения F(x) и эмпирической Fn(x) по критерию Мизеса (критерию - w 2 ) выступает средний квадрат отклонений по всем значениям аргумента x +¥ w n2 = ò [Fn ( x ) - F ( x )] dF ( x ) 2 (13) -¥ Статистика критерия 68 n é i - 0 ,5 ù 1 + å ê F ( xi ) nw = n úû 12n i =1 ë 2 2 n (14) При неограниченном увеличении n существует предельное распределение статистики nwn2 . Задав значение вероятности a можно определить критические значения nwn2 (a ) . Проверка гипотезы о законе распределения осуществляется обычным образом: если фактическое значение nwn2 окажется больше критического или равно ему, то согласно критерию Мизеса с уровнем значимости a гипотеза НО о том, что закон распределения генеральной совокупности соответствует F(x), должна быть отвергнута. Пример. Проверить с помощью критерия Мизеса гипотезу о том, что ЭД, представленные в предыдущем примере, подчиняются нормальному распределению при уровне значимости a = 0,1. Решение. Исходные данные и результаты вычислений представлены в табл. 4. Таблица 4 i xi 1 25,79 2 25,98 3 25,98 4 26,12 5 26,13 6 26,49 7 26,52 8 26,60 9 26,66 10 26,69 11 26,74 Fn(xi ) F(xi ) Di х10-3 i xi Fn(xi ) F(xi ) Di х10-3 i xi Fn(xi ) F(xi ) Di х10-3 i xi Fn(xi ) F(xi ) Di х10-3 0,011 0,036 0,618 12 26,85 0,261 0,246 0,231 23 27,40 0,511 0,456 3,103 34 28,37 0,761 0,817 3,090 0,034 0,055 0,429 13 26,90 0,284 0,263 0,439 24 27,49 0,534 0,492 1,765 35 28,38 0,784 0,819 1,230 0,057 0,055 0,003 14 26,91 0,307 0,267 1,572 25 27,64 0,557 0,555 0,003 36 28,50 0,807 0,851 1,908 0,080 0,073 0,047 15 26,96 0,330 0,284 2,071 26 27,66 0,580 0,561 0,332 37 28,63 0,830 0,879 2,461 0,102 0,075 0,726 16 27,02 0,352 0,305 2,243 27 27,71 0,602 0,583 0,374 38 28,67 0,852 0,888 1,271 0,125 0,144 0,378 17 27,11 0,375 0,337 1,467 28 27,78 0,625 0,610 0,216 39 28,90 0,875 0,928 2,791 0,148 0,151 0,009 18 27,19 0,398 0,371 0,717 29 27,89 0,648 0,656 0,063 40 28,99 0,898 0,939 1,737 0,171 0,170 0,000 19 27,21 0,421 0,378 1,790 30 27,89 0,671 0,656 0,213 41 28,99 0,921 0,940 0,381 0,193 0,188 0,025 20 27,28 0,443 0,406 1,391 31 28,01 0,693 0,701 0,067 42 29,03 0,943 0,944 0,001 0,216 0,196 0,409 21 27,30 0,466 0,412 2,866 32 28,10 0,716 0,731 0,238 43 29,12 0,966 0,954 0,149 0,237 0,211 0,742 22 27,38 0,489 0,447 1,755 33 28,11 0,739 0,735 0,013 44 29,28 0,989 0,968 0,432 В этой таблице: Fn(xi)=(i–0,5)/44 – значение эмпирической функции распределения; F(xi) – значение теоретической функции распределения, соответствует значению функции нормального распределения в точке xi; 2 Di = [Fn ( xi ) - F ( xi )] . Критическое значение статистики критерия Мизеса при заданном уровне значимости равно 0,347 (таблица распределения Мизеса).Фактическое значение статистики 44 1 nwn2 = + å Di = 0 ,044 , 12 × 44 i=1 69 что меньше критического значения. Следовательно, гипотеза Н0 не противоречит имеющимся данным. Достоинством критерия Мизеса является быстрая сходимость к предельному закону, для этого достаточно не менее 40 наблюдений в области часто используемых на практике больших значений nwn2 (а не несколько сот, как для критерия хиквадрат). Сопоставляя возможности различных критериев, необходимо отметить следующие особенности: · Критерий Пирсона устойчив к отдельным случайным ошибкам в ЭД. Однако его применение требует группирования данных по интервалам, выбор которых относительно произволен и подвержен противоречивым рекомендациям. · Критерий Колмогорова слабо чувствителен к виду закона распределения и подвержен влиянию помех в исходной выборке, но прост в применении. · Критерий Мизеса имеет ряд общих свойств с критерием Колмогорова: оба основаны непосредственно на результатах наблюдения и не требуют построения статистического ряда, что повышает объективность выводов; оба не учитывают уменьшение числа степеней свободы при определении параметров распределения по выборке, а это ведет к риску принятия ошибочной гипотезы. Их предпочтительно применять в тех случаях, когда параметры закона распределения известны априори, например, при проверке датчиков случайных чисел. При проверке гипотез о законе распределения следует помнить, что слишком хорошее совпадение с выбранным законом распределения может быть обусловлено некачественным экспериментом («подчистка» ЭД) или предвзятой предварительной обработкой результатов (некоторые результаты отбрасываются или округляются). Выбор критерия проверки гипотезы относительно произволен. Разные критерии могут давать различные выводы о справедливости гипотезы, окончательное заключение в таком случае принимается на основе неформальных соображений. Точно также нет однозначных рекомендаций по выбору уровня значимости. Рассмотренный подход к проверке гипотез, основанный на применении специальных таблиц критических точек распределения, сложился в эпоху "ручной" обработки ЭД, когда наличие таких таблиц существенно снижало трудоемкость вычислений. В настоящее время математические пакеты включают процедуры вычисления стандартных функций распределений, что позволяет отказаться от использования таблиц, но может потребовать изменения правил проверки. Лекция 7 ПРОВЕРКА СТАТИСТИЧЕСКИХ ГИПОТЕЗ 1. Сущность задачи проверки статистических гипотез Статистическая гипотеза представляет собой некоторое предположение о законе распределения случайной величины или о параметрах этого закона, формулируемое на основе выборки. 70 Примерами статистических гипотез являются предположения: генеральная совокупность распределена по экспоненциальному закону; математические ожидания двух экспоненциально распределенных выборок равны друг другу. В первой из них высказано предположение о виде закона распределения, а во второй – о параметрах двух распределений. Гипотезы, в основе которых нет никаких допущений о конкретном виде закона распределения, называют непараметрическими, в противном случае – параметрическими. Гипотезу, утверждающую, что различие между сравниваемыми характеристиками отсутствует, а наблюдаемые отклонения объясняются лишь случайными колебаниями в выборках, на основании которых производится сравнение, называют нулевой (основной) гипотезой и обозначают Н0. Наряду с основной гипотезой рассматривают и альтернативную (конкурирующую, противоречащую) ей гипотезу Н1. И если нулевая гипотеза будет отвергнута, то будет иметь место альтернативная гипотеза. Различают простые и сложные гипотезы. Гипотезу называют простой, если она однозначно характеризует параметр распределения случайной величины. Например, если l является параметром экспоненциального распределения, то гипотеза Н0 о равенстве l=10 – простая гипотеза. Сложной называют гипотезу, которая состоит из конечного или бесконечного множества простых гипотез. Сложная гипотеза Н0 о неравенстве l>10 состоит из бесконечного множества простых гипотез Н0 о равенстве l=bi, где bi – любое число, большее 10. Гипотеза Н0 о том, что математическое ожидание нормального распределения равно двум при неизвестной дисперсии, тоже является сложной. Сложной гипотезой будет предположение о распределении случайной величины Х по нормальному закону, если не фиксируются конкретные значения математического ожидания и дисперсии. Проверка гипотезы основывается на вычислении некоторой случайной величины – критерия, точное или приближенное распределение которого известно. Обозначим эту величину через z, ее значение является функцией от элементов выборки z=z(x1, x2, …, xn). Процедура проверки гипотезы предписывает каждому значению критерия одно из двух решений – принять или отвергнуть гипотезу. Тем самым все выборочное пространство и соответственно множество значений критерия делятся на два непересекающихся подмножества S0 и S1. Если значение критерия z попадает в область S0, то гипотеза принимается, а если в область S1, то гипотеза отклоняется. Множество S0 называется областью принятия гипотезы или областью допустимых значений, а множество S1 – областью отклонения гипотезы или критической областью. Выбор одной области однозначно определяет и другую область. Принятие или отклонение гипотезы Н0 по случайной выборке соответствует истине с некоторой вероятностью и, соответственно, возможны два рода ошибок. Ошибка первого рода возникает с вероятностью a тогда, когда отвергается верная гипотеза Н0 и принимается конкурирующая гипотеза Н1. Ошибка второго рода возникает с вероятностью b в том случае, когда принимается неверная гипотеза Н0, в то время как справедлива конкурирующая гипотеза Н1. 71 Доверительная вероятность – это вероятность не совершить ошибку первого рода и принять верную гипотезу Н0. Вероятность отвергнуть ложную гипотезу Н0 называется мощностью критерия. Следовательно, при проверке гипотезы возможны четыре варианта исходов, табл. 1. Таблица 1. Гипотеза Н0 Верна Неверна Решение Вероятность Примечание Принимается 1- a Доверительная вероятность Отвергается a Вероятность ошибки первого рода Принимается b Отвергается 1-b Вероятность ошибки второго рода Мощность критерия Например, рассмотрим случай, когда некоторая несмещенная оценка параметра q вычислена по выборке объема n, и эта оценка имеет плотность распределения f( q ), рис. 1. Рис. 1. Области принятия и отклонения гипотезы Предположим, что истинное значение оцениваемого параметра равно Т. Если рассматривать гипотезу Н0 о равенстве q =Т, то насколько велико должно быть различие между q и Т, чтобы эту гипотезу отвергнуть. Ответить на данный вопрос можно в статистическом смысле, рассматривая вероятность достижения некоторой заданной разности между q и Т на основе выборочного распределения параметра q . Целесообразно полагать одинаковыми значения вероятности выхода параметра q за нижний и верхний пределы интервала. Такое допущение во многих случаях позволяет минимизировать доверительный интервал, т.е. повысить мощность критерия проверки. Суммарная вероятность выхода параметра q за пределы интервала с границами q 1–a/2 и q a/2, составляет величину a . Эту величину следует выбирать такой, чтобы выход за пределы интервала был маловероятен. Если оценка параметра попала в заданный интервал, то нет основания подвергать сомнению проверяемую гипотезу. Гипотезу равенства q =Т можно принять. Но если после получения выборки окажется, что оценка выходит за установленные пределы, то в этом случае есть 72 серьезные основания отвергнуть гипотезу Н0. Отсюда следует, что вероятность допустить ошибку первого рода равна a (равна уровню значимости критерия). Если предположить, например, что истинное значение параметра в действительности равно Т+d , то согласно гипотезе Н0 о равенстве q =Т – вероятность того, что оценка параметра q попадет в область принятия гипотезы, составит b , рис. 1. При заданном объеме выборки вероятность совершения ошибки первого рода можно уменьшить, снижая уровень значимости a . Однако при этом увеличивается вероятность ошибки второго рода b (снижается мощность критерия). Аналогичные рассуждения можно провести для случая, когда истинное значение параметра равно Т–d. Единственный способ уменьшить обе вероятности состоит в увеличении объема выборки (плотность распределения оценки параметра при этом становится более "узкой"). При выборе критической области руководствуются правилом Неймана – Пирсона: следует так выбирать критическую область, чтобы вероятность a была мала, если гипотеза верна, и велика в противном случае. Однако выбор конкретного значения a относительно произволен. Употребительные значения лежат в пределах от 0,001 до 0,2. В целях упрощения ручных расчетов составлены таблицы интервалов с границами q 1–a/2 и q a/2 для типовых значений a и различных способов построения критерия. При выборе уровня значимости необходимо учитывать мощность критерия при альтернативной гипотезе. Иногда большая мощность критерия оказывается существеннее малого уровня значимости, и его значение выбирают относительно большим, например 0,2. Такой выбор оправдан, если последствия ошибок второго рода более существенны, чем ошибок первого рода. Например, если отвергнуто правильное решение "продолжить работу пользователей с текущими паролями", то ошибка первого рода приведет к некоторой задержке в нормальном функционировании системы, связанной со сменой паролей. Если же принято решения не менять пароли, несмотря на опасность несанкционированного доступа посторонних лиц к информации, то эта ошибка повлечет более серьезные последствия. В зависимости от сущности проверяемой гипотезы и используемых мер расхождения оценки характеристики от ее теоретического значения применяют различные критерии. К числу наиболее часто применяемых критериев для проверки гипотез о законах распределения относят критерии хи-квадрат Пирсона, Колмогорова, Мизеса, Вилкоксона, о значениях параметров – критерии Фишера, Стьюдента. Лекция 9 МЕТОДЫ ОЦЕНКИ ПАРАМЕТРОВ РАСПРЕДЕЛЕНИЯ 1. Точечная оценка параметров распределения Сущность задачи точечного оценивания параметров 73 Точечная оценка предполагает нахождение единственной числовой величины, которая и принимается за значение параметра. Такую оценку целесообразно определять в тех случаях, когда объем ЭД достаточно велик. Причем не существует единого понятия о достаточном объеме ЭД, его значение зависит от вида оцениваемого параметра (к этому вопросу вернёмся при изучении методов интервальной оценки параметров, а предварительно будем считать достаточной выборку, содержащую не менее чем 10 значений). При малом объеме ЭД точечные оценки могут значительно отличаться от истинных значений параметров, что делает их непригодными для использования. Задача точечной оценки параметров в типовом варианте постановки состоит в следующем: Имеется выборка наблюдений (x1, x2, …, xп) за случайной величиной Х. Выборка должна быть представительной. Объем выборки n фиксирован. Известен вид закона распределения величины Х, например, в форме плотности распределения f(T, x), где T – неизвестный (в общем случае векторный) параметр распределения. Параметр является неслучайной величиной. Требуется найти оценку q параметра T закона распределения. Существует несколько методов решения задачи точечной оценки параметров, самыми распространёнными из которых являются методы максимального (наибольшего) правдоподобия, моментов и квантилей. Метод максимального правдоподобия Метод был предложен Р. Фишером в 1912 г. Метод основан на исследовании вероятности получения выборки наблюдений ( x1 , x2 ,....xn ) за некоторым количественным признаком Х. Значения выборки случайные, поэтому будем рассматривать их как независимые случайные величины X 1 , X 2 ,...X n , имеющие одинаковые распределения. Если количественный признак Х непрерывный, то случайные величины X 1 , X 2 ,...X n имеют одинаковые плотности распределения вероятностей, зависящих от х и от вектора параметров распределения q = (q1 ,q 2 ,...q n ) , т.е. плотности f (x ,q ). Суть метода максимального правдоподобия заключается в свойстве случайной величины реализовывать в эксперименте в основном те свои значения x1 , x2 ,....xn , вероятность которых максимальна. В этом случае совместная плотность распределения вероятностей f (x1 , x2 ,...xn ,q ) случайных величин X 1 , X 2 ,...X n должна быть максимальной для непрерывного признака или совместное распределение вероятностей P ( X 1 = x1 , X 2 = x 2 ,... X n = x n ,) случайных величин X 1 , X 2 ,...X n должно быть максимальным для дискретного признака. Ввиду независимости случайных величин X 1 , X 2 ,...X n их совместная плотность распределения будет равна произведению плотностей распределения, т.е. f (x1 , x2 ,...xn ,q ) = f (x1 ,q ) × f (x2 ,q ) × ... × f (xn ,q ), 74 а совместная вероятность равна произведению вероятностей событий P ( X 1 = x1 , X 2 = x2 ,...X n = xn ,q ) = p x1 (q ) × ...× p xn (q ). Функцией правдоподобия L(x1 , x2 ,...xn ,q ) называется функция, которая в случае непрерывного количественного признака определяется по формулам: для непрерывного признака L(x1 , x2 ,...xn ,q ) = f (x1 ,q ) × f (x2 ,q ) × ... × f (xn ,q ). для дискретного признака L (x1 , x 2 ,... x n ,q ) = p x1 (q )× ... × p x n (q ). ~ ~ В качестве оценки q 1 ,.... q n неизвестных параметров распределения q1 ,...q n берут те значения, при которых функция правдоподобия достигает максимума. Технически задача поиска максимального значения функции правдоподобия облегчается, если рассмотреть не саму функцию, а натуральный логарифм от неё, т.е. функцию ln L(x1 , x2 ,...xn ,q ) . За оценки неизвестного вектора параметров q берётся ~ ~ решение q1 ,....q n уравнения правдоподобия ¶ ln L(x1 , x2 ,...xn ,q ) = 0 ; k = 1, n . ¶q k Для проверки того, что точка оптимума соответствует максимуму функции правдоподобия, необходимо найти вторую производную от этой функции. И если вторая производная в точке оптимума отрицательна, то найденные значения параметров максимизируют функцию. Итак, нахождение оценок максимального правдоподобия включает следующие этапы: построение функции правдоподобия (ее натурального логарифма); дифференцирование функции по искомым параметрам и составление системы уравнений; решение системы уравнений для нахождения оценок; определение второй производной функции, проверку ее знака в точке оптимума первой производной и формирование выводов. Пример1. Время до выхода из строя прибора определяется показательным распределением ìl e - l x , x ³ 0 ; f (x ) = í x < 0. î 0, По выборке 2, 3, 5, 6, 8 найти оценку параметра l . 75 Решение: Функция правдоподобия имеет вид: L(2 ,3,5,6 ,8,l ) = le - l 2 × le - l 3 × le - l 5 × le - l 6 × le - l 8 = l5e - l 24 . Тогда ln L(2 ,3,5,6 ,8,l ) = 5 ln l - 24l . Уравнение правдоподобия будет иметь вид: ¶ ln L 5 ~ 5 = - 24 = 0 . Отсюда получаем оценку l = . ¶l l 24 Метод максимального правдоподобия позволяет получить состоятельные, эффективные (если таковые существуют, то полученное решение даст эффективные оценки), достаточные, асимптотически нормально распределенные оценки. Этот метод может давать как смещенные, так и несмещенные оценки. Смещение удается устранить введением поправок. Метод особенно полезен при малых выборках. Если функция максимального правдоподобия имеет несколько максимумов, то из них выбирают глобальный. Метод моментов Метод предложен К. Пирсоном в 1894 г. Будем считать, что вид функции распределения изучаемого количественного признака известен, но параметры этого распределения неизвестны. Нужно оценить эти параметры по выборке. Сущность метода моментов заключается в том, что по конкретной выборке x1 , x2 ,....xn генеральной совокупности Х, распределение которой известно с точностью до параметров q1 ,...q m выбирается столько эмпирических моментов, сколько требуется оценить неизвестных параметров распределения. Желательно применять моменты младших порядков, так как погрешности вычисления оценок резко возрастают с увеличением порядка момента. Вычисленные по ЭД оценки моментов приравниваются к теоретическим моментам; параметры распределения определяются через моменты, и составляются уравнения, выражающие зависимость параметров от моментов, в результате получается система уравнений. Решение этой системы дает оценки параметров распределения генеральной совокупности. Пример 2. Время до выхода из строя прибора определяется показательным распределением ìl e - l x , x ³ 0 ; f (x ) = í x < 0. î 0, 76 По выборке x1 , x2 ,....xn найдём точечную оценку неизвестного параметра распределения l . Так как для показательного распределения нужно оценить единственный параметр l , то для его оценивания нужно одно уравнение. Приравняем, например, первые начальные моменты – теоретический и выборочный. Первый начальный теоретический момент равен ¥ m1 = ò x × le -lx dx = 1 . l Первый начальный выборочный момент вычисляется по формуле xi . i =1 n n m1B = å Приравнивая их получим 1 1 n xi ~ = å . Отсюда l = n . x l i=1 n å i i =1 n 2. Интервальная оценка параметров распределения Сущность задачи интервального оценивания параметров Интервальный метод оценивания параметров распределения случайных величин заключается в определении интервала (а не единичного значения), в котором с заданной степенью достоверности будет заключено значение оцениваемого параметра. Интервальная оценка характеризуется двумя числами – концами интервала, внутри которого предположительно находится истинное значение параметра. Иначе говоря, вместо отдельной точки для оцениваемого параметра можно установить интервал значений, одна из точек которого является своего рода "лучшей" оценкой. Интервальные оценки являются более полными и надежными по сравнению с точечными, они применяются как для больших, так и для малых выборок. Совокупность методов определения промежутка, в котором лежит значение параметра Т, получила название методов интервального оценивания. К их числу принадлежит метод Неймана. Постановка задачи интервальной оценки параметров заключается в следующем: Имеется выборка наблюдений x1 , x2 ,....xn за случайной величиной Х. Объем выборки ~ n фиксирован. Предположим, что статистическая характеристика q , рассчитанная по данным выборки, является выборочной оценкой неизвестного параметра q гене- 77 ~ ральной совокупности, причём q - это постоянное число. Оценка q характеризует ~ параметр q тем точнее, чем меньше абсолютная величина разности q - q . ~ ~ Если q - q < d , где d > 0 , то число d называется точностью оценки q . ~ Нельзя утверждать, что оценка q абсолютно удовлетворяет неравенству ~ q - q < d . Однако можно задать вероятность g , с которой это неравенство осуществляется. ~ Надёжность или доверительная вероятность оценки q - это вероят~ ность g , с которой осуществляется неравенство q - q < d . Принято надёжность оценки задавать перед процессом оценивания параметра генеральной совокупности. В качестве вероятности g берут число, близкое к 1 (от 0,95 до 0,999). ( ) ~ ~ Пусть P q - q < d = g . Если неравенство q - q < d заменить равносильным ему двойным неравенством: ~ ~ ~ - d < q - q < d , или q - d < q < q + d , ( ) ~ ~ то получим: P q - d < q + d = g . Вероятностный смысл данного отношения таков: вероятность того, что интервал ~ ~ ~ q - d ;q + d покрывает неизвестный параметр q , равна g . ( ) ( ) ~ ~ Интервал q - d ;q + d , который покрывает оцениваемый неизвестный параметр ~ q , с заданной вероятностью g называется доверительным интервалом. ~ ~ Границы q - d и q + d называются доверительными границами интервала. Они определяются на основе выборочных данных и являются функциями от случайных величин X 1 , X 2 ,...X n , и, следовательно, сами являются случайными величинами. Общий метод построения доверительных интервалов Доверительный интервал для оценки математического ожидания нормального распределения при известном среднем квадратичном отклонении 78 Будем считать, что количественный признак Х генеральной совокупности подчиняется нормальному закону распределения с параметрами a и s . Среднее квадратичное данного распределения s считается известным. Нужно оценить неизвестное математическое ожидание a генеральной совокупности по выборочной средней x , т.е. необходимо определить доверительный интервал, который покрывает параметр a с заданной надёжностью g . Имеет место следующее утверждение: если случайная величина Х подчиняется нормальному закону распределения, то выборочная средняя X , найденная по независимым наблюдениям, также подчиняется нормальному закону распределения с параметрами: M (X ) = a ; s ( X ) = s . n Для того, чтобы определить доверительный интервал, необходимо, чтобы выполнялось равенство P ( X - a < d ) = g . Для нормально распределённой случайной величины Х имеет место равенство: æ d P ( X - a < d ) = 2Фçç ès X ( ) X имеет нормальное распределение и s ( X ) = ö ÷÷. ø s . Поэтому n æd n ö d n ÷÷ = 2Ф(t ) ; t = P ( X - a < d ) = 2Фçç . s s è ø Из равенства t = d n s выразим d = t и получим: s n s ö æ Pç X - a < t ÷ = 2Ф(t ). nø è Так как вероятность Р задана и равна g , окончательно получаем: s s ö æ ln A , то q = q 2 с вероятностью ошибки b ; если же L t < ln B , то q = q1 с вероятностью ошибки a. Принято рассматривать в описанной ситуации две величины - n и Dgt , определяемые так: то минимальное значение t, при котором выполняется одно из неравенств L t > ln A или L t < ln B - является случайной величиной, зависящей от той или иной реализации случайной величины x 1 , x 2 , K, x k ,K ; именно это минимальное значение t и является случайной величиной n ; учитывая смысл только что введенного обозначения n , введем величину f (x i , q 2 ) 1 n Dgt = å ln ; n i =1 f (xi ,q1 ) очевидно, и эта величина - также случайная. На практике часто имеют смысл математические ожидания последних двух случайных величин; их находят как обычные средние арифметические значения реализаций этих величин, конструируемых по различным реализациям x 1 , x 2 , K, x k ,K Лекция 10 Интервальные оценки при экспоненциальном законе распределения В качестве примера рассмотрим определение интервальных оценок для показателей надёжности. Вероятность безотказной работы определяется формулой 1 P(t ) = e -lt , l = , T = mt . T Доверительный интервал для средней наработки до отказа Т при равных вероa ятностях выхода за правую и левую границы 2 2Tå r 2Tå r < T < , c a2 2 c12-a 2 где ca2 2 , c12-a 2 находятся по таблице при числе степеней свободы k = 2r , где r число зарегистрированных отказов, для соответствующих вероятностей. Когда вычисляется только нижняя граница, 2T r T > å2 . ca Если k > 100 , то доверительный интервал вычисляется по формулам для нормального усечённого распределения. Для вероятности безотказной работы двухсторонний доверительный интервал вычисляется по формуле 81 é ca2 2ti ù é c12-a 2ti ù exp êú < P(ti ) < exp êú. T T 2 2 êë êë år ú år ú û û В случае односторонней границы é ca2 ti ù P(ti ) > exp êú. T 2 ú år û ëê Лекция 10 10.1 Планирование испытаний и обработка экспериментальных данных В соответствии с требованиями ГОСТ 27.002-83 планирование испытаний предусматривает ряд предварительных условий, обеспечивающих эффективность испытаний. Вводятся условные обозначения различных планов в виде совокупности трех символов, первый из которых указывает число испытываемых объектов (устройств) N, второй - наличие (R) или отсутствие (U) замены (восстановления) объектов, отказавших во время испытаний, третий - длительность испытаний (r или Т). Таким образом, для испытаний N объектов без замены отказавших элементов, имеем следующие три плана: (N, U, r) - испытания до r-го отказа, r< N; (N, U, T) - испытания длительностью Т; [N, U, (r, T)] - испытание длительностью, равной min (t r или T ) , где t r - момент r-го отказа, а Т - заведомо заданный ресурс (время, или величина пробега, или число циклов и т.д.) Аналогично вводятся обозначения для планов с заменой (восстановлением) отказавших устройств: (N, R, r); (N, R, T); [N, R, (r, T)]. В плане (N, R, r) в отличие от (N, U, r) число r может быть больше, чем N (где, в частности, допустимо N = 1). Здесь, приведено 6 наиболее распространенных типов испытаний. ГОСТ 27.001-83 предусматривает 16 планов испытаний, где учтены кроме названных условий и такие как M - восстановление объектов при испытаниях в случае их отказов; S - решение об окончании испытаний (о приемке или браковке) восстанавливаемых объектов (основывается на суммарном времени испытаний). Таблица 1. План испытаний Суммарная наработка Tå r , ч r (NUr) Tå r = å t j + ( N - r ) × t r (NUT) Tå r = å t j + ( N - r ) × T NU(rT) при t r < T , Tå r = å t j + ( N - r ) × t r ; j =1 r j =1 r j =1 82 r при t r ³ T , Tå r = å t j + ( N - r ) × T j =1 (NRr) Tå r = Ntr (NRT) Tå r = NT NR(rT) при t r < T , Tå r = Nt r ; при t r ³ T , Tå r = NT Результаты статистической обработки испытаний существенно зависят от вероятностных моделей, то есть от априорных (теоретических) распределений интервалов безотказной работы и восстановлений. Эти результаты могут приводить к заведомо ошибочным выводам, если модель не отражает реальные процессы возникновения отказов и механизмы восстановления. Поэтому до решения основных задач апостериорного (на основе опыта) анализа надежности целесообразно сначала проверить, с помощью статистического критерия согласия, на соответствие выбранного априорного распределения эмпирическому распределению, построенному на основании данных проведенных испытаний. Исходными данными (случайными величинами), которые подвергаются обработке, являются время наработки на отказ, время наработки на восстановление и число отказов однотипных элементов. После того, как такой материал собран, его обработка позволяет установить законы распределения показателей надежности: вероятность безотказной работы, интенсивность отказов, среднее время наработки на отказ и др. Знание законов распределения дает возможность определить все остальные количественные показатели надежности. Таким образом, основная задача статистической обработки состоит в определении одного из законов распределения исходных экспериментальных величин. В ряде случаев вид закона распределения известен заранее, до опыта. Например, как уже отмечалось выше, для электронной аппаратуры средств автоматики и релейной защиты справедлив экспоненциальный закон распределения показателей надежности. Это подтверждается многочисленными опытными данными, полученными в условиях эксплуатации. При определении или подтверждении закона распределения целесообразен следующий порядок: подготовка опытных данных; построение гистограмм оцениваемого количественного показателя надежности; аппроксимация гистограмм теоретическим законом распределения и определение его параметров; проверка допустимости предполагаемого закона распределения на основе ис2 пользования критериев согласия. Наиболее часто используется критерий c или критерий Колмогорова. Для получения достаточно точных результатов число наблюдений случайной величины (отказов) должно быть не менее 40-50. 10.2. Интервальная оценка показателей надежности Количество статистических данных для оценки надежности, полученных в процессе эксплуатации, принципиально ограничено. Полученные по ограниченному объему информации точечные оценки могут оказаться весьма приблизительными. Причем отклонения этих оценок от истинного значе83 ния оцениваемого параметра являются величинами случайными. Очевидно, что с увеличением числа наблюдений (отказов) случайная ошибка оценки показателей уменьшается. На основе опытных данных используется специальная методика оценки показателей надежности в определенном интервале возможных их значений. Предположим, что истинное значение средней наработки до отказа составляет Т0, а средняя наработка до отказа определена по полученным отказам: T̂ = 1 n å ti , n i=1 где n - количество отказов за время испытаний, ti - наработка до i-го отказа. Чем меньше n тем больше расхождение между T0 и T̂ , то есть существует интервал расхождения. Найти точные границы, в пределах которых находится истинное значение искомой величины, не представляется возможным. Однако можно определить интервал ее возможных значений с некоторой доверительной вероятностью b . При этом, чем больше доверительная вероятность b , тем шире границы интервала и наоборот. В общем виде, эта зависимость имеет запись b = P(TH £ T0 £ TB ) (1) где Тн и Тв - соответственно нижняя и верхняя границы средней наработки до отказа, где лежат T̂ и Т0 . Вероятность того, что значение Т0 выйдет за заданный интервал называется уровнем значимости: a = P(TH > T0 > TB ) = 1 - b . (2) Значения доверительных вероятностей b обычно принимают равными 0,9; 0,95; 0,99. Соответствующие им уровни значимости составят 0,1; 0,05; 0,01. Доверительная вероятность b , определяемая выражением (1), характеризует степень достоверности результатов двусторонней (то есть с определением верхней и нижней границ) оценки. Доверительный интервал для средней наработки до отказа при равных вероятностях a выхода за 2 правую (верхнюю) и левую (нижнюю) границы для экспоненциального распределения определяется по выражению TH = 2Tå r c a2 2 где c a2 2 ;2r и c 12-a 2 ;2r - значения < T̂ < ,2 r 2Tå r c2a = TB , 1- ,2 r 2 c 2 (хи-квадрат) при параметрах степеней свободы, для вероятностей P = (3) a a и 1- ; 2 2 2r = k - число a a ; P = 1 - соответственно. 2 2 Когда вычисляется только нижняя граница, то 84 T̂ > 2Tå r (4) ca2 ,2 r В выражениях (3) и (4) Tå r - суммарная наработка до отказа по отказам, зафиксированным во время эксперимента. Значения c a2 2 деления и ;2r c 12-a 2 ;2r определяются по таблице П-1 квантилей распре- c 2 (хи-квадрат). Таким образом, для заданных уровней значимости a и числа степеней свободы k по таблице на2 ходят соответствующие значения c , подставляют в выражение (3) и находят Tн и Tв. Величина a задается в зависимости от требований, предъявляемых к анализируемой системе. 10.3. Статистические оценки генеральной совокупности. Задача об оценке качества по выборке Пусть x - случайная величина и x 1 ,...,x n - ее значения, совокупность ко-торых принято называть генеральной совокупностью. Напомним две основные статистические оценки генеральной совокупности: математическое ожидание - 1 n 2 Dx = å (x i - M x ) . n i= 1 Часто вместо символов число Mx = 1 n åxi , n i= 1 и дисперсия – число M x ,Dx пишут M(x), D(x); кроме того, во многих случаях применяется обозначение s(x) для числа D(x ) . Если наряду со случайной величиной x имеется вторая случайная величина h и набор ее значений того же объема - h1 ,..., h n , то рассматривают кор-реляционный момент пары случайных величин x,h - число kxh = 1 n å (x i - M(x ))(h i - M (h )) , n i =1 а также коэффициент корреляции rxh = kxh s (x )s (h ) . Если о паре случайных величин x,h известно, что между ними существует некоторая зависимость h = h x , то имеется методика приближенного отыскания этой зависимости, основанная () на методе наименьших квадратов. Если требуется среди многочленов степени, не превосходящей некоторого m, найти тот, который имитирует эту зависимость лучше всего, то достаточно решить несложную систему линейных алгебраических уравнений и при m < n получить определенный m многочлен Q( x ) = q0 + q1x +...+ qm x ; доказывается, что при наличии зависимости h = h(x ) дисперсия D(h) не превосходит числа 1 n 2 å (h i - Q(x i )) , n - m i=1 которое принято называть оценкой дисперсии по методу наименьших квадра-тов. 85 Рассмотрим классическую задачу о статистической оценке качества. Имеется партия некоторого товара (например, партия автомобилей). Требуется оценить качество этой партии, сделав некоторую выборку из нее (например, отобрав несколько автомобилей, оценить их и по этим оценкам сделать выводы обо всей партии). В такой постановке речь идет об оценке качества генеральной совокупности по выборке. Излагаемая ниже методика решения этой задачи предложена А. Вальдом. Введем какой-либо числовой параметр, характеризующий качество партии. Например, если речь идет о партии автомобилей, то таким параметром может быть доля (процент) автомобилей бракованных. Если этот параметр меньше некоторого числа, то партия будет считаться хорошей, а если этот параметр больше некоторого другого числа, то партия будет считаться плохой. Пусть a - вероятность забраковать хорошую партию (т.е. вероятность того, что анализ выборки покажет, что партия плохая в то время, как на самом деле, партия хорошая). Пусть b - вероятность принять плохую партию (т.е. вероятность того, что анализ выборки покажет, что партия хорошая в то время, как в действительности - партия плохая). Пусть N - объем партии и n0 - объем выборки. Пусть, далее, n - приемочное число, а т - количество дефектных изделий в выборке, т.е. если в выборке m £ n , то партия принимается (или партия - хорошая), а если в выборке т> n , то вся партия не принимается (или партия - плохая). Наконец, обозначим через L количество дефектных изделий во всей партии и пусть l0 ,l1 , l0 < l1 - границы контролируемого параметра L, т.е. если L £ l0 , то партия - хорошая, а если L ³ l1 , то партия - плохая. В этих условиях вероятности a, b можно вычислить явно с помощью комбинаторных формул. Напишем лишь окончательный результат, используя стандартное обозначение P(a) для вероятности наступления события a: n Clm0 C Nn0--l0m m=0 C Nn0 a = P( M > n | L = l0 ) = 1 - å n Cl1m C Nn0--l1m m =0 C Nn0 b = P( M £ n | L = l1 ) = å Если n0 £ N , то a, b приближенно равны: 10 n n , . a @ 1 - å Cnm0 p0m (1 - p0 ) 0 , b @ å Cnm0 p1m (1 - p1 ) 0 , m=0 где p0 = l0 , N p1 = n -m n -m m=0 l1 . Если p0 , p1 < 01 . , то возможны дальнейшие упроще- ния формул. N Числа a, b называют, соответственно, риском поставщика и риском потребителя. При составлении плана контроля партии по выборке стремятся так выбрать границы l0 < l1 , чтобы эти риски были, по возможности, одинаковыми. Примечание 1. Практическое вычисление выражений типа n C m C n0 - m n C m C n0 - m l0 N - l0 l N -l a = 1- å , b = å 1 n0 1 n0 CN CN m= 0 m= 0 требует предварительного преобразования их к виду, удобному для реализации вычислений на компьютере; без такой подготовки непосредственное вычисление факториалов создает ситуацию «переполнения» памяти. Вот основное в этой связи наблюдение. k Выражение Cn представляет собой следующее произведение n - k +1 n n -1 n - 2 × × ×K × . k k -1 k - 2 1 86 Здесь ровно k сомножителей. Теперь рассмотрим выражение Clm0 C Nn0--l0m C Nn0 ; распишем его в соответствии с указанным только что правилом: l0 l0 - 1 l - m + 1 N - l0 N - l0 - 1 N - l0 - n0 + m + 1 × ×K× 0 × × ×K× m m -1 n0 - m n0 - m - 1 1 1 . N N -1 N - n0 + 1 ×K × × n0 n0 - 1 1 В знаменателе этого выражения указано n0 дробей; в числителе количество дробей есть m + ( n0 - m ) , то есть – их столько же, сколько в знаменателе. Примечание 2. В практических задачах, связанных с статистической оценкой генеральной совокупности, потребуется следующая справочная информация (приводится здесь исключительно для удобства читателя). 1) Случайная величина с равномерным распределением в интервале (a,b) имеет плотность 0, x < a , ì ï 1 , a £ x £ b, f (( x ) = í b a ï 0 , x > b. î 2 a+b ( b - a) и ее дисперсия D = . Ее математическое ожидание M = 2 12 Для имитации на компьютере случайной величины, равномерно распределен-ной в интервале (0,1), в каждом языке программирования предусмотрена специальная команда; будем условно обозначать ее RND (то есть символ RND обозначает число из совокупности равномерно распределенных чисел в интервале (0,1)). Таким образом, для имитации чисел, равномерно распределенных в интервале (a,b), нужно применять команду ( b - a )RND . 2) Случайная величина с экспоненциальным распределением имеет плотность ì 0, x < 0, f ( x ) = í -lx île , x ³ 0. 1 1 Ее математическое ожидание M = , дисперсия D = 2 . Можно доказать, что имитация экспоl l ненциально распределенной случайной величины осуществляется по схеме 1 = x - ln( 1 - RND ) . l 3) Случайная величина с нормальным распределением имеет плотность æ ( x - a )2 ö 1 ÷. f( x)= expçç 2b ÷ø 2pb è Ее математическое ожидание M = a , дисперсия D = b . Можно доказать, что имитация нормально распределенной случайной величины осуществляется по схеме x = 12 b æç å RND - 6 ö÷ + a . è1 ø 4) Случайная величина с логнормальным распределением имеет плот87 æ (ln( x ) - a )2 ö 1 æ1 ö ÷÷. Ее математическое ожидание M = expç b + a ÷ , expçç ность f ( x ) = 2b x 2pb è2 ø è ø дисперсия D = (exp( b + 2a ))(exp( b ) - 1) . Для имитации логнормально-распределенной случайной величины с параметрами a ,b , указанными в вышеприведенной формуле для плотности, следует имитировать нормально распределенную случайную величину с математическим ожиданием a и дисперсией b, а затем взять экспоненциальную функцию от каждого из этих чисел. 10.4. Случайные процессы как модели порождения данных. Линейные модели случайных процессов. Существует много задач, данные для которых возникают случайно. Например, задача о прогнозировании землетрясений опирается на данные, собираемые специальными датчиками и рассматриваемые современной наукой как случайные числа. При исследовании случайных данных всегда возникает необходимость в их имитации на уровне уже распознанных закономерностей тогда исследования можно проводить без дальнейшего сбора информации на уровне датчиков. В этой связи к настоящему времени сложилось несколько моделей порождения данных. Одной из таких моделей является случайный процесс. Случайным процессом называется функция вещественной переменной, значениями которой являются случайные величины. Если z( t ) - случайный процесс, то аргумент t удобно рассматривать как время. Мы будем предполагать, что для всех рассматриваемых ниже случайных процессов z ( t ) при t<0 обязательно z( t ) =0. Поскольку для каждого t величина z( t ) - случайная, постольку для каждого t имеется функция (или плотность) распределения этой случайной величины. В теории случайных процессов принято особым образом выделять такой случайный процесс z ( t ) , для которого каждая случайная величина z ( t ) является нормально распределенной с математическим ожиданием 0 и дисперсией 1 (для всех t³0); этот случайный процесс называется белым шумом и всюду ниже будет обозначаться символом a(t). По традиции, аргумент случайного процесса записывают в виде индекса, т.е. вместо z( t ) пишут zt , вместо a(t) пишут at и т.д. Условимся, что всюду в дальнейшем время будет меняться дискретно, только по целым числам, т.е. t = 0,±1,±2,... . Среди случайных процессов особое место занимают процессы линейные, которые определяются так: случайный процесс zt называется линейным, если имеет место равенство: zt = at + p 1at -1 + p 2 at -2 + ... + p s at -s + ..., где p 1 ,p 2 ,...,p s ,... - числовая последовательность, а суммирование в правой части последнего равенства происходит фактически на конечном числе слагаемых, так как, напомним, для всех рассматриваемых случайных процессов z( t ) при t<0 обязательно z ( t ) =0. При описании линейных случайных процессов удобно использовать технику формальных степенных рядов по следующей схеме. Будем обозначать символом B так называемый оператор сдвига, т.е. символ, позволяющий формально сделать запись: zt -1 = Bzt . Такая запись означает, что zt - 2 = Bzt -1 , а вместе эти две записи означают, что zt - 2 = Bzt -1 = B 2 zt . Продолжая этот формализм, получаем: zt - k = B k zt . Поэтому выражение zt = at + p 1at -1 + p 2 at - 2 +...+p sa t - s +... переписывается так: 88 zt = at + p 1Bat + p 2 B 2at +...+p k B k at +.... ; последнее выражение можно переписать в виде zt = p ( B )at , где p ( B ) = 1 + p 1 B + p 2 B 2 + ... + p k B k + ... - формальный степенной ряд. Учитывая, что у этого формального степенного ряда свободный член отличен от нуля, можно найти обратный к нему такой ряд y ( B) , что y ( B)p ( B) = 1 ; при этом y ( B) будет иметь следующий вид: y (B ) = 1 + y 1 B + y 2 B 2 + .. + y k B k + .. , где y 1 ,y 2 ,...,y k ,... - числовая последовательность. С помощью последнего формального степенного ряда линейный случайный процесс zt= at - y 1 zt -1 - y 2 z t - 2 -...-y k zt - k -... . переписывается в следующем виде: Таким образом, всякий линейный случайный процесс возникает по указанному выше правилу из белого шума и может задаваться через самого себя рекуррентно. Остановимся на вопросе компьютерной имитации линейного случайного процесса. Начнем с имитации белого шума. Из математики известно, что если взять двенадцать случайных чисел, распределенных равномерно в интервале (0,1) - s1 ,..., s12 - и с их помощью построить число s = æç å si ö÷ - 6 , то числа s , многократно повторенные с помощью такой процедуры, è i= 1 ø 12 будут нормально распределены с математическим ожиданием 0 и дисперсией 1. Отсюда ясно, как осуществить имитацию белого шума: если требуется воспроизвести значения a0 , a1 ,..., a n , то нужно построить с помощью датчика случайных чисел, равномерно распределенных в интервале (0,1), наборы si ,1 , si , 2 ,..., si ,12 , i = 0,1,..., n , а затем положить æ 12 ö ai = ç å si, k ÷ - 6, i = 0,1,..., n . è k= 1 ø Заметим, что если бы требовалось имитировать нормально распределенные случайные числа, но не с нулевым математическим ожиданием и единичной дисперсией, а с произвольно заданными математическим ожиданием m и дисперсией D, то несложное обобщение вышеупомянутой математической теоремы привело бы нас к следующей формуле: 12 ai = D æç æç å si ,k ö÷ - 6 ö÷ + m ,i = 0 ,1,...,n . è è k =1 ø ø В случае же произвольного линейного случайного процесса zt = at + p 1at -1 + p 2 at -2 + ... + p s at -s + ... последовательность действий для его имитации такова: 1) фиксируется число n, для которого будут имитироваться числа z0 , z1 ,..., z n ; 2) формируется имитация белого шума a0 , a1 ,...,an в соответствии с описанным выше правилом; 3) формируются числа z0 = a0 ; z1 = a1 + p 1a0 ; z 2 = a2 + p 1a1 + p 2 a0 ; ................................. z n = an + p 1an-1 + p 2 an-2 + ... + p n a0 . 89 90