Исследование максиминно эффективных планов для модели Михаэлиса—Ментена Текст научной статьи по специальности «Математика»

УДК 519.24 Вестник СПбГУ. Сер. 1, 2006, вып. 2

В. Б. Мелас, Ю. М. Старосельский

ИССЛЕДОВАНИЕ МАКСИМИННО ЭФФЕКТИВНЫХ ПЛАНОВ ДЛЯ МОДЕЛИ МИХАЭЛИСА—МЕНТЕНА*

1. Введение. Регрессионная модель Михаэлиса—Ментена

E[Y\t] = ^-b, i€[0,io], (1)

где a и b — оцениваемые параметры, широко используется для описания физических и биологических феноменов (см. например, [1] и [2]). Задача построения планов эксперимента, локально оптимальных в смысле различных критериев, изучалась в работах [3], [4] и многих других. В настоящей работе рассматриваются стандартизированные максминно эффективные E-оптимальные планы (далее для краткости такие планы называются СММЭ-планами). Такие планы были численно построены в работе [5] на основе аналитического решения промежуточной задачи — построения стандартизированных локально E-оптимальных планов. Преимущество СММЭ-планов по сравнению с локально оптимальными планами заключается в использовании не точечных, а интервальных оценок для нелинейно входящего параметра b.

В настоящей работе мы исследуем СММЭ-планы на основе функционального подхода, введенного в работах [6], [7] и развитого в монографии [8]. Этот подход позволяет аппроксимировать опорные точки и весовые коэффициенты СММЭ-планов отрезками рядов Тейлора.

В следующем разделе мы сформулируем задачу и представим предварительные результаты. В разделе 3 формулируются основные теоретические результаты, а в разделе 4 приводятся результаты численного исследования. Доказательство основных теоретических результатов дано в приложении.

2. Постановка задачи. Модель (1) представляет собой специальный случай общей нелинейной регресионной модели

E(Y\t)_ ri(t,e),

где в _ (в\,.. .,вт)Т —вектор неизвестных параметров. Примем стандартное предположение о том, что ошибки измерений — независимые одинаково распределенные случайные величины. Под планом эксперимента будем понимать дискретную вероятностную меру

£ _ I Х1 ... ХП

У W1 ... Wn

где ti _ tj (i _ j) —точки, в которых проводятся наблюдения,а Wi —весовые коэффици-

n

енты: Wi > 0, wi _ 1. Пусть N — общее количество экспериментов. Будем говорить,

i=i

что эксперимент проводится в соответствии с планом если в точках ti проведено r

* Работа выполнена при частичной финансовой поддержке РФФИ (грант №06-01-00284). © В. Б. Мелас, Ю. М. Старосельский, 2006

экспериментов, где т^ = ^.¿N1 или т^ = ^.¿N1 + 1 так, что ^ т^ = N. Пусть оценивание

i=1

параметров производится по методу наименьших квадратов,

N ri

1

в

ê(N ) = - V(ti,ê))2 ,

i=ij=i

где У^ —результат ] -го измерения в точке г = 1,...,п; ] = 1

Обозначим через в* истинное значение параметра в. Хорошо известно (см., например, [9]), что при некоторых условиях регулярности (эти условия нетрудно проверить для модели (1)) при N —> оо вектор уПУ(в^ — в*) имеет асимптотически нормальное распределение N(0, а2О), где а2 —дисперсия ошибки измерения, а О = М-1(£,в) \в=в*, где

— информационная матрица Фишера.

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

Рассмотрим регрессионную модель (1). Информационная матрица Фишера для этой модели имеет вид

M ({, a, b) = DaJ f (t, b)f (t, b)T {(dt) Da,

где

Da = diag{l,a}, /М) = (ïfb, - (^2) . Следуя работе [10], введем стандартизированный критерий ^-оптимальности. Через A-будем обозначать матрицу, обобщенно обратную для матрицы A. Положим dT A-d = ж, если вектор d не принадлежит подпространству, натянутому на столбцы матрицы A.

Обозначим eT = (1,0), eT = (0,1). Через (a,b), j = 1,2, обозначим локально ej-оптимальный план эксперимента, т.е. план, минимизирующий величину

eT M - (С, a, b)ej

при фиксированных значениях a и b, являющихся начальными приближениями для истинных значений a*, b*. Положим

Ка,ь = diag {(efM-(£î, a, 6)ei)-*, (е^М-(£2*, а, Ь)е2)-^} .

Заметим, что матрица

C (i,b) = (KlbM-({, a, b)Ka,b)-1

не зависит от a. План, максимизирующий величину Amin(Cb)), будем называть стандартизированным локально ^-оптимальным планом.

Для преодоления зависимости критерия от начального приближения для параметра b, следуя работе [5], введем следующее определение.

Определение 1. План, .максимизирующий величину

. Аmin(C(g, Ъ)) mm --—--——,

ье[ь1,Ъ2] max Amin(C(п, b)) п

где [bi,b2] — интервал возможных значений параметра b, будем называть стандартизированным максиминно эффективным E-оптимальным планом (или, для краткости, СММЭ-планом).

В работе [5] стандартизированные локально E-оптимальные планы для модели (1) найдены в явном виде. Кроме того, доказано, что при достаточно малой величине b2 -bi СММЭ-план имеет вид

Ч 1 1 - w) <2>

и максимизирует величину

min aAmin(C(£, bi)) + (1 - a)Amin(C(£, b2)).

ae[0,i]

В следующем разделе мы исследуем СММЭ-планы с помощью функционального подхода, описанного в [8].

3. Аналитические свойства СММЭ-планов. Фиксируем одну из границ bi и будем изучать зависимость величин t и w, входящих в формулу (2), от границы z = b2. Аналогичные результаты получаются, если зафиксировать границу b2 и z = bi. Пусть

Ф(а, z) = aAmin(C(£, bi)) + (1 - a)Amin(C(£, z)).

t to

Обозначим и = (м1,м2,иэ) = (¿,и),а), — ( ^ ^ ^

Ъ(и,г) = Ф(а,£и,%).

Рассмотрим систему уравнений

д

— Щщг)= 0, г = 1,2,3. (3)

дщ

Пусть и(г) = (1(г), й(г), а (г)) —некоторое решение этой системы. Через . (г) обозначим матрицу Якоби системы (3), оцененную для этого решения:

J{z)=fd2*b*r3

\ дщдип /•• л

\ 1 j У i,j = i

. (4)

u=U(z)

Справедлива следующая теорема

Теорема 1. Если 0 < г — &1 <6, где 6 —достаточно малая положительная величина, то система (3) имеет единственное решение такое, что 0 < и1 < ¿о, 0 < щ < 1,г = 1, 2, причем матрица Якоби (4) обратима, а функции %Ь, и а, сопоставляющие г единственное решение системы (3), являются вещественными аналитическими функциями. Кроме того, план

£ = £и = ( Ъ 1- Ъ

где и = и{г) = {Ь{г),Ъ{г),а{г)), есть СММЭ-план для .модели (1) на отрезке [61,62] = [6иг].

Доказательство этой теоремы будет дано в разделе 5.

Хотя теорема утверждает только факты для малых 6, численные расчеты показывают, что при 61 = 1 и 6 < 7 матрица J{г) обратима. Отсюда вытекает, что все утверждения теоремы при 61 = 1 справедливы для 6 < 7. Далее для упрощения обозначений знак «волна» у функции и{г) будем опускать.

Пусть и{г) = и{г), ио = и {го), и для любой вещественно аналитической функции ¥{г)

1

в

в > 1,

где ¿о > 61 некоторая фиксированная точка.

Разложение функции и{г) в ряд Тейлора в окрестности этой точки в наших обозначениях принимает вид

и{г) = П(г){г - го). i=о

Коэффициент и (о) может быть найден численно, путем решения системы (3) при 2 = го. Остальные коэффициенты могут быть найдены численно с помощью общих рекуррентных формул, введенных в работе [11] (см. также [8, гл. 2]). Для системы (3) эти формулы принимают вид

и(в) = Jo1{го){g{u<s-l>{г), г))(s), (5)

где

' \ ди1 ' ди2 ' диз

к-1

и<к-1> {г) = ^2 и(г){г - го)1. г=о

Формула (5) и теорема 1 позволяют предложить следующий алгоритм построения и исследования СММЭ-планов.

1о Численно решаем систему уравнений (3) для некоторого значения г = го >61. Полагаем ко =0.

2о С помощью рекуррентных формул (5) находим и(к) при к = ко + 1.

3о Проверяем для плана £ = £и, и = и<к> {г) выполнение с заданной точностью условий теоремы 3.3 из работы [5] при интересующих нас значениях г. Если условия не выполняются, то полагаем ко ■= ко + 1 и возвращаемся к шагу 2о.

В следующем разделе мы продемонстрируем работу этого алгоритма.

4. Численные результаты. Рассмотрим случай 61 = 1, г = 62,Ьо = 10,го = 5 (для других случаев мы получили аналогичные результаты, которые опускаются ради краткости). С помощью пакета Ма^аЬ находим

и(о) = {11757, 0.4550, 0.5575).

г=г0

После трех итераций алгоритма получаем коэффициенты и^), 3 = 1, 2, 3. Значения этих коэффициентов приведены в таблице 1. Планы, построенные с помощью этих коэффициентов для г =1, 2, 3, 4, 6, 7, 8, приведены в таблице 2. Все эти планы удовлетворяют необходимым и достаточным условиям, накладываемым на СММЭ-планы с точностью до 10~5. Отметим, что они также совпадают с планами, численно найденными в работе [5]. Заметим, что при г > 8 СММЭ-планы уже не имеют вид (2), а число их опорных точек равно 3. Однако планы, построенные с помощью описанного алгоритма, и при г > 8 имеют минимальную эффективность

„ • А тш(С(£,Ь)) ей: = тт

ье[ь1,ь2] тах Атт(С(п, &))' п

близкую к минимальной эффективности СММЭ-планов, найденных численно в работе [5], см. табл. 3.

Таблица 1. Коэффициены и(о), М(х), «(2)> и(3) разложения и* (г) = и(в)(г — го)в при го = 5

з=0

и(0) «(1) «(2) и(3)

1.1757 0.4550 0.5575 0.0816 -0.0064 0.0015 -0.0084 0.0008 -0.0005 0.0009 0.0001 0.0001

Таблица 2. Опорные точки и веса СММЭ-планов для различных г

г 1 2 3 4 6 7 8

хх 0.6045 0.8169 0.9693 1.0847 1.2497 1.3111 1.363

Х2 10 10 10 10 10 10 10

11)1 0.4844 0.5120 0.5274 0.5375 0.5506 0.552 0.5587

и>2 0.5156 0.4880 0.4726 0.4625 0.4494 0.448 0.4413

Заметим, что на практике обычно используются планы в равноотстоящих точках. Численное исследование показывает, что минимум эффективности таких планов достигается на левой границе при Ь = Эта минимальная эффективность растет с увеличением числа опорных точек, но довольно медленно. Для плана в 15-и равноотстоящих точках она равна 0.3685. Таким образом, применение СММЭ-планов позволяет уменьшить число экспериментов в 2 и более раза в зависимости от величины Ь2.

Таблица 3. Минимальные эффективности СММЭ-планов найденных численно (£*) и планов полученных с помощью рядов Тейлора (£2)

г 2 4 6 8 10 12 14 16

ей(Г) 0.9544 0.8451 0.7733 0.7270 0.7034 0.6906 0.6829 0.6779

0.9544 0.8451 0.7733 0.7253 0.6911 0.6657 0.6460 0.6302

Отметим, что используя коэффициенты, указанные в табл. 2, читатель может построить вручную план для интервала вида [Ь1,Ь2] = [1, г] при ¿о = 10. Первая опорная точка и ее вес вычисляются по формулам

х1 = £ « 1.1757 + 0.0816(2 - 5) - 0.0084(2 - 5)2 + 0.0009(г - 5)3, £1 = £ « 0.4550 - 0.0064(2 - 5)+0.0008(2 - 5)2 + 0.0001(г - 5)3.

Вторая опорная точка совпадает с ¿о = 10, а ее вес равен 1 - £1.

5. Приложение: доказательство теоремы 1. Решающий пункт доказательства заключается в проверке невырожденности матрицы Якоби.

Пусть сначала г — произвольное число, удовлетворяющее неравенству г > 61.

Обозначим т = {г,ъ), £Т = ( ^ ] ,

^ ' ^ ъ 1 — ъ)

Н{т, 6) = \тт{С{&,6)).

Непосредственное вычисление показывает, что матрица J = J{г) имеет следующую

структуру:

А й

J = V йТ 0

где

А= (дЧ{щг)

\ дщдиз л

\ 1 3 / ij=1

¿т = (¿№> ЬО - я(Т, г)], ¿[Я(г, 60 - Н(т, г)]) ,

и = и{г) = {и1, и2, из) —некоторое решение системы (3), т = {и1, и).

В силу известных формул Фробениуса (см., например, [9]), если det А = 0, й = {0, 0)Т, то

det J = — ^ А)йТА-1й = 0.

Покажем, что й = {0, 0)Т. Пусть, напротив, й = {0,0)Т. Тогда, как показывает прямое вычисление, производные по г и ъ функций

и3Н{т, 61) + {1 — и3)Н{т, г) и Н{т, 61) — Н{т, г) равны нулю при т = {и1,и2). Отсюда вытекает, что

Н {т,6)=0,Н'ш {т,6)=0 (6)

при 6 = 61 и 6 = г. В работе [5] доказано, что единственным решением уравнений (6) по т является вектор т{6) = {г{6),ъ{6)) такой, что план

Щ го

ъ{6) 1 — Ъ{6)

есть стандартизированный локально ^-оптимальный план. Из явной формулы [5, формула (3.4)] вытекает, что г{61) = г{г). Таким образом, выполнение равенства (6) для одного и того же вектора т невозможно. Полученное противоречие доказывает, что 6 = {0, 0)Т.

Пусть теперь 0 < г — 61 <6, где 6 достаточно мало. В силу аргумента непрерывности при 6 ^ 0 матрица А стремится к матрице Якоби для системы уравнений

дд -яМо = о,-яМ1) = о,

которые однозначно определяют единственный стандартизированный локально Е-оп-тимальный план. Обозначим эту матрицу ,/(Ь1):

J(b)

о2 \ 2

Я(т,Ь)

/ i,j=1

dridrj

t(b) to w(b) 1 - w(b) Обозначим

ir = (t(b),w(b))

— стандартизированный локально Я-оптимальный план.

Q(T,P2,b)

PT C(£r (b),b)P

PT p

где p = (1,p2)T, 0 < b < ж.

Пусть т* (b) = (t* (b),w*(b)) таково, что £r*(b) —стандартизированный локально ^-оптимальный план, (1,p2(b))T —собственный вектор матрицы C(£r*(b),b), соответствующий ее минимальному собственному числу. В силу результатов [5] т* (b) и p2(b) определены однозначно.

В силу свойств минимального собственного числа имеем

Q(t *(b),p*(b), b) = Amin (C (£r *(b), b)), д

— Q(T*(b),p2,b) = 0 при P2=P*2(b),

e

r

д

—Q(T,p*2(b),b)=0 при т = т*(Ь), г =1,2. (7)

Обозначим

д

4b)= ^(т*(Ь),р2,Ь) dp2

2

TQ(r*(b),p2,b)

p2=p*(b) p p

-jr- [e^C*e2 - A]

где

p = (1, p2(b))T, C = C(£r *(b),b), Amin = Amin(C), e2 = (0,1)T.

В силу единственности минимального собственного числа матрицы C (см. [5]), получаем, что d(b) > 0. Пусть

G(b)=(-s$F-Q(T,P2,bj) т = т*(Ъ), Р2=Р*2(Ь),

\ г j 2 ф)= (lSrMT>P^)i=1> Т = Р2=Р*2(Ъ),

92(b)= (q(b))2 .

Нам потребуется следующая Лемма. При любом b, 0 < b < ж

J(b) = G(b) - qib)qT{b)

d(b)

Доказательство этой леммы аналогично доказательству теоремы 2.4.2 в работе [8]. Непосредственное вычисление показывает, что

С{6) =

Оц {6) 0

det ,1{6) = —С11 {6)

Ыъ <1(Ъ)

где

РТ к-/{г*)

д2{6) = —2{К-1)

122-

РТ Р

г*

+

го

К** + 6)2 {го + 6)2\

г* = г*{6), р ={1,р22{6))Т.

Из (7) следует, что

г=г

[рТ к-1/{г*)]2

-'тр- = Ат;п(С(£т*(Ь), 6))

рТ р

и, значит, ц_2 {6) > 0.

Докажем, что Оц {6) = 0. Из соотношения (7) следует

9'{г* {6))=0, (8)

где д{г) = рТк-1/{г), р = {1,р2^{6))Т.

Пусть С1 = {К-1)11, С2 = {Кфцр*{6). Тогда имеем

С1г2 + С16г + С2 сф + С162 — 2с2

»(*) = (¿ + 6)2 ' 9 (*) = -^з-.

„( _ С1Ъ(г + Ъ)-Цсф + сф2 -2с2) 9 [ ) - (г + 6)4

Из (8) получаем соотношение

С1г* {6)6 + С162 — 2с2 =0

и, значит,

Отсюда вытекает, что О11{6) =0 и det Л{6) = 0.

Следовательно, матрица .1 {г) обратима, а система (3) имеет единственное решение для достаточно малых 6. Аналитичность функций г{г), Ъ{г), а{г) при достаточно малых 6 вытекает из теоремы о неявной функции (см., например, [12]).

Summary

V. B. Melas, Yu. M. Staroselsky. The maximin efficient designs for the Michaelis—Menten model.

The standardized maximin efficient B-optimal experimental designs are studied by the functional approach. This approach allows one to approximate the support points and weight factors of designs by the Taylor series.

Литература

1. Beverton R. J. H., Holt S. J. On the Dynamics of Exploited Fish Populations. London, 1957.

2. Johansen S. Functional relations, Random Coefficients and Nonlinear Regression, with application to Kinetic Data. New York: Springer, 1984.

3. Duggleby R. G. Experimental designs for estimating the kinetic parameters for enzymecatal-ysed reactions // J. Theoret. Biology, 1979. Vol.81. P. 671-684.

4. Dunn G. Optimal designs for drug, neurotransmitter and hormone receptor assays // Statist. Medicine, 1988. Vol. 7. P. 805-815.

5. Dette H., Melas V., Pepelyshev A. Standardized B-optimal designs for the Michaelis—Menten model // Statistica Sinica, 2003. Vol. 13. P. 1147-1163.

6. Melas V. B. Optimal designs for exponential regression // Math. Operationsforsh. Statist. 1978. Vol. 9. P. 45-59.

7. Мелас В. Б. Оптимальное планирование эксперимента для экспоненциальной регрессии // Математические методы планирования эксперимента / Под ред. В. В. Пененко. Новосибирск: Наука, 1981. C. 174-198.

8. Melas V. B. Functional Approach to Optimal Experimental Design. Heidelberg: Springer, 2006. 336 p.

9. Ермаков С. М., Жиглявский А. А. Математическая теория оптимального эксперимента. М.: Наука, 1987. 320 c.

10. Dette H. B-optimal designs for regression models with quantitive factors — a reasonable choice // Canad. J. Statist. 1997. Vol.25. P. 531-543.

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

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

Статья поступила в редакцию 12 февраля 2006 г.