CC BY
70
П. В. Шпилев
АЛГОРИТМ ПОСТРОЕНИЯ ОПТИМАЛЬНЫХ ПЛАНОВ ДЛЯ ТРИГОНОМЕТРИЧЕСКОЙ РЕГРЕССИИ
1. Постановка задачи
Тригонометрические регрессионные модели широко используются для описания периодических явлений в различных областях. Традиционными сферами применения являются медицина и биология (см. [1]).
Хорошо известно, что использование оптимального плана во многих случаях позволяет существенно повысить точность оценок параметров регрессионной модели. Для оценки параметров обычно используется метод наименьших квадратов. В большинстве случаев исследователи оценивают сразу все параметры модели, при этом оптимальными оказываются планы, сосредоточенные в равноотстоящих точках с равными весами. Однако во многих биологических задачах большое значение играют отдельные параметры, которые необходимо оценить с наибольшей точностью. Недавно (в работах [1], [2]) была решена задача нахождения оптимальных планов для оценки индивидуальных коэффициентов на полном интервале планирования [- п, п]. Однако на практике часто бывает невозможно проводить наблюдения на всем интервале [- п,п]. Данная статья посвящена численному методу нахождения планов, оптимальных для оценивания индивидуальных коэффициентов тригонометрической регрессии на интервале планирования [-а, а] при 0 < а < п. Метод основан на использовании функционального подхода (см. [3]).
Рассмотрим тригонометрическую регрессионную модель
У
= во + л2 @2j -1 sin(jt) + л2 @2j cos(jt) + £, t £ [ - а, а], а £ (0 п]. (1.1)
Определим в = (во, в1,..., в2т)Т как вектор параметров и
f (t) = (1, sint, cost,..., sin(mt), cos(mt))T = (fo(t),..., f2m(t))T
как вектор регрессионных функций. Следуя подходу, предложенному в [2], будем искать оптимальный план на множестве симметричных планов. Определим симметричный план £(а, п) как дискретную вероятностную меру
£( )= f-tn(0) -tn-Ла) ... - А(а) to^) ... tn-Ла) tn( а)\ л2)
, \ли(а) шп-1 (а) ... шо(а) шо(а) ... шп-1 (а) ип(а)] , .
где 0 < to(a) < ... < tn( а) < а, щ(а) > 0л^ п=о wi(a) = 1/2. Под информационной матрицей, как обычно, понимается матрица
M (€ ( а )) = [ f (t)fT (№((а, n)).
- a
Будем говорить, что параметр в к оцениваем для плана £(а, п), если существует линейная несмещенная оценка этого параметра, основанная на результатах наблюдений
П. В. Шпилев, 2007
©
в точках носителя этого плана. Пусть вк £ К2т+Х определяет единичный вектор (к = 0,..., 2 т) и А- есть обобщенная обратная матрица для матрицы А £ Е2т+1х2т+1. Тогда план £*(а, п) называется вк-оптимальным или оптимальным для оценивания коэффициента в к, если £*(a,n) минимизирует функцию Фк(£(а, п)) = вТМ-(£(а, Л)вк на множестве всех планов £ (а,п) вида (1.2), таких, что в к оцениваем для плана £(а,п).
2. Метод численного нахождения оптимальных планов на интервале [-а, а], а £ (0, п]
В основе метода лежит функциональный подход, описанный в монографии [3]. Его идея заключается в следующем.
1. Определяем число точек в оптимальном плане, которое не меняется в некоторой заданной области значений а £ (а*, Ъ*] С (0, п].
2. Выбираем вектор-функцию т (а,п), состоящую из нетривиальных точек и весов оптимального плана ( т(а, п) = ,'о(а),..., ЬпЛ(а), шо(а),..., Шг1 -1 (а))Т).
3. Составляем систему уравнений, которым удовлетворяет оптимальная вектор- функция т* (а,п). В качестве таких уравнений можно взять необходимые условия экстремума.
4. Проверяем невырожденность матрицы Якоби системы уравнений.
5. Выбираем точку ао, в которой матрица Якоби невырожденная, и находим оптимальный план в этой точке. В силу теоремы о неявном отображении существует аналитическая вектор-функция т* (а,п), являющаяся решением системы, которая может быть разложена в сходящийся ряд Тейлора в окрестности точки ао. Для нахождения такого выражения используется рекуррентная формула из работы [3].
Метод численного нахождения еА-оптимального плана
1. На первом шаге полагаем р : = 0, пр := т, Ьпр (а) : = а, ар : = 0. В окрестности нуля значения вектор-функции т* (а,пр) найдены [1]. 1.1. Находим точку а*+1:
Коли ГХ,-=п (Л) полагалм j О. пр | , пр и пороходим к путткту 3.
2.1. Если1 11 = 0 и*(а) = 0, то существует набор из 1р + 1 (0 < 1р < пр) различных величин 0 < вл < ... <
Бур < пр, для которых выполнены равенства А.$^) (ар+1) = 0, I = 0,..., 1р. Полагаем Пр+1 := пр, 7 : = 1р + 1.
2.2. Если 7 > 0, полагаем £{ (а,пр) : = £{ (а,пр) 1 {±Ь*(р) (а) ,ш* р (а)}, пр+1 := пр+1 - 1,
в] в]
7 :=7 — 1. Если7 = 0, повторяем процедуру, описанную в пункте 2.2. 134
3. ЕслиП "=1 \Ь*(а)-Ь*-1(а)1 = 0, то полагаем 7 := 7 + 1. Если 1Ь* (ар+1) —Ь*-1(ар+1)\ = 0, то и(а,пр) : = £*(а,пр) \ {±Ц-1{а), и *—1{а)}, и](а> := и,|(а) + А]-1(а>’ пр+1 := пр+1 — 1. Если 7 < пр, повторяем процедуру, описанную в пункте 3.
4. После перенумерации точек и весов плана £k(a,np) полагаем £1(а,пр+г) := £1 (а, пр), р := р +1. Если для функции т* (а, пр) матрица Якоби соответствующей системы не вырождена в точке ар, строим с помощью функционального подхода разложение функции т* (а, пр) в ряд Тейлора в окрестности точки ар . Переходим к пункту 1.1. Повторяя описанную процедуру, за конечное число шагов (р < 4т/3+1) численно находим оптимальный план.
Замечание. Применять вышеописанный метод для нахождения оптимального плана можно только в том случае, если якобианы соответствующих систем не равны нулю во всех точках а1,..., ар. Можно доказать, что это верно по крайне мере для е21 и б2г -1-оптимальных планов при
Были численно найдены вк-оптимальные планы для тригонометрических регрессионных моделей различных порядков. В связи с громоздкостью записи этих планов ограничимся двумя примерами.
Пример 1. Рассмотрим тригонометрическую модель третьего порядка (m = 3.)
1. а, G (0, А]. Оптимальный план £| (а, 3) имеет вид (1.2) и найден в [1]. B этом случае точки плана определяются по формуле t*(a) = arccos а-------ggAM cog (ш.А + i+cos(a)~j а а _
n 0 */ \ 1+2cos( a ) * / \ 5+7cos( a ) * / \ 7+5cos(a) */ \
U,...,3, а веса ujn(a) = igig—V4, wna = oa oa —V4, u>o(a) = oa oa —V4, ujn(a) =
' ' ' 0 V / 18+18cos(a)' 1V/ 36+36cos(a)' 2 1 j 36+36cos(a)' 3V j
2+cos(a) 18+18cos(a) •
2. a О [Ар,агееов(-|)]. Оптимальный план Q(a,2) имеет вид (1.2), где t2(a) = а. Веса плана в этом случае будут ш*(а) = 1/12, А*(а,) = 1/4, ш*(а) = 1/6. С помощью функционального подхода найдены разложения в ряд Тейлора функций 10 (а) и И (а) в точке а* = А-. Коэффициенты разложения этих функций в ряд указаны в таблице 1.
Таблица 1.
i 0 1 2 3 4 5 6 7
*o(e) 0.8957 -1.3868 -0.3695 -0.6236 -1.0123 1.8402 -3.4751 -6.8496
t*(a) 1.6961 0.2182 -0.5999 -0.3629 0.3038 0.1714 0.0005 0.00002
3. а О [агееов(-|), 7г]. Оптимальный план £|(а, 2) имеет вид (1.2), tg(a) = 0, t2(a) = а, ш1(а) = 1/4. С помощью функционального подхода найдены разложения в ряд Тейлора функций Шд(а) и t*(a) в точке а2 = агееов(-|). Коэффициенты разложения этих функций в ряд указаны в таблице 2. Кроме того, выполнено равенство ш* (а) = 1/4 — ш*(а).
Т аблица 2.
i 0 1 2 3 4 5 6 7
w О < 1/6 -0.1389 0.1389 -0.054 0.0334 -0.014 0.0071 -0.0031
t*(a) 1.7721 -0.6124 0.4465 0.008 0.0679 -0.0682 0.0349 -0.035
Пример 2. Проиллюстрируем на другом примере эффективность вк-оптимального плана. На практике эксперимент обычно проводится в равноотстоящих точках. Рас-
is- hx а2
\ о о о о
1G- О О *+-»■+-
14- О О е 5> О Ф О □ □ □ □
12- 0 о в
10- о О 0
8- о о с» 0
6 о
4j 2- 1 1 1 1 _ о ■ 1 1 ■ 1 1 ■ ■ 1 1 ■ 1 1 ■ 1 1 ■ 1 1 ■ > о 0 1 ■ 1 1 1 ■ ■ 1 1 '
и 0,5 1 1.5 2 2.5 3
Рис.1. Поведение функций эффективности плана £Q (а,п) по сравнению с планами £г(а), I = 0, 1,2,3 на интервале а О (0^, aQ = агееов(-3/5), ^ « 2.4712 .
смотрим тригонометрическую модель четвертого порядка (т = 4) и планы / _п па а п а па
е(а) = ( J+- _-г_ пГ 1
(2/1+1) (2п+1) (2п+1) (2п + 1) (2п+1) ''' (2п + 1) (2п + 1)
где п = 1 + 3, 1 = 0,1, 2, 3.
Функция эффективности численно найденного плана £6 (а,1 + 3) по сравнению с планом £‘ (а), I = 0,1, 2, 3 на интервале а £ (0, п] определяется отношением
Результаты, представленные на рис. 1, показывают, что построенный план более чем в 2 раза эффективнее планов, обычно используемых на практике. Это означает, что для достижения той же точности нам потребуется в два раза меньшее число экспериментов при использовании плана £6(а). Аналогичные результаты были получены и для других случаев. Необходимо отметить, что планы, оптимальные для оценки младших коэффициентов (т.е. е21 и еа-1 -оптимальные планы при I < т/3), менее эффективны. Это связано со спецификой этих планов (см. [2]).
Summary
P. V. Shpilev. The algorithm for constructing the optimal designs in the trigonometric regression model.
In the common trigonometric regression model the optimal design problem for the estimation of the individual coefficients on the partial symmetric intervals is investigated.
Литература
1. Dette H., Melas V. B. Optimal designs for estimating individual coefficients in Fourier regression models // Ann. Statist. 2003. Vol.31. N5 . P. 1669-1692.
2. Dette H., Melas V. B., Shpilev P. V. Optimal designs for estimating the coefficients of the lower frequencies in trigonometric regression models. Preprint Ruhr-Univ. Bochum, 2005. http:/www.ruch-uni-bochum.de/mathematik3/preprint.htm
3. Мелас В. Б. Общая теория функционального подхода к оптимальному планированию эксперимента. Изд-во СПбГУ, 1999.
Статья поступила в редакцию 17 мая 2007 г.
