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

УДК 519.24

В. Б. Мелас, Л. А. Крылова

Вестник СПбГУ. Сер. 1, 2003, вып. 4 (№25)

О ФУНКЦИОНАЛЬНОМ ПОДХОДЕ К ИССЛЕДОВАНИЮ Е-ОПТИМАЛЬНЫХ ПЛАНОВ ДЛЯ ПОЛИНОМИАЛЬНОЙ РЕГРЕССИИ*

1. Введение

Проблема построения Е-оптимальных планов для полиномиальной регрессионной модели впервые рассматривалась Ковригиным [1]. Он нашел решение этой задачи на отрезке [—1,1]: в случае степени полинома большей или равной двум существует единственный Е-оптимальный план, который сосредоточен в экстремальных точках многочлена Чебышева первого рода. Пукельсгейм и Стадден [2] переоткрыли это решение в несколько более общем виде( рассматривались также усеченные Е-оптимальные планы). Хейлигерс [3] нашел, что это решение может быть распространено на симметричные отрезки, длина которых меньше некоторого критического значения, и на знакоопределенные отрезки путем естественного линейного преобразования точек плана. Случай симметричных отрезков, длина которых больше критического значения, был изучен в работах Меласа [4, 5] и Меласа и Крыловой [6, 7] на основе функционального подхода. Идея этого подхода была впервые сформулирована в работе [8] и заключается в исследовании точек и весов оптимальных планов как функций некоторых вспомогательных величин, неявно заданных условиями оптимальности. В работе [9] были введены общие рекуррентные формулы для разложения этих функций в ряд Тейлора. Эти формулы оказались более эффективными, чем формулы из работы [5].

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

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

2. Задача нахождения Е-оптимального плана и двойственная задача

Пусть результаты измерений описываются соотношением

где , ■■■, вт — оцениваемые параметры, е — случайная величина (ошибка эксперимента) такая, что Ее = 0, Ее2 = а2 > 0, причем ошибки измерений в различных экспериментах некоррелированы. Модель (1) называется моделью полиномиальной регресии на симметричном отрезке.

* Работа выполнена при финансовой поддержке Министерства образования РФ (грант РД02-

т

(1)

1.1-38)

© В. Б. Мелас, Л. А. Крылова, 2003

Под планом эксперимента будем, как обычно, (см., например, [10]) понимать дискретную вероятностную меру, задаваемую таблицей

. = ҐХ! ... Хп

І ші ... шп

где Хі = Хз (І = Д Хі Є [-г, г], Ші > 0, £ Ші = 1.

План £ называется ^-оптимальным планом (см. [4]), если он максимизирует величину

Ашіп(М (£)),

где Атіп(М) — минимальное собственное число матрицы М, М(£) — информационная матрица плана £,

/ п \т

М (£) = (53 Їі(хі )3 (хі )шА .

\1 = 1 ' і,3 = 1

Здесь и далее будем рассматривать случай т = 2к + 1, к > 1 (случай т = 2к формулируется аналогично, см. [4, 5]).

Как показано в [4] для модели (1) всегда существует единственный ^-оптимальный план и этот план (при т = 2к +1) имеет вид:

¡Е = £(г) = I-г - Х2 (г) ... -Хк(г) 0 Хк (г) ... Х2(г) г

\ Ші Ш2 ... Шк Шк+1 Шк ... Ш2 Ші

Пусть Н(у) — произвольный многочлен степени 2к с единичным старшим коэффициентом неотрицательный на интервале [0,г2]. Такой многочлен всегда может быть представлен в виде суммы квадратов двух многочленов (см. [11]):

Ну) = <А(У)+ УÍA(У),

где ¥1(У) = РТ(1,У,...,Ук), Р2(у) = ЯТ(1,y,...,yk-l), РТ = (Р0,...,Рк), ЯТ =

(qo,...,qk-l), Рк = 1

Рассмотрим задачу: найти векторы р и q, доставляющие минимум величины

Н(о а] _ тах У1Ы+УУ2Ы (2]

1И2 + 1Ы12 ’ (2)

среди всех векторов p Є Rk+1, q Є Rk, pk = 1.

Пусть p = p(r) = (p0,pi/r2,... ,pk/r2k)T, p* = (p*, ... ,p*k)T — вектор нормированных

ненулевых коэффициентов многочлена Чебышева первого рода степени 2k:

p*T(1, х1, ..., x2k)T = cos(2k arccosx),

0 < xk < x*k_ 1 < ■■■ < x*l < 1 — неотрицательные экстремальные точки этого многочлена.

Обозначим через £ch = £ch(r) следующий план эксперимента:

, , (—r -rx* ... -rx*k 0 rx*k ... rxl r

clh \ W1 W2 ... Wk Wk + 1 Wk ... Ul W1

где вектор w = w(r) = (wk+1(r), Wk(r),..., W1(r))T равен

w(r) = NFr 1p(r)/\\p(r)\\z,

N = diag{l, 1/2,. - -, 1/2}, Fr = ((rx*fc+2^2i~2(-l)J'+1),,ti,

V

k

x* = i, xk+i = °-Далее мы используем результаты работы [5].

При r < 1 минимальное собственное число матрицы M (ich(r)) является простым. Пусть г* — минимальное положительное г, для которого кратность этого числа больше или равна 2 (фактически она всегда равна 1 или 2). Тогда при r < r*, m > 2 существует единственный ^-оптимальный план

i(r) = £ch (r),

а задача минимизации величины H(p,q) имеет единственное решение p = r2kp(r)/pk, q = 0, причем при этих значениях p и q максимум в (2) достигается тогда и только тогда, когда y = r2(x*)2, i = i,...,k +1 и равен

Amin(M (£ch (r))) =i/\\p(r)\\2.

Можно проверить [4,5], что r* < ж при любом m.

Пусть теперь r > r*. Через Pr обозначим множество собственных векторов матрицы M(i(r)), соответствующих ее минимальному собственному числу.

При r > r* задача нахождения ^-оптимального плана и задача минимизации величины H(p, q) имеют единственное решение и также связаны соотношением двойственности [4]:

max Xmin(M(i)) = min H(p, q).

5 pq

Причем, если минимум в правой части достигается для (p,q), то p и q G Pr. Кроме того, максимум в (2) достигается для y = x2 (r), i = i,...,k +1, Xk+i =0, Xi = r. Для нахождения решения обеих задач при r > r* может быть использован функциональный подход, который мы опишем в следующем разделе.

3. О функциональном подходе к решению экстремальных задач

Пусть функция ф(в,г), где в G Kd, z G R задана, непрерывна по z и непрерывно дифференцируема по в при в G О, z G I, О — замкнутое множество в Rd, I — конечный интервал. Пусть эта функция является вещественно аналитической на множестве 1пШ х IntI и при z = zo G I, в = в(о) G IntQ удовлетворяет системе уравнений

д

—ф(в,г) = 0, i = l, (3)

При некотором дополнительном условии мы можем построить вектор-функцию e(z), являющуюся решением системы (3) в некоторой окрестности точки zo.

Для любой (скалярной, векторной или матричной) функции ^(z), вещественно аналитической в некоторой окрестности точки zo, обозначим

<Р(о) = ¥{s) = ^(s)(^o), s = 1,2,...,

s

V(s) (z) = ^Z ¥(t)(z - z0)t.

t=0

Пусть 1 (в, г) — матрица Якоби для системы (3),

(І

1(в, г)

d

1(и) = 1 (во,и), д(в, г)

Теорема 3.1. Предположим, что функция ф(в,х) обладает перечисленными свойствами и =0, г = 0, 1,...,в — 1, det ,1(3) = 0. Тогда справедливы следующее утвер-

ждения:

(I) В некоторой окрестности (V) точки хо существует вещественно аналитическая функция в(х) такая, что в(хо) = в(о),

(III) Если при х = хо решение в = в(о) системы уравнений (3) единственно, то при любом х € и решение этой системы единственно.

Утверждения (I) и (III) теоремы вытекают из теоремы о неявном отображении [12], более подробное доказательство можно найти в работе [9]. В этой же работе дано доказательство утверждения (II).

Для применения теоремы 3.1 к исследованию оптимальных планов нужно записать необходимые условия оптимальности в виде уравнения (3), решить (аналитически или численно) систему уравнений (3) в некоторой точке хо и проверить условие невырожденности матрицы Якоби.

Для случая ^-оптимальных планов полиномиальной регрессии на симметричном отрезке эти шаги были осуществлены в работах [4,5] и будут описаны в следующем разделе. Основная идея состоит в том, чтобы включить в вектор в параметры решения обеих задач — прямой и двойственной.

при любом фиксированном г Є и и в = в(г) и ряд

Ж

в(4)(г — гоУ

4=0

абсолютно сходится при г Є и.

(II) Коэффициенты в(4) могут быть вычислены по формуле

где

в<г-1>(г) = ^2 в{л(г - го).

з=о

4. Основное уравнение для В-оптимальных планов

Пусть х = 1/г2. Определим вектор в следующим образом:

вт = (гкр0, гк~1р1, грк-1, ^гк~1д0, ф,гк~2дъ ф,дк-\, 2шк/х,...,2ш\/х, ххХк, ...,хх2).

Пусть в(х) = в = в, где рг = рг, С = С, г = 0,...,к — 1, Хг = хг, г = 2,...,к, Шг = Шг г = 1,...,к, причем (р , д) — решение двойственной задачи.

Обозначим

-г -Х2 . . . -Хк 0 Хк ... Х2 г

Ші Ш2 ... Шк Шк + 1 Шк ... Ш2 Ші

Шк + 1 = 1 — 2 Ші

і=1

Л,„ , ртмт)р+ятм(т)я ,

'Ч#, г) = -----—■—-------, рк = 1,

Р1 р + я1 я

Р = (ро, 0,Р1,.. .,0,рк)Т, Я = (0,до, 0,Я1,..., дк-1, 0)т.

Нетрудно проверить, что

дв.,

■Х(в,г)=0, і = 1,...,4к — 1,

при в = в(х), если х достаточно мало.

Как показано в работах [4, 5] при

ф(в, х) = Х(в, х), хо = 0

справедливы следующие условия:

'(о) = 0, det '^(1) = 0

причем решение системы уравнений

единственно и найдено в явном виде.

Таким образом, к решению задачи нахождения ^-оптимального плана и двойственной к ней задачи можно применить теорему 3.1.

5. Построение и исследование рядов Тейлора

Для т < 4 задача нахождения ^-оптимального плана полностью решена в работе [5]. Однако, уже при т = 5 ряды, построенные в этой работе, не позволяют находить решение задачи для некоторых промежутков. В связи с этим мы построим некоторое переразложение этих рядов.

д

Пусть т = 5. Величина г* = 1/г*“2 определяется из соотношения

г;* = а^М |,г = —; Лтт{М(^ск(г))) имеет кратность 1| .

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

64г5 - 32г4 + 96г3 - 44г2 + 33г - 12 = 0,

« 0.396787. При г < г* = 1/л/г* имеем

((г) = Ыг) = {~Г -Г/^ 0 Г

I ^1 1^2 1^3 1^2 ^1

_ 4г2(1 + 2г2)

М1 " 64^4 + 64г2 + 1’ М2

М^(г))

16г2(1+г2)

1

64г4 + 64г2 + 1'

Исследуем поведение _Е-оптимальных планов при г < г* (г > 1 /%/**)• Точное значение вектора в (о), полученное с помощью теоремы 6.8.1 из [5] равно

в(0) = (1/4, -5/4, -3/4, 1, 8,1,1/4)Т,

Хо = Иш Атш(М(£(г))) = 1.

Т^Ж \ /

Матрица ,1(о) есть нулевая матрица, а матрица .1(1) равна

( 9 3 0 0 - 1/4 -1/4 -6

3 3/2 0 0 0 0 - 3/2

0 0 2 3/2 1/8 1/4 -2

0 0 3/2 9/8 - 1/32 1/4 - 3/2

-1/4 0 -1/8 -1/32 0 0 0

-1/4 0 1/4 1/4 0 0 0

-6 -3/2 -2 -3/2 0 0 - 3/2

к — ~ гкро, Р1 = гк- 1р1 1§ !> 7 к г >1 = гк 2^[хд\1 т =

VI = 2&1/г.

Первые коэффициенты разложения вектор-функции в (г) в ряд Тейлора в окрестности точки г = 0 представлены в таблице 1. Там же даны коэффициенты разложения функции А(в(г), г). Уже из этой таблицы видно, что коэффициенты рядов для весов (г>1 и ¿>2) и для минимального собственного числа (А) растут значительно быстрее остальных. Численное исследование показывает, что радиус сходимости для рядов V1 и V2 приблизительно равен 0.16, а радиус сходимости для остальных элементов вектора в приблизительно равен 0.35.

А

2

Коэффициенты разложений в точке zo = 0

0 1 2 3 4

Ро 1/4 -5/12 5/18 7/27 -146/81

Р1 -5/4 5/6 -5/9 -26/27 532/81

до -3/4 5/4 5/24 -31/72 4717/864

<71 1 -5/3 -5/18 31/54 -151/24

V2 8 -688/9 45632/81 -2691136/729 149801216/6561

VI 1 -41/9 1888/81 -96560/729 4945072/6561

У2 1/4 5/6 -5/9 -2/27 52/81

А 1 -9 56 -320 1808

Заметим, что величины Л, ) и щ можно выразить через остальные коэффициенты вектора д. Непосредственные вычисления показывают, что

Л =

к

ilPi

i=0

/k—i

+ I 2 qi

\i=0

^ pizii + — q2z2i i=0 i=0

Pk =

(4)

а при k = 2 (для нашего случая)

Vi =

Л(до - qiz2)

Vi(i - yi)(Vo + qiyi):

\(qiz2 - <703/2)

(1 - S/2)(<7o + <71 ) '

(5)

Таким образом, можно находить величины ро, р1, )о, ц\, У2 с помощью разложения функции д(г) в ряды Тейлора, а величины Л, и 1)2 по формулам (4)-(5).

Такой способ позволяет с высокой точностью находить точки и веса _Е-оптимального плана при г < 0.3 (г > 1/а/0.3).

Чтобы решить задачу в промежуточном случае 0.3 < г < 0.396, построим пере-разложение функции д(г) в окрестности точки г\ = 0.2. Описанным выше способом, используя 15 членов разложения функции в нуле, мы получим значение вектора д(г\), приведенное в первом столбце табл. 2. Остальные столбцы табл. 2 построены с помощью рекуррентных формул из теоремы 3.1.

Таблица 2

Коэффициенты разложений в точке zi = 0.2

i

i

0 1 2 3 4 5

Ро 0.17804 -0.30698 0.22962 -0.20355 -0.02974 0.36925

Pi -1.10688 0.60901 -0.43501 0.55512 -0.10179 -1.11856

до -0.48819 1.41862 1.00063 3.72679 11.81191 40.91730

<71 0.65139 -1.88521 -1.31832 -5.06216 -16.14959 -54.73347

V2 2.21396 -9.66086 36.65304 -121.46442 374.99630 -1108.795

VI 0.52635 -1.33370 2.49681 -6.18865 16.62481 -45.22451

У2 0.39473 0.61899 -0.48325 0.25700 0.20932 -0.34837

А 0.23768 -1.44590 5.21544 -15.39351 41.72046 -109.9943

Коэффициенты, приведенные в табл. 2, позволяют предположить, что ряды для P0, pii и yi имеют быструю сходимость. Численная проверка показывает, что при z = 0 и z = z* использование первых четырех коэффициентов позволяет находить векторы p(0) =

(1/4, —3/4,1), p(z*) = (g, —1,0) и точки уг(0) = 0.25, уг(-г*) = 0.50 с погрешностью, не превышающей 5 • 10~6. В частности, для y = y2 получаем

y{0)(z*) = 0.39473, y{i)(z * ) = 0.51654,

V(i)(z*) = 0.49782, y{3)(z * ) = 0.49978,

q(4)(z*) = 0.50009, y{5)(z * ) = 0.49999,

V(6) (z* ) = 0.50000, у(7) (z* ) = 0.50000,

тогда как истинное значение yi(z*) = 0.5. Поведение величин y<j>(0), j =0,1,...,7, а также величин pi<j> (z), z = 0,z*, i = 0,1, j = 0,1,...,7, вполне аналогично, что читатель может легко проверить.

Из соотношения МР = ЛР, где М = M(£(r)),

Р = (po, 0,pi, 0,1), Л = Лт1п(М(Ç(r)),

получаем, что веса оптимального плана и величина Л могут быть вычислены через po, pii и yi по формулам

2u)i = щг = A Ai, i = 1,2, uj3 = ХА3, А = 1

з :

ZA

j=i

А =(Л1,Л2,Лз) = Г diag{ 1, г ,г }р,

где Г = (рзУ1-1)3з,1=1, & = (1, У'2)р5, р = 1)Т, у1 = 1, уз = °.

Погрешность вычисления этих величин при использовании четырех первых коэффициентов из табл. 2 при г = 0, г = г* также не превышает 5 • 10-6. Хотя это не гарантирует такой же точности для всех остальных значений 0 < г < г *, выборочная проверка подтверждает высокую точность вычислений.

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

Таблица 3

Коэффициенты разложения величин Р0, р1 и У2 в ряд Тейлора в окрестности точки ^ = 0.2

0 1 2 3

ро 0.17804 -0.307 0.230 -0.204

Р1 -1.10688 0.609 -0.435 0.555

У2 0.39473 0.619 -0.483 0.257

Аналогичные таблицы можно построить и для случая т > 5.

6. Заключительные замечания

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

Summary

Melas V. B., Krylova L. A. On a functional approach to Я-optimal design for polinomial regression

The paper is devoted to a brief description of a functional approach to optimal experimental designs and its application to Я-optimal designs for polynomial regression models. The model of four degree is thoroughly studied.

Литература

1. Ковригин А. Б. Построение оптимальных Я-планов. Деп. в ВИНИТИ, N 3544-79.

2. Pukelsheim, F., Studden, W. Я-optimal designs for polynomial regression // Ann. Statist. 1993. Vol. 21. N 1. P. 402-415.

3. Heiligers B. Я-optimal designs in weighted polynomial regression // Ann. Statist. 1994. Vol. 22. P. 917-929.

4. Мелас В. Б. Я-оптимальные планы эксперимента. СПб., 1997.

5. Melas V.B. Analytical theory of Я-optimal designs for polynomial regression // Advances

in Stochastic Simulation Methods, 2000. P. 85-115.

6. Мелас В. Б., Крылова Л. А. Я-оптимальные планы для кубической регрессии на симметричном отрезке // Вестн. С.-Петерб. ун-та. Сер. 1. 1996. Вып. 3 (№15). С. 26-30.

7. Melas V. B., Krylova L. A. On Taylor series for points and weights of Я-optimal design for polynomial regression // Proceedings of the 3rd St. Petersburg workshop on simulation. Ed. by S.M. Ermakov, Y. N. Kashtanov and V.B. Melas. St.Petersburg University Press. 1998. P. 191-196.

8. Melas V. B. Optimal designs for exponential regression // Math. Oper. Forsch. Stat., Ser. Statistics, 9. 1978. P. 45-59.

9. Мелас В.Б., Пепелышев А. Н. Степенные разложения неявных функций и локальнооптимальные планы эксперимента // Статистические модели с приложениями в эконометрике и смежных областях. СПб., 1999. C. 108-117.

10. Федоров В. В. Теория оптимального эксперимента. М., 1971.

11. Карлин С., Стадден В. Чебышевские системы и их применение в анализе и статистике.

М., 1976.

12. Ганнинг Р., Росси Х. Аналитические функции многих комплексных переменных. М., 1969.

Статья поступила в редакцию 16 января 2003 г.