Раздел 4.  Метод наименьших квадратов (продолжение)

На простейшем примере рассмотрим порядок решения задачи с использованием МНК.

При разработке нового образца газоразрядного источника света обнаружено, что основное влияние на его светоотдачу оказывают две компоненты наполнения, которые мы условно назовем добавка А и добавка В. Прежде чем приступить к постановке опытов по изучению влияния добавок, необходимо оценить дисперсию воспроизводимости результатов, т.е. решить вопрос - какой разброс значений параметра возможен не за счет влияния факторов, а за счет неконтролируемых причин. Поэтому изготовили пять опытных образцов с одинаковым наполнением и измерили их светоотдачу (в относительных единицах). Получены следующие результаты: y = 62.0, 63.6, 63.5, 61.0, 62.2. По этим данным оценена дисперсия S_у² = 1.2. Затем выполнена серия опытов при различных количествах добавок (в мг), результаты представлены в таблице 4.1.

Таблица 4.1. Светоотдача ламп у (в относительных единицах) в зависимости от количества добавок (количество добавки А в мг обозначено х₁, количество добавки В - х₂).

y = b₁ + b₂x₁ + b₃x₂ + b₄x₁x₂                          (4.8)

52.0 = b₁ + 40b₂ + 2b₃x₂ + 80b₄

58.1 = b₁ + 60b₂ + 2b₃x₂ + 120b₄

* * * * * * * * * * * *

74 = b₁ + 80b₂ + 6b₃x₂ + 480b₄

,

в ней 9 строк и 4 столбца. С такой матрицей решить задачу по формулам (4.4, 4.5) с помощью ЭВМ не представляет труда. Однако, сделав замену переменных, мы не только решим задачу "вручную", но продемонстрируем полезное свойство данного плана эксперимента. Перейдем к новым переменным по принципу , где - средний уровень данного фактора ( у нас для x₁ - это 60, а для x₂ - 4), а h _x - выбирается произвольно, но так, чтобы значения z были равны одной или нескольким единицам. Нам удобно взять h_x1 = 20, h_x2 = 2.   

Итак,
,              (4.9)

y = а₁ + а₂z₁ + a₃z₂ + a₄z₁z₂.                     (4.10)

          

и коэффициенты а определятся по формулам:

а₁ = (52 + 58.1 + 65.5 + 53.8 + 61.3 + 70 + 56.2 + 65.9 + 74) × 1 / 9 = 61.87

а₂ = (-52 + 65.5 - 53.8 + 70 - 56.2 + 74) × 1 / 6 = 7.92

а₃ = (-52 - 58.1 - 65.5 + 56.2 + 65.9 + 74) × 1 / 6 = 3.42

а₄ = (52 -65.5 - 56.2 + 74) × 1 / 4 = 1.075.

, ,   ,

t₁ = 169.4, t₂ = 17.7, t₃ = 7.6, t₄ = 1.96.

Задав уровень значимости, например α = 0.1, находим по таблице t-распределения при ν = 4 t_q = 2.13. Следовательно, коэффициент а₄ - незначимый. Поскольку матрица А^ТА диагональна, коэффициенты определены независимо и незначимый можно просто отбросить.

Итак, зависимость y от z₁, z₂ имеет вид:

y = 61.87 + 7.92z₁ + 3.42z₂.                        (4.11)

Поскольку z₁, z₂ принимают очень простые значения (-1, 0, 1), то адекватность модели проще проверить на этом этапе, до перехода к "старым" переменным х₁, х₂. Результаты сравнения эксперимента с расчетом по модели (4.11) - в таблице 4.2.

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

Номер опыта	1	2	3	4	5
yj эксп.	52	58.1	65.5	53.8	61.3
расч.	50.53	58.45	66.37	53.95	61.87
	2.152	0.123	0.752	0.022	0.321
Номер опыта	6	7	8	9
yj эксп.	70	56.2	65.9	74
расч.	69.78	57.37	65.28	73.20
	0.0467	1.362	0.379	0.638

В нашем случае число опытов - 9, число значимых коэффициентов модели - 3, число степеней свободы для вычисления дисперсии адекватности согласно (4.7) - 6; поделив на 6 сумму элементов последней стоки, получаем S_a² = 0.966.

В этом случае очевидно, что модель адекватна, т.к. S_a² < S_y², если бы S_a² оказалась больше S_y² следовало бы сравнивать их отношение с критическим значением F_q при числе степеней свободы к₁ = 6, к₂ = 4.

Теперь можно вернутся к модели (4.8). Подставляя (4.9) в (4.10.), имеем:

y = 61.87 + 7.92(x₁ - 60) / 20 + 3.42(x₂ - 4) / 2 = 31.3 +0.396x₁ + 1.71x₂ ,

т.е. b₁ = 31.28, b₂ = 0.396, b₃ = 1.71.

Если нас интересуют оценки дисперсий коэффициентов b, то их тоже легко найти, пользуясь свойством дисперсии и независимостью коэффициентов а.

,
,
.

4.2. Метод наименьших квадратов - общий случай

Кратко остановимся на случае, когда нет уверенности, что все погрешности независимы и имеют одинаковую дисперсию. В этом случае нужно проделать большую работу, выполнив многократные (р раз) измерения y_j во всех точках плана эксперимента и построить оценку ковариационной матрицы D погрешностей измерений по следующему правилу:

      (4.12)

Здесь i и j - номера экспериментальных точек, k - номера повторных изменений. - средние по к измерениям в i-ой ( j-ой) точке.

Видно, что D_ii - оценка дисперсии измерения в i-ой точке (р >> 1, поэтому р - 1» р). Обозначив W = D^-1 можно получить [4] следующие формулы вместо (4.4 и 4.5).

b = (A^TWA)^-1A^TW y,                                   (4.13)

S_bi² = [(A^TWA)^-1] _ii                             (4.14)

Вообще же матрица (A^TWA)^-1 является ковариационной матрицей оценок b, т.е. ее недиагональные элементы - оценки коэффициентов ковариации b.

Проверка значимости коэффициентов проводится так же, как в предыдущем случае (по 4.6).

Проверку адекватности модели в этом случае можно организовать как проверку гипотезы о том, что "истинное значение" в каждой экспериментальной точке равно рассчитанному по модели (см. приложение 1.3) по критерию Т².

Легко видеть, что если погрешности y независимы и одинаковы, то D = S_y²×I (I - единичная матрица) и (4.13, 4.14) переходят в (4.4, 4.5).

Подводя итог этому разделу, перечислим основные этапы решения задачи методом наименьших квадратов.

1. Формулируем модель, т.е. вид зависимости параметра от факторов.

2. Выбираем диапазон варьирования факторов (область плана эксперимента).

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

- обязательно делать опыты на границах варьирования факторов;

- обязательно делать опыты в центре плана (за исключением случая, когда мы уверены, что в модель входят только линейные по факторам члены, и это не нуждается в экспериментальной проверке);

- располагать экспериментальные точки симметрично относительно центра плана. В общем случае, имея в своем распоряжении ЭВМ, можно всегда подобрать набор экспериментальных точек, обеспечивающих максимальное значение определителя матрицы (A^ТA), что в свою очередь, минимизирует в совокупности [(A^ТA)^-1]_ii, при всех i.

4. Оцениваем дисперсию (ковариационную матрицу) результатов опытов путем повторных измерений при одинаковых условиях.

5. Проверяем гипотезу о ее диагональности и об однородности дисперсий y (см. приложение 1). Если сами условия эксперимента создают уверенность в независимости погрешностей у и одинаковости их дисперсий, этот пункт можно опустить.

6. Проводим эксперимент и оцениваем коэффициенты модели по формулам (4.4) или (4.13).

7. Оцениваем погрешности коэффициентов по (4.5) или (4.14).

Расчет по формулам (4.4, 4.13, 4.5, 4.14) даст правильные результаты и в том случае, если параметр не является нормально распределенной случайной величиной, т.к. при выводе этих формул такое предположение не иcпользовалось, а вот проверки значимости коэффициентов и адекватности модели при "не нормальном" распределении параметра выполнить нельзя, поскольку в этих проверках используются критерии, распределенные по законам,"производным от нормального".

8. Проверяем значимость коэффициентов по t-критерию (4.6).

9. Проверяем адекватность модели по F-критерию (4.7).

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

далее ...