Построение оптимальных и локально оптимальных планов для регрессионных моделей

Корчажников Федор Васильевич

Построение оптимальных и локально оптимальных планов для регрессионных моделей

Корчажников Федор Васильевич. Построение оптимальных и локально оптимальных планов для регрессионных моделей. Научный руководитель: д. ф.-м. н., профессор В. Б. Мелас. Направление математика, механика, кафедра статистического моделирования. В работе рассматривается задача поиска L-оптимальных и локально D-оптимальных планов для регрессионных моделей. Найдены L-оптимальные планы для квадратичной и кубической модели, локально D-оптимальный план для обобщенной модели Михаэлиса-Ментен. Рассмотрен алгебраический подход в задаче поиска локально D-оптимального плана для дробно рациональной модели с четыремя параметрами и достаточно больших промежутков планирования, который был продолжен на случай малых промежутков. Предложен численный метод построения плана для этого случая, реализован алгоритм на языке R. Количество использованных источников: 14. Корчажников, Ф.В. Построение оптимальных и локально оптимальных планов для регрессионных моделей: бакалаврская работа / Корчажников Федор Васильевич. - СПб., 2016. - 35 с. - Библиогр.: с. 30-31.

Математика

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36425f1be77c40d58acc

UUID: b60931bb-e115-4e61-a6c2-52d36c09f8e9

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 245

Корчажников Федор Васильевич

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 553280 bytes

Поделиться работой

Санкт-Петербургский государственный университет Прикладная математика и информатика Вычислительная стохастика и статистические модели Корчажников Федор Васильевич Построение оптимальных и локально оптимальных планов для регрессионных моделей Бакалаврская работа Научный руководитель: д. ф.-м. н., профессор В. Б. Мелас Рецензент: к. ф.-м. н., доцент П. В. Шпилев Санкт-Петербург 2016

Saint Petersburg State University Applied Mathematics and Computer Science Computational Stochastics and Statistical Models Korchazhnikov Fedor Vasilyevich Constructing optimal and locally optimal designs for regression models Bachelor’s Thesis Scientific Supervisor: Professor V. B. Melas Reviewer: Associate Professor P. V. Shpilev Saint Petersburg 2016

Содержание Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. 4 Основные понятия . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1. План эксперимента . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2. Критерии оптимальности . . . . . . . . . . . . . . . . . . . . . . . 6 Глава 1. Нахождение L-оптимальных планов для полиномиальных мо делей . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.1. Случай квадратичной модели . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2. Случай кубической модели . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Глава 2. Нахождение локально D-оптимальных планов для нелинейных моделей . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1. Обобщенная модель Михаэлиса-Ментен . . . . . . . . . . . . . . . . . . . 11 2.1.1. Число точек плана . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.2. Весовые коэффициенты плана . . . . . . . . . . . . . . . . . . . . 12 2.1.3. Опорные точки плана . . . . . . . . . . . . . . . . . . . . . . . . . 13 Дробно-рациональная модель с четырьмя параметрами . . . . . . . . . . 15 2.2.1. Дифференциальное уравнение . . . . . . . . . . . . . . . . . . . . 15 2.2.2. Алгебраическое уравнение . . . . . . . . . . . . . . . . . . . . . . 17 2.2.3. Решение уравнения для достаточно большого промежутка . . . . 17 2.2.4. Решение уравнения для малого промежутка . . . . . . . . . . . . 21 Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Литература . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Приложение А. 32 2.2. Реализация алгоритма . . . . . . . . . . . . . . . . . . . . . 3

Введение Теория планирования — раздел математической статистики, занимающийся опти мальным планированием условий эксперимента. Основание этой теории было заложено Р. Фишером в 1935 году в его работе [1]. Результат эксперимента зависит от некоторых параметров и часто стоит задача нахождения оценки этих параметров или проверки некоторой гипотезы относительно них. Когда проведение опыта связано со значитель ными временными или материальными затратами, требуется осуществлять рациональ ный выбор плана эксперимента, чтобы за наименьшее число опытов получить наиболее точную оценку. Существуют множество критериев оптимальности, наиболее исследованным из ко торых является D-критерий, минимизирующий объем доверительного эллипсоида. Так же популярным является критерий L-оптимальности, который позволяет минимизиро вать среднюю дисперсию оценок параметров (определения соответствующих критериев будут даны в пункте 1.2). Целью данной работы является построение L-оптимального плана для квадрати ческой и кубической модели, а также D-оптимального плана для обощенной модели Михаэлиса-Ментен и дробно-рациональной модели с четырьмя параметрами (в слу чае большого и малого промежутков). Работа осуществлена в рамках темы СПбГУ 6.38.435.2015. Сначала будут рассмотрены L-оптимальные планы на отрезке [−1, 1] для некото рых полиномиальных моделей (D-оптимальные планы для них уже были найдены, см. [2, с. 64]). Далее в работе будут рассмотрены локально D-оптимальные планы для некоторых нелинейных по параметрам моделей. Например, для обобщенной модели Михаэлиса Ментен, которая находит широкое применение во многих областях науки: например, сельское хозяйство [3], биология сохранения живой природы [4], безопасность окружа ющей среды [5] — лишь некоторые из них. В биохимии с помощью этой модели опи сываются уравнения ферментативной кинетики — зависимости скорости протекания химической реакции от концентрации субстрата. Наконец, мы рассмотрим дробно-рациональную модель с четырьмя параметрами, для которой будет более подробно описан метод нахождения локально D-оптимального 4

плана, представленный в работе [6], а также найдены некоторые новые результаты, позволяющие численно находить локально D-оптимальный план в том случае, когда выразить его явно не представляется возможным. Данная модель также представляет практический интерес, поскольку применяется в химической кинетике [7], [8], а также в биологии и сельском хозяйстве [9], [10]. 1. Основные понятия В данном разделе мы ознакомимся с основными элементами теории, которые будут использованы далее при получении результатов. Они подробно рассмотрены в работах [2], [11]. 1.1. План эксперимента Результатом эксперимента будем считать набор вещественных значений 𝑦1 , . . . , 𝑦𝑁 . На практике этот набор часто может быть представлен в виде следующей модели: 𝑦𝑖 = 𝜂(𝑥𝑖 , 𝜃) + 𝜖𝑖 , 𝑖 = 1, . . . , 𝑁, (1) где 𝜂(𝑥𝑖 , 𝜃) — вещественная функция, зависящая от параметров 𝜃 = (𝜃1 , . . . , 𝜃𝑚 )T , 𝜖1 , . . . , 𝜖𝑁 — случайные величины, соответствующие ошибкам наблюдений, 𝑥1 , . . . , 𝑥𝑁 — условия проведения эксперимента из множества планирования X. Непрерывным (приближенным) планом эксперимента называют вероятностную меру ⎛ 𝜉=⎝ 𝑥1 . . . 𝑥 𝑛 ⎞ ⎠ , 𝑥𝑖 ∈ X, 𝑖 = 1, 2, . . . , 𝑛, (2) 𝜔1 . . . 𝜔𝑛 ∑︀ причем 𝑥𝑖 ̸= 𝑥𝑗 при 𝑖 ̸= 𝑗, 𝜔𝑖 ≥ 0, 𝑛𝑖=1 𝜔𝑖 = 1, 𝜔𝑖 — весовые коэффициенты, а 𝑛 — число попарно различных точек в плане. Под информационной матрицей приближенного плана будем понимать матрицу Z 𝑀 (𝜉) = 𝑓 (𝑥)𝑓 T (𝑥)𝜉(𝑑𝑥), (3) X где 𝑓 (𝑥) = (𝑓0 (𝑥), 𝑓1 (𝑥), . . . , 𝑓𝑚 (𝑥))T — вектор частных производных по параметру, 𝑓𝑖 (𝑥) = 𝜕𝜂(𝑥, 𝜃) . 𝜕𝜃𝑖 5

1.2. Критерии оптимальности План назовем невырожденным, если невырождена его информационная матрица. В таком случае для нее будет существовать обратная матрица, которую мы назовем дисперсионной 𝐷(𝜉) = 𝑀 (𝜉)−1 . Критерий D-оптимальности имеет вид log det 𝑀 (𝜉) → sup , (4) 𝜉∈Ξ𝐻 или log det 𝐷(𝜉) → inf , 𝜉∈Ξ𝐻 (5) где Ξ𝐻 - множество невырожденных непрерывных планов. Данный критерий соответ ствует минимизации объема доверительного эллипсоида: ̂︀ T 𝐷(𝜃̃︀ − 𝜃) ̂︀ ≤ 𝛼}. {𝜃̃︀ : (𝜃̃︀ − 𝜃) Критерий L-оптимальности имеет вид 𝑡𝑟𝐿𝐷(𝜉) → inf , 𝜉∈Ξ𝐻 (6) где 𝐿 - фиксированная неотрицательно определенная матрица. В задаче поиска L-оптимального плана для проверки найденного решения мы будем пользоваться теоремой эквивалентности [11, с. 112]: Теорема 1. Пусть множество информационных матриц M компактно и существу ет невырожденный L-оптимальный план 𝜉 * . Тогда следующие утверждения эквива лентны: а) 𝜉 * = 𝑎𝑟𝑔 𝑚𝑖𝑛 𝑡𝑟 𝐿𝐷(𝜉); 𝜉∈Ξ б) 𝑚𝑎𝑥 𝜙(𝑥, 𝜉 * ) = 𝑡𝑟 𝐿𝐷(𝜉 * ), где 𝜙(𝑥, 𝜉 * ) = 𝑓 𝑇 (𝑥)𝐷(𝜉)𝐿𝐷(𝜉)𝑓 (𝑥) 𝑥∈X При этом в точках 𝑥𝑖 L-оптимального плана 𝜉 * имеет место равенство 𝜙(𝑥𝑖 , 𝜉 * ) = 𝑡𝑟 𝐿𝐷(𝜉 * ). 6 (7)

Глава 1 Нахождение L-оптимальных планов для полиномиальных моделей 1.1. Случай квадратичной модели Начнем с поиска L-оптимального плана на отрезке. Рассмотрим один из простей ших случаев: полиномиальную модель с количеством параметров 𝑚 = 3, т.е. 𝑓 (𝑥) = (1, 𝑥, 𝑥2 )𝑇 . Матрица 𝐿 = 𝐼 — единичной матрице. Если в качестве интервала планирова ния взять X = [−1, 1], тогда из соображений симметрии ясно, что оптимальным будет план вида ⎛ 𝜉2 = ⎝ −1 1−𝜔 2 ⎞ 0 1 𝜔 1−𝜔 2 ⎠. Для того, чтобы найти 𝜔, выпишем информационную матрицу: ⎛ ⎞ 1 0 1−𝜔 ⎜ ⎟ ⎜ ⎟ 𝑀 (𝜉2 ) = ⎜ 0 1−𝜔 0 ⎟, ⎝ ⎠ 1−𝜔 0 1−𝜔 и обратную к ней 𝐷(𝜉2 ) = 𝑀 −1 (𝜉2 ), дисперсионную матрицу: ⎛ ⎞ 1 1 0 −𝜔 ⎜ 𝜔 ⎟ ⎜ ⎟ 1 𝐷(𝜉2 ) = ⎜ 0 1−𝜔 0 ⎟ ⎝ ⎠ 1 1 − 𝜔 0 𝜔(1−𝜔) След дисперсионной матрицы равен 𝑡𝑟 𝐿𝐷(𝜉2 ) = 𝑡𝑟 𝐷(𝜉2 ) = 2 . 𝜔(1 − 𝜔) Согласно критерию оптимальности, находим inf этой функции — он достигается при 𝜔 = 21 . Следовательно, L-оптимальный план для данной задачи выглядит так: ⎛ ⎞ −1 0 1 ⎠. 𝜉2* = ⎝ 1 4 1 2 1 4 План 𝜉2* был найден, исходя из некоторого предположения о симметрии располо жения точек 𝑡𝑖 ∈ X. Поэтому, проверим наш результат с помощью теоремы эквивалент ности. 7

В нашем случае, при 𝐿 = 𝐼, 𝜙(𝑥, 𝜉2* ) = 𝑓 𝑇 (𝑥)𝐷(𝜉2* )𝐿𝐷(𝜉2* )𝑓 (𝑥) = 𝑓 𝑇 (𝑥)𝐷2 (𝜉2* )𝑓 (𝑥). Для плана 𝜉2* дисперсионная матрица и её след равны: ⎛ ⎞ 2 0 −2 ⎜ ⎟ ⎜ ⎟ * 𝐷(𝜉2 ) = ⎜ 0 2 0 ⎟ , 𝑡𝑟 𝐷(𝜉2* ) = 8. ⎝ ⎠ −2 0 4 Тогда функция 𝜙(𝑥, 𝜉2* ) будет выглядеть как ⎛ ⎞2 2 0 −2 ⎜ ⎟ ⎜ ⎟ * 2 𝑇 𝜙(𝑥, 𝜉2 ) = (1, 𝑥, 𝑥 ) ⎜ 0 2 0 ⎟ (1, 𝑥, 𝑥2 ) = 8 + 20𝑥2 (𝑥2 − 1). ⎝ ⎠ −2 0 4 Учитывая, что 𝑥 ∈ [−1, 1], нетрудно видеть, что 𝑚𝑎𝑥 𝜙(𝑥, 𝜉2* ) = 8 = 𝑡𝑟 𝐷(𝜉2* ). При 𝑥∈X этом, он достигается в точках 𝑥𝑖 = {−1, 0, 1}. Следовательно, по теореме эквивалентности план 𝜉2* действительно является L-опти мальным. 1.2. Случай кубической модели Теперь рассмотрим чуть более сложную, полиномиальную модель с четырьмя па раметрами: 𝑓 (𝑥) = (1, 𝑥, 𝑥2 , 𝑥3 )𝑇 , матрица L и интервал планирования прежние. В плане будет 4 точки и также из соображений симметрии, подобно предыдущему пункту, мы будем искать его в виде ⎛ 𝜉3 = ⎝ −1 −𝑥 𝑥 1−𝜔 2 𝜔 2 𝜔 2 1 1−𝜔 2 ⎞ ⎠. Информационная матрица будет выглядеть так: ⎛ ⎞ 1 0 1 + 𝜔(𝑥2 − 1) 0 ⎜ ⎟ ⎜ ⎟ 2 4 ⎜ 0 1 + 𝜔(𝑥 − 1) 0 1 + 𝜔(𝑥 − 1)⎟ ⎜ ⎟, 𝑀 (𝜉3 ) = ⎜ ⎟ 2 4 ⎜1 + 𝜔(𝑥 − 1) ⎟ 0 1 + 𝜔(𝑥 − 1) 0 ⎝ ⎠ 0 1 + 𝜔(𝑥4 − 1) 0 1 + 𝜔(𝑥6 − 1) а дисперсионная 𝐷(𝜉3 ) = 𝑀 −1 (𝜉3 ): ⎛ 𝑏 ⎜ 𝑏−𝑎2 ⎜ ⎜ 0 𝐷(𝜉3 ) = ⎜ ⎜ 𝑎 ⎜ 𝑎2 −𝑏 ⎝ 0 0 𝑎 𝑎2 −𝑏 𝑐 𝑎𝑐−𝑏2 0 0 1 𝑏−𝑎2 𝑏 𝑏2 −𝑎𝑐 0 8 0 ⎞ ⎟ 𝑏 ⎟ ⎟ 𝑏2 −𝑎𝑐 ⎟ ⎟, 0 ⎟ ⎠ 𝑎 𝑎𝑐−𝑏2

где 𝑎 = 1 + 𝜔(𝑥2 − 1), 𝑏 = 1 + 𝜔(𝑥4 − 1), 𝑐 = 1 + 𝜔(𝑥6 − 1). След этой матрицы равен 𝑡𝑟 𝐷(𝜉3 ) = 𝜔(𝑥4 − 1) + 2 + 1 + 𝜔(𝑥4 − 1) − (1 + 𝜔(𝑥2 − 1))2 + 𝜔(𝑥6 + 𝑥2 − 2) + 2 (1.1) (1 + 𝜔(𝑥2 − 1))(1 + 𝜔(𝑥6 − 1)) − (1 + 𝜔(𝑥4 − 1))2 Чтобы найти inf этой функции, решим систему, полученную из необходимых усло вий экстремума: ⎧ ⎪ ⎨𝑡𝑟 𝐷(𝜉3 )′𝑥 = 0, ⎪ ⎩𝑡𝑟 𝐷(𝜉3 )′ = 0. 𝜔 𝑡𝑟 𝐷(𝜉3 )′𝑥 = 4(2𝑥6 𝜔 − 3𝑥2 𝜔 + 𝜔 + 2𝑥4 + 3𝑥2 − 1) = 0. (𝜔 − 1)𝜔𝑥3 (𝑥2 − 1)3 (1.2) 2𝑥6 𝜔 2 + 𝑥2 (4𝜔 − 2) − 2(𝜔 − 1)2 = 0, 𝑥2 (𝑥2 − 1)2 (𝜔 − 1)2 𝜔 2 (1.3) 𝑡𝑟 𝐷(𝜉3 )′𝜔 = Из уравнений (1.2), (1.3) получаем ограничения на 𝑥 и 𝜔: 𝜔 ̸= 0, 𝜔 ̸= 1, 𝑥 ̸= 0, 𝑥 ̸= ±1. (1.4) Из (1.2) выразим 𝜔: 𝜔= 1 − 3𝑥2 − 3𝑥4 . 2𝑥6 − 3𝑥2 + 1 (1.5) Теперь подставим выражение в (1.3), получив уравнение от 𝑥: 2𝑥 6 (︂ 1 − 3𝑥2 − 3𝑥4 2𝑥6 − 3𝑥2 + 1 )︂2 )︂ (︂ )︂2 (︂ 1 − 3𝑥2 − 3𝑥4 1 − 3𝑥2 − 3𝑥4 +𝑥 4· 6 −2 −2 − 1 = 0. 2𝑥 − 3𝑥2 + 1 2𝑥6 − 3𝑥2 + 1 2 Перепишем это уравнение в альтернативной форме: 2𝑥2 (5𝑥10 + 7𝑥8 − 2𝑥6 + 𝑥4 + 5𝑥2 − 1) = 0. (𝑥 − 1)(𝑥 + 1)(2𝑥4 + 2𝑥2 − 1)2 Решим уравнение численно, получим корни: 𝑥1 = −0.43955, 𝑥2 = 0, 𝑥3 = 0.43955 Заметим, согласно ограничениям (1.4), что 𝑥2 = 0 — посторонний корень. 9

Найдем 𝜔 из (1.5): 𝜔 = 0.7093. Таким образом, мы получаем план: ⎞ ⎛ −1 −0.43955 0.43955 1 ⎠. 𝜉3* = ⎝ 0.14535 0.35465 0.35465 0.14535 Также убедимся, что план является оптимальным, проверив его по теореме эк вивалентности (1). Вычисления будут аналогичны предыдущему пункту (для 𝑚 = 3), поэтому не будем приводить их здесь. 10

Глава 2 Нахождение локально D-оптимальных планов для нелинейных моделей 2.1. Обобщенная модель Михаэлиса-Ментен Рассмотрим модель с тремя параметрами и функцией регрессии 𝜂(𝑥, 𝜃) = 𝜃1 + 𝜃2 𝑥 , 𝑥 + 𝜃3 (2.1) где 𝜃1 , 𝜃2 , 𝜃3 ∈ R, 𝜃3 > 0. В качестве множества планирования эксперимента возьмем отрезок X = [0, 𝑑], где 𝑑 — некоторое заданное положительное число. Будем искать локально D-оптималь ный план. Сформулируем и докажем ряд предложений, которые помогут определить вид искомого плана. 2.1.1. Число точек плана Предложение 2.1.1. Число точек локально D-оптимального плана для модели (2.1) равно числу параметров 𝑚 = 3. Доказательство. Пусть план ⎛ 𝜉=⎝ 𝑥*3 ⎞ 𝜔1* , 𝜔2* , 𝜔3* ⎠ 𝑥*1 , 𝑥*2 , является D-оптимальным. Будем считать, без ограничения общности, что точки пере нумерованы в порядке возрастания: 0 ≤ 𝑥*1 < 𝑥*2 < 𝑥*3 ≤ 𝑑. (2.2) Для модели (2.1) вектор функций (︂ 𝑓 (𝑥) = 1, 𝑥 −𝜃2 𝑥 , 𝑥 + 𝜃3 (𝑥 + 𝜃3 )2 11 )︂T .

По теореме Кифера-Вольфовица [2, с. 53] имеем 𝑓 T (𝑥)𝑀 −1 (𝜉)𝑓 (𝑥) ≤ 3, 𝑥 ∈ [0, 𝑑], 𝑓 T (𝑥*𝑖 )𝑀 −1 (𝜉)𝑓 (𝑥*𝑖 ) = 3, 𝑖 = 1, 3. Обозначим 𝑔(𝑥) = 𝑓 T (𝑥*𝑖 )𝑀 −1 (𝜉)𝑓 (𝑥*𝑖 ) − 3. . Функция 𝑔(𝑥) представляет собой сумму дробей со знаменателями вида (𝑥+𝜃3 )𝑠 , 𝑠 = 0, 4. Приводя эти дроби к общему знаменателю, получим, что функция 𝑔(𝑥) выражается как 𝑔(𝑥) = 𝑃 (𝑥) , (𝑥 + 𝜃3 )4 (2.3) где 𝑃 (𝑥) — многочлен 4 степени. Заметим, что число точек в оптимальном плане 𝑛 ≥ 𝑚 = 3, иначе det 𝑀 (𝜉) = 0. Это следует из теоремы об информационных матрицах [2, с. 45]. С другой стороны, если 𝑛 > 3, то число нулей 𝑔(𝑥) с учетом кратности не менее, чем (𝑛 − 2) × 2 + 2 > 4 (будем иметь нули не менее, чем второй кратности в (𝑛 − 2) внутренних точках и 2 нуля не менее чем первой кратности в крайних точках). Получаем противоречие, поскольку 𝑔(𝑥) может иметь не более 4 нулей, с учетом их кратности, согласно его виду (2.3). Таким образом, приходим к выводу о том, что число опорных точек искомого D-оптимального плана 𝑛 = 3. 2.1.2. Весовые коэффициенты плана Предложение 2.1.2. Весовые коэффициенты локально D-оптимального плана для модели (2.1) равны во всех точках: 1 𝜔𝑖 = , 𝑖 = 1, 𝑛. 3 Доказательство. Информационная матрица * 𝑀 (𝜉 ) = 𝑚 ∑︁ 𝑓 (𝑥*𝑖 )𝑓 T (𝑥*𝑖 )𝜔𝑖 𝑖=1 12

может быть представлена [2, с. 64] в виде произведения трех матриц 𝑀 (𝜉) = 𝐹 T 𝑊 𝐹, где 𝐹 = (𝑓𝑗 (𝑥*𝑗 ))𝑛,𝑚 𝑖,𝑗=1 , a 𝑊 — диагональная матрица весовых коэффициентов 𝜔𝑖 , 𝑖 = 1, 𝑛. Из предложения 2.1.1 следует, что эти матрицы квадратные и одного порядка (𝑚 = 3). Запишем определитель матрицы 𝑀 (𝜉) в виде произведения: T det 𝑀 (𝜉) = det(𝐹 𝑊 𝐹 ) = (det 𝐹 ) 2 3 ∏︁ 𝜔𝑖* . (2.4) 𝑖=1 Отсюда, в силу неравенства Коши: 𝑛 ∏︁ (︂ ∑︀𝑛 𝑖=1 𝜔𝑖 ≤ 𝜔𝑖 𝑛 𝑖=1 )︂𝑛 , получаем, что максимум произведения (2.4) достигается, когда все три весовых коэф фициента равны 31 . Таким образом, мы перешли к задаче максимизации (det 𝐹 )2 . Определитель 𝐹 для плана, сосредоточенного в трех точках на отрезке [0, 𝑑] выглядит следующим образом: det 𝐹 = 𝐶 × (𝑥2 − 𝑥1 )(𝑥3 − 𝑥1 )(𝑥3 − 𝑥2 ) , (𝜃3 + 𝑥1 )2 (𝜃3 + 𝑥2 )2 (𝜃3 + 𝑥3 )2 (2.5) где 𝐶 — соответствующая константа. При опорных точках, упорядоченных по возрас танию (2.2), определитель det 𝐹 > 0, поэтому далее мы будем максимизировать det 𝐹 . 2.1.3. Опорные точки плана Для того, чтобы вычислить определитель матрицы 𝐹 , нужно знать, в каких точках сосредоточен оптимальный план. Сформулируем и докажем следующее предложение. Предложение 2.1.3. Две точки локально D-оптимального плана для модели (2.1) принадлежат границам отрезка X = [0, 𝑑], т.е. план имеет вид ⎛ ⎞ 0 𝑥 𝑑 ⎠, 𝜉* = ⎝ 1 3 1 3 1 3 где 𝑥 ∈ (0, 𝑑). Доказательство. Докажем последовательно два пункта. 13 (2.6)

а) В первую очередь покажем, что 𝑥1 = 0. Рассмотрим формулу (2.5) для определи теля 𝐹 , который соответствует плану, сосредоточенному в трех точках на отрезке [0, 𝑑]: (𝑥2 − 𝑥1 )(𝑥3 − 𝑥1 )(𝑥3 − 𝑥2 ) . (𝜃3 + 𝑥1 )2 (𝜃3 + 𝑥2 )2 (𝜃3 + 𝑥3 )2 Считаем, что точки отсортированы по возрастанию: det 𝐹 = 𝐶 × 0 ≤ 𝑥1 < 𝑥2 < 𝑥3 ≤ 𝑑. После замены 𝑥𝑖 → 𝑥𝑖 − 𝑥1 , 𝑖 = 1, 2, 3, значение числителя (2.5) не изменится, а знаменатель станет меньше. Значит, для максимизации det 𝐹 одну из точек нужно брать на левой границе отрезка, 𝑥1 = 0. б) Далее нужно показать, что 𝑥3 = 𝑑. Докажем от противного: пусть 𝑥3 < 𝑑. Тогда по теореме Кифера-Вольфовица, подобно тому, как мы делали в предложении 2.1.1, рассмотрим нули функции 𝑔(𝑥) = 𝑓 T (𝑥)𝑀 −1 (𝜉)𝑓 (𝑥) − 3. Эта функция в точке 𝑥1 = 0 имеет нуль 1 кратности, во внутренней точке 𝑥1 < 𝑥2 < 𝑥3 — нуль 2 кратности. Тогда при 𝑥3 < 𝑑 в данной точке будет также нуль 2 кратности и мы получим, что число нулей (с учетом их кратности) равно 5, что больше, чем допускает вид функции (2.3). Это противоречие. Значит, точка 𝑥3 должна быть на правой границе промежутка, то есть 𝑥3 = 𝑑. Таким образом, мы доказали, что из трех опорных точек оптимального плана две являются границами отрезка X = [0, 𝑑]. Значит, неизвестной остается только точка 0 < 𝑥2 < 𝑑. Напомним, что (︂ )︂ 𝑥 −𝜃2 𝑥 𝑓 (𝑥) = 1, , . 𝑥 + 𝜃3 (𝑥 + 𝜃3 )2 Тогда, пользуясь видом плана (2.6), получаем формулу для определителя: det 𝐹 = 𝜃2 𝑑𝑥(𝑑 − 𝑥) . (𝜃3 + 𝑑)2 (𝜃3 + 𝑥)2 (2.7) Для того, чтобы найти единственную неизвестную точку 𝑥, приравняем производ ную определителя к нулю, найдём точку экстремума det 𝐹 . Получим: 𝑥= 𝜃3 𝑑 . 2𝜃3 + 𝑑 14

2.2. Дробно-рациональная модель с четырьмя параметрами В этой части мы рассмотрим дробно-рациональную модель с четырьмя парамет рами: 𝜂(𝑥, 𝜃) = 𝜃1 𝜃3 + , 𝑥 + 𝜃2 𝑥 + 𝜃4 𝜃2 , 𝜃4 > 0. (2.8) Будем искать для нее локально D-оптимальный план на отрезке X = [0, 𝑑], где 𝑑 — заданный положительный параметр. Для этой модели сформулируем предложение, объединяющее предложения 2.1.1, 2.1.2, 2.1.3 предыдущего пункта: Предложение 2.2.1. Локально D-оптимальный план для модели (2.8) сосредоточен в 4 точках, с равными весами, т.е. имеет вид: ⎛ ⎞ * * * 0 𝑥2 𝑥3 𝑥4 ⎠, 𝜉* = ⎝ 1 4 1 4 1 4 1 4 (2.9) где 𝑥 ∈ (0, 𝑑). Доказательство строится аналогично доказательствам вышеприведенных предло жений, поэтому опустим его здесь. Для модели (2.8) вектор функций 𝑓𝑖 (𝑥) = 𝜕𝜂(𝑥,𝜃) , 𝑖 = 1, . . . , 4 выглядит так: 𝜃𝑖 (︂ )︂T −𝜃1 −𝜃3 1 1 , , , 𝑓 (𝑥) = . 𝑥 + 𝜃2 (𝑥 + 𝜃2 )2 𝑥 + 𝜃4 (𝑥 + 𝜃4 )2 Заметим, что локально D-оптимальный план не зависит от линейно входящих па раметров 𝜃1 , 𝜃3 , поскольку они не влияют на максимизацию определителя 𝑀 (𝜉). После этого, так же, как и в случае обобщенной модели Михаэлиса-Ментен максимизацию det 𝑀 (𝜉) мы сводим к максимизации определителя матрицы 𝐹 = (𝑓𝑗 (𝑥*𝑗 ))𝑛,𝑚 𝑖,𝑗=1 . Для на хождения неизвестных 𝑥2 , 𝑥3 и 𝑥4 будем максимизировать определитель det 𝐹 = 𝐶 × (𝑥2 + 𝜃2 )2 (𝑥 𝑥2 𝑥3 𝑥4 (𝑥4 − 𝑥3 )(𝑥4 − 𝑥2 )(𝑥3 − 𝑥2 ) , 2 2 2 2 2 2 + 𝜃4 ) (𝑥3 + 𝜃2 ) (𝑥3 + 𝜃4 ) (𝑥4 + 𝜃2 ) (𝑥4 + 𝜃4 ) (2.10) где 𝐶 — соответствующая константа. 2.2.1. Дифференциальное уравнение Идея поиска максимума определителя (2.10) c помощью дифференциальных урав нений развивает подход, предложенный Стильтьесом для максимизации выражений 15

подобного вида [12, с. 159]. Мы будем следовать методу нахождения D-оптимального плана, описанному в работе [6, с. 39]. Аналогичный функционально-алгебраический под ход также применялся для взвешенной полиномиальной регрессионной модели в работе [13]. Обозначим через 𝜓(𝑡) многочлен 𝜓(𝑡) = (𝑡 − 𝑥*2 )(𝑡 − 𝑥*3 )(𝑡 − 𝑥*4 ), (2.11) а коэффициенты этого многочлена обозначим через 𝜓0 , 𝜓1 , . . . , , 𝜓3 : 𝜓(𝑡) = 3 ∑︁ 𝜓𝑖 𝑡3−𝑖 , 𝜓0 = 1. (2.12) 𝑖=0 Производная функции (2.10) по 𝑥2 , 𝑥3 , 𝑥4 в силу необходимого условия экстремума должна обращаться в нуль при 𝑥𝑖 = 𝑥*𝑖 , 𝑖 = 2, 3, 4. Следовательно, получаем систему уравнений: ∑︁ 𝑄′(𝑥𝑖 ) 1 1 + −2 = 0, 𝑥𝑖 𝑥 𝑄(𝑥 𝑖 − 𝑥𝑗 𝑖) 𝑖̸=𝑗, 𝑖 = 2, 3, 4, (2.13) 2≤𝑗≤4 где 𝑄(𝑥) = (𝑥 + 𝜃2 )(𝑥 + 𝜃4 ). Используя формулу 1 1 ∑︁ 𝜓′′(𝑥𝑖 ) , = 2 𝑖̸=𝑗 𝑥𝑖 − 𝑥𝑗 𝜓′(𝑥𝑖 ) (2.14) которая была доказана в работе [6, c. 41], заметим, что функция ℎ(𝑥) = 𝜓′′(𝑥)𝑥𝑄(𝑥) + 2𝜓′(𝑥)[𝑄(𝑥) − 2𝑥𝑄′(𝑥)] (2.15) обращается в нуль при 𝑥 = 𝑥*2 , 𝑥*3 , 𝑥*4 . Следовательно, эта функция представляется в виде 𝜓(𝑥)𝜆(𝑥), где 𝜆(𝑥) = 𝜆0 𝑥 + 𝜆1 (так как ℎ(𝑥) является многочленом 4 степени, 𝜓(𝑥) — многочленом 3 степени). Итак, мы получили уравнение 𝜓′′(𝑥)𝑥𝑄(𝑥) + 2𝜓′(𝑥)[𝑄(𝑥) − 𝑥𝑄′(𝑥)] = 𝜆(𝑥)𝜓(𝑥), (2.16) где 𝑄(𝑥) = (𝑥 + 𝜃2 )(𝑥 + 𝜃4 ). Преобразуем его в алгебраическую форму и найдем решения — неизвестные опор ные точки локально D-оптимального плана. 16

2.2.2. Алгебраическое уравнение В работе [6, с. 39] было доказано, что дифференциальное уравнение вида (2.16) может быть представлено как алгебраическое: 𝜙T (𝑥)𝐴𝜓 = 𝜙T (𝑥)𝐶𝜆 𝜓, (2.17) где 𝜙(𝑥) = (𝑥4 , 𝑥3 , 𝑥2 , 𝑥, 1), 𝜓 = (1, 𝜓1 , 𝜓2 ) — коэффициенты разложения квадратного трехчлена 𝜓(𝑥), 𝐴 и 𝐶𝜆 — матрицы порядка 5 × 3. При решении данного уравнения мы будем рассматривать два случая для модели (2.8) и промежутка планирования X = [0, 𝑑], 𝑑 > 0 : а) d — достаточно большое б) d — малое значение Стоит уточнить, что мы понимаем под «большими» и «малыми» значениями d. Найдем D-оптимальный план для модели (2.8) и множества планирования X = [0, 𝑑]. Для достаточно большого значения 𝑑, наибольшая из опорных точек оптимального плана 𝑥*𝑛 < 𝑑. При уменьшении 𝑑, начиная с некоторого критического значения 𝑑* , наибольшая из точек плана 𝑥*𝑛 = 𝑑. Малыми значениями 𝑑 мы будем называть такие 𝑑 ≤ 𝑑* , для которых 𝑥*𝑛 = 𝑑. Для больших промежутков найдено решение, которое явно выражает опорные точки локально D-оптимального плана [6, с. 43]. Мы подробнее рассмотрим этот метод решения, с тем чтобы продолжить его на симметричный случай малых промежутков. Для малых значений 𝑑, когда не удается найти решение явно, будет предложен числен ный метод нахождения опорных точек, который будет реализован в виде алгоритма на языке 𝑅. 2.2.3. Решение уравнения для достаточно большого промежутка В случае достаточно большого 𝑑 опорные точки локально D-оптимального пла на 𝑥*𝑖 < 𝑑. При этом, уравнение (2.17) удается решить аналитически, получив явное выражение для 𝑥*2 , 𝑥*3 и 𝑥*4 . Покажем это. 𝑄(𝑥) = (𝑥 + 𝜃2 )(𝑥 + 𝜃4 ) = 𝑥2 + 𝑎𝑥 + 𝑏, 17

√ √ где 𝑎 = 𝜃2 + 𝜃4 , 𝑏 = 𝜃2 𝜃4 . Обозначим 𝑥˜ = 𝑥/ 𝑏, 𝑎 ˜ = 𝑎/ 𝑏. Тогда 𝑄(𝑥) = 𝑏(˜ 𝑥2 + 𝑎 ˜𝑥˜ + 1). Достаточно решить уравнение (2.17) для 𝑏 = 1. Решение для произвольного 𝑏 ∈ R √ √ находится с помощью обратной замены 𝑥 = 𝑥˜ 𝑏, 𝑎 = 𝑎 ˜ 𝑏. Знак волны далее будем опускать. В случае 𝑏 = 1 уравнение (2.16) принимает вид (6𝑥 + 2𝜓1 )𝑥(𝑥2 + 𝑎𝑥 + 1) + 2(3𝑥2 + 2𝜓1 𝑥 + 𝜓2 )(−3𝑥2 − 𝑎𝑥 + 1) = = (𝜆0 𝑥 + 𝜆1 )(𝑥3 + 𝜓1 𝑥2 + 𝜓2 𝑥 + 𝜓3 ). (2.18) Приводя подобные члены в левой части, получим −12𝑥4 − 10𝜓1 𝑥3 + (12 − 2𝑎𝜓1 − 6𝜓2 )𝑥2 + (6𝜓1 − 2𝑎)𝑥 + 2𝜓2 = ⎞ ⎛ −12 0 0 0 ⎟⎛ ⎞ ⎜ ⎟ ⎜ ⎜ 0 −10 0 0⎟ 1 ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ = (𝑥4 , 𝑥3 , 𝑥2 , 𝑥, 1) ⎜ 12 −2𝑎 −6 0⎟ ⎜𝜓1 ⎟ . ⎟⎝ ⎠ ⎜ ⎟ ⎜ ⎜ 0 6 −2𝑎 0⎟ 𝜓2 ⎠ ⎝ 0 0 2 0 (2.19) В правой части (2.18) имеем 𝜆(𝑥)𝜓(𝑥) = (𝑥4 , 𝑥3 , 𝑥2 , 𝑥, 1)𝐶𝜆 (1, 𝜓1 , 𝜓2 , 𝜓3 )T , (2.20) где ⎛ 1 ⎜ ⎜ ⎜0 ⎜ ⎜ 𝐶𝜆 = 𝜆0 ⎜0 ⎜ ⎜ ⎜0 ⎝ 0 ⎞ ⎛ 0 0 0 0 ⎜ ⎟ ⎜ ⎟ ⎜1 1 0 0⎟ ⎟ ⎜ ⎟ ⎜ 0 1 0⎟ + 𝜆1 ⎜0 ⎟ ⎜ ⎜ ⎟ ⎜0 0 0 1⎟ ⎝ ⎠ 0 0 0 0 ⎞ 0 0 0 ⎟ ⎟ 0 0 0⎟ ⎟ ⎟ 1 0 0⎟ . ⎟ ⎟ 0 1 0⎟ ⎠ 0 0 1 Отсюда сразу находим старший коэффициент 𝜆0 = −12. Для коэффициента 𝜆1 получаем уравнение: (𝐴 − 𝜆0 𝐸0 − 𝜆1 𝐸1 )𝜓 = 0, (2.21) где 𝐸0 = (𝐼4 𝑂1 ), 𝐸1 = (𝑂1 𝐼4 ). 𝐼4 — единичная матрица 4 × 4, 𝑂𝑖 — нулевая матрица 𝑖 × 4. 18

Таким образом, 𝜆1 является корнем уравнения det(𝐵−𝜆𝐼) = 0, где 𝐵 — квадратная матрица, получаемая из 𝐴 − 𝜆0 𝐸0 вычеркиванием первой строки (которая состоит из нулей). ⎞ ⎛ −𝜆 −2 0 0 ⎟ ⎜ ⎟ ⎜ ⎜ 12 −2𝑎 − 𝜆 6 0 ⎟ ⎟. 𝐵 − 𝜆𝐼 = ⎜ ⎟ ⎜ ⎜ 0 6 −2𝑎 − 𝜆 12 ⎟ ⎠ ⎝ 0 0 2 −𝜆 (2.22) Представляя определитель матрицы в виде суммы произведений определителей меньшего порядка, получаем, что он равен ⎞ ⎛ ⎛ ⎞ det ⎝ −𝜆 −2 12 −2𝑎 − 𝜆 ⎠ det ⎝ −2𝑎 − 𝜆 12 2 −𝜆 ⎛ ⎠ − det ⎝ −𝜆 0 12 ⎞ ⎠ det ⎝ 6 ⎞ ⎛ 6 12 ⎠= 0 −𝜆 = (𝜆(2𝑎 + 𝜆) − 24)2 − 36𝜆2 . (2.23) Получаем уравнение (𝜆2 + (2𝑎 − 6)𝜆 − 24)(𝜆2 + (2𝑎 + 6)𝜆 − 24) = 0. (2.24) Отсюда получаем, что возможные значения 𝜆1 равны −(𝑎 + 3) ± √︀ (𝑎 + 3)2 + 24, −(𝑎 − 3) ± √︀ (𝑎 + 3)2 + 24. Заметим, что вектор 𝜓 является решением уравнения ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 0 2 0 0 1 ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜𝜓 ⎟ ⎜12 −2𝑎 6 0 ⎟ ⎜𝜓1 ⎟ ⎜ ⎟ ⎜ ⎟ = 𝜆1 ⎜ 1 ⎟ . ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜𝜓2 ⎟ ⎜0 6 −2𝑎 12⎟ ⎜𝜓2 ⎟ ⎝ ⎠⎝ ⎠ ⎝ ⎠ 𝜓3 0 0 2 0 𝜓3 (2.25) Для нас представляют интерес только такие векторы 𝜓, которые соответствуют многочлену с положительными корнями. Раскрыв скобки в (2.11), получим 𝜓1 = −𝑥*2 − 𝑥*3 − 𝑥*4 , 𝜓2 = 𝑥*2 𝑥*3 + 𝑥*3 𝑥*4 , 𝜓3 = −𝑥*2 𝑥*3 𝑥*4 , 19

где 𝑥*𝑖 > 0, 𝑖 = 2, 3, 4. Отсюда получаем следующие ограничения: 𝜓1 < 0, 𝜓2 > 0, 𝜓3 < 0. (2.26) Теперь из первого уравнения (2.25) имеем 2𝜓1 = 𝜆1 , 𝜓1 = 𝜆1 . 2 Из второго уравнения: 12 − 2𝑎𝜓1 + 6𝜓2 = 𝜆1 𝜓1 . Подставив сюда 𝜓1 = 𝜆1 , 2 получаем: 12 − 𝑎𝜆1 + 6𝜓2 = 𝜆21 , 2 или, умножив на 2 и перенеся слагаемые в левую часть, 𝜆21 + 2𝑎𝜆1 − 12𝜓2 − 24 = 0. Так как 𝜆21 + 2𝑎𝜆1 ± 6𝜆1 − 24 = 0 из уравнения (2.24), а 𝜓2 должен быть положительным из ограничений (2.26), получаем 𝜆1 < 0, 𝜓2 = −𝜆1 /2. Из последнего уравнения получаем 2𝜓2 = 𝜆1 𝜓3 , 𝜓3 = −1. Так как √︀ (𝑎 + 3)2 + 24 > |𝑎+3|, а 𝜆1 должно быть отрицательным и являться решением уравнения 𝜆21 + 2𝑎𝜆1 + 6𝜆1 − 24 = 0, единственным возможным решением будет 𝜆1 = −(𝑎 + 3) − √︀ (𝑎 + 3)2 + 24. Заметим, что при 𝜃2 , 𝜃4 > 0, 𝑎 = 𝜃2 + 𝜃4 > 0. Следовательно, 𝜆1 < −3 − √ 33. Теперь получаем 𝜓(𝑥) = 𝑥3 + 𝜆1 2 𝜆1 𝜆1 𝑥 − 𝑥 − 1 = (𝑥 − 1)(𝑥2 + 𝑥(1 + ) + 1). 2 2 2 20

Точки оптимального плана являются корнями 𝜓(𝑥) и следовательно (︃ )︃ √︂ 𝜆 𝜆 1 1 1 1+ ± (1 + )2 − 4 , 𝑥*3 = 1. (2.27) 𝑥*2,4 = − 2 2 2 √ При 𝜃2 , 𝜃4 > 0, 𝜆1 < −3 − 33 < −6. Следовательно, выражение под корнем положительно и формула (2.27) задана корректно. Таким образом, для модели (2.8) в случае достаточно большого промежутка мы получили явное выражение для опор ных точек (2.27). Тогда, согласно утверждению 2.2.1 о виде плана, мы нашли локально D-оптимальный план. Отметим, что решение для больших промежутков было получе но в работе [6], но без детального изложения, которое нам здесь потребовалось для распространения метода на случай малых промежутков. 2.2.4. Решение уравнения для малого промежутка В случае малого промежутка (когда значение 𝑑 достаточно мало, 𝑥*4 = 𝑑) мат ричное уравнения (2.17) не допускает простого аналитического решения, поэтому явно выразить опорные точки, подобно пункту 2.2.3 не удается. Мы будем искать численное решение, пользуясь теоремой из работы [6, с. 41], которую сформулируем ниже. При 𝑥*4 = 𝑑 уравнение (2.16) примет вид 𝑥(𝑑 − 𝑥)𝑄(𝑥) + 𝜓′(𝑥)[𝑄(𝑥)(𝑑 − 2𝑥) − 2𝑥(𝑑 − 𝑥)𝑄′(𝑥)] = 𝜆(𝑥)𝜓(𝑥), где 𝑄(𝑥) = (𝑥 + 𝜃2 )(𝑥 + 𝜃4 ) = 𝑥2 + 𝑎𝑥 + 𝑏, 𝜓(𝑡) = (𝑡 − 𝑥*2 )(𝑡 − 𝑥*3 ), 𝜆(𝑥) = 2 ∑︀ (2.28) 𝜆𝑖 𝑥𝑖 . 𝑖=0 Так же, как в предыдущем пункте, будем решать уравнение только для 𝑏 = 1. После приведения подобных членов в левой части (2.28) получим: 3𝑥4 + (−𝑎 − 5𝑑 + 2𝜓1 )𝑥3 + (−𝑎𝑑 − 5 − 3𝑑𝜓1 ))𝑥2 + + (3𝑑 + 𝜓1 (−𝑎𝑑 − 2))𝑥 + 𝜓1 𝑑 = ⎛ ⎞ 3 0 0 ⎜ ⎟⎛ ⎞ ⎜ ⎟ ⎜−𝑎 − 5𝑑 2 0⎟ 1 ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ = (𝑥4 , 𝑥3 , 𝑥2 , 𝑥, 1) ⎜−𝑎𝑑 − 5 (2.29) −3𝑑 0⎟ ⎜𝜓1 ⎟ . ⎜ ⎟⎝ ⎠ ⎜ ⎟ ⎜ 3𝑑 −𝑎𝑑 − 2 0⎟ 𝜓2 ⎝ ⎠ 0 𝑑 0 В правой части (2.28) имеем: 𝜆(𝑥)𝜓(𝑥) = (𝑥4 , 𝑥3 , 𝑥2 , 𝑥, 1)𝐶𝜆 (1, 𝜓1 , 𝜓2 )T , 21

где ⎛ ⎛ 0 ⎟ ⎜ ⎟ ⎜ ⎜1 0⎟ ⎟ ⎜ ⎟ ⎜ 1⎟ + 𝜆1 ⎜0 ⎟ ⎜ ⎟ ⎜ ⎜0 0⎟ ⎠ ⎝ 0 0 1 0 0 ⎜ ⎜ ⎜0 ⎜ ⎜ 𝐶𝜆 = 𝜆 0 ⎜ 0 ⎜ ⎜ ⎜0 ⎝ 0 1 0 0 0 ⎞ ⎞ ⎛ 0 0 0 ⎟ ⎜ ⎟ ⎜ ⎜0 0 0⎟ ⎟ ⎜ ⎟ ⎜ 1 0⎟ + 𝜆2 ⎜1 ⎟ ⎜ ⎟ ⎜ ⎜0 0 1⎟ ⎠ ⎝ 0 0 0 ⎞ 0 0 ⎟ ⎟ 0 0⎟ ⎟ ⎟ 0 0⎟ . ⎟ ⎟ 1 0⎟ ⎠ 0 1 (2.30) Отсюда находим старший коэффициент 𝜆0 = 3. Для коэффициентов 𝜆1 и 𝜆2 полу чаем уравнение: (𝐴 − 𝜆0 𝐸0 − 𝜆1 𝐸1 − 𝜆2 𝐸2 )𝜓 = 0, (2.31) где 𝐸0 = (𝐼3 𝑂2 ), 𝐸1 = (𝑂1 𝐼3 𝑂1 ), 𝐸2 = (𝑂2 𝐼3 ). 𝐼3 — единичная матрица 3×3, 𝑂𝑖 — нулевая матрица 𝑖 × 3. Введем матрицу 𝐵 =𝐴− 2 ∑︁ 𝜆𝑖 𝐸𝑖 . (2.32) 𝑖=0 Элементы матрицы обозначим через 𝑏𝑖𝑗 . Пусть 𝐵(𝑖) матрица, составленная из (𝑖 + 1), (𝑖 + 2), (𝑖 + 3)-ей строк матрицы B. 𝐵(𝑖) = 𝐵(𝑖) (𝜆), 𝜆 = (𝜆1 , 𝜆2 ). Положим 𝑇𝑖 = 𝑇𝑖 (𝜆) = det 𝐵(𝑖) (𝜆) В работе [6, с. 41] была доказана теорема, которая помогает решить матричное уравнение (2.31). Для нашего случая соответствующая теорема будет звучать так: Теорема 2. Существует и притом единственное решение задачи 2 ∑︁ 𝑇𝑖2 (𝜆) → min2 (2.33) 𝜆∈R 𝑖=1 такое, что точки локально D-оптимального плана для функции регрессии (2.8), не совпадающие с 0 и 𝑑, являются корнями многочлена 𝜓(𝑥) = 𝜓0 𝑥2 + 𝜓1 𝑥 + 𝜓2 , (2.34) причем коэффициенты 𝜓0 , 𝜓1 , 𝜓2 могут быть вычислены при 𝜆, на котором достига ется минимум, по рекуррентным формулам 𝜓0 = 1, 𝜓𝑠+1 = − 𝑠 ∑︁ 𝑏𝑠+1,𝑗 𝜓𝑗 /𝑏𝑠+1,𝑠+1 , 𝑠 = 0, 1. 𝑗=0 22 (2.35)

Применив данную теорему, мы найдем значение 𝜆, на котором достигается ми нимум. После этого, подставив полученное значение в правую часть (2.28) получим уравнение относительно 𝑥, и далее, решив его, найдем значения двух неизвестных то чек локально D-оптимального плана, которые находятся внутри промежутка (0, 𝑑). Составлена программа на языке 𝑅, которая использует теорему 2 для нахождения неизвестных значений 𝑥*2 и 𝑥*3 . Опишем алгоритм работы программы. На вход программа принимает параметры модели — значения 𝑎 = 𝜃2 + 𝜃4 и 𝑑. По ним строится матрица 𝐴. Далее по этой матрице мы строим 𝑇𝑖 (𝜆) = det 𝐵(𝑖) (𝜆), 𝑖 = 1, 2 как функцию от входного параметра 𝜆. Тогда получаем сумму из утверждения теоремы 2 как функцию от вектора (𝜆1 , 𝜆2 ) ∈ R2 (считаем известным 𝜆0 = 3). Далее мы должны минимизировать данную функцию на плоскости. Для того, чтобы запустить функцию минимизации, необходимо задать некоторые ограничения на 𝜆1 и 𝜆2 , а также начальное приближение (𝜆01 , 𝜆02 ), в окрестности которого мы ищем экстремальное значение [14]. Удалось найти ограничения и метод поиска начального приближения, при которых функция минимизации находит значение вектора 𝜆, соответствующее единственному, согласно теореме 2, решению задачи. Ограничения на 𝜆1 и 𝜆2 были получены из условия на опорные точки: 𝑥*2 , 𝑥*3 ∈ (0, 𝑑). Тогда из 𝜓(𝑥) = (𝑥 − 𝑥*2 )(𝑥 − 𝑥*3 ) = 𝑥2 − (𝑥*2 + 𝑥*3 )𝑥 + 𝑥*2 𝑥*3 = 𝜓0 𝑥2 + 𝜓1 𝑥 + 𝜓2 получаем ограничения на коэффициенты 𝜓1 = −𝑥*2 − 𝑥*3 , 𝜓2 = 𝑥*2 𝑥*3 : −2𝑑 < 𝜓1 < 0, (2.36) 0 < 𝜓2 < 𝑑2 . Из матричного уравнения (2.31) ⎛ ⎞ 3 − 𝜆0 0 0 ⎜ ⎟⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜−𝑎 − 5𝑑 − 𝜆1 2 − 𝜆0 0 ⎟ 1 0 ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ (𝐴 − 𝜆0 𝐸0 − 𝜆1 𝐸1 − 𝜆2 𝐸2 )𝜓 = ⎜−𝑎𝑑 − 5 − 𝜆2 −3𝑑 − 𝜆1 − 𝜆0 0 ⎟ ⎜𝜓1 ⎟ = ⎜0⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠ ⎜ ⎟ ⎜ 3𝑑 −𝑎𝑑 − 2 − 𝜆2 −𝜆1 ⎟ 𝜓2 0 ⎝ ⎠ 0 𝑑 −𝜆2 выражаем из второй строки 𝜆1 = −𝑎 − 5𝑑 − 𝜓1 , 23 (2.37)

из третьей строки: 𝜆2 = 2𝑎𝑑 − 5 + 15𝑑2 + 𝜆1 (𝑎 + 8𝑑) + 𝜆21 − 3𝜓2 . (2.38) Тогда из (2.36) получаем ограничения на 𝜆1 , 𝜆2 : −𝑎 − 5𝑑 < 𝜆1 < −𝑎 − 3𝑑, (2.39) −10𝑎𝑑 − 38𝑑2 − 10 < 𝜆2 < 2𝑎𝑑 + 32𝑑2 − 10. Начальное приближение (𝜆01 , 𝜆02 ) мы будем строить по следующей пошаговой схеме. Сначала найдем решение для большого промежутка по формулам (2.27). Это решение соответствует всякому 𝑑 >= 𝑑* . Мы будем последовательно приближать значение 𝑑 к нашему 𝑑0 < 𝑑* , для которого ищем оптимальный план. Это означает, что мы будем делать следующее: зная оптимальный план для некоторого 𝑑−𝑘 , по его опорным точкам −𝑘 −𝑘 будем находить 𝜆−𝑘 как 1 , 𝜆2 по формулам (2.37), (2.38). Далее использовать вектор 𝜆 начальное приближение для решения следующей задачи с 𝑑 = 𝑑−𝑘+1 . Рис. 2.1. Иллюстрация схемы поиска начального приближения Таким образом, мы нашли ограничения (2.39) на область, по которой ведется ми нимизация, а также вектор начального приближения (𝜆01 , 𝜆02 ). Значит, мы находим точку экстремума функции из теоремы 2, 𝜆 = (𝜆1 , 𝜆2 ). При известном 𝜆 по формулам (2.35) находим коэффициенты многочлена 𝜓(𝑥), корни которого являются опорными точками 24

плана, не совпадающими с 0 и 𝑑. Находим 𝑥*2 и 𝑥*3 — опорные точки локально D-опти мального плана для малых значений 𝑑. Методы, реализующие вышеописанный алгоритм поиска численного решения за дачи приведены в приложении. Посмотрим на результат работы программы для опре деленного примера. Пример 1. Рассмотрим пример со значениями 𝜃2 = 51 , 𝜃4 = 5. Найдем локально D-опти мальный план по формулам (2.27) для случая достаточно большого промежутка: √︀ 𝜆1 = −(𝑎 + 3) − (𝑎 + 3)2 + 24 = −17.75196, √︂ 1 𝜆1 𝜆1 𝑥 ̃︀2 = − (1 + + (1 + )2 − 4) = 0.12908, 2 2 2 (2.40) 𝑥 ̃︀3 = 1, √︂ 𝜆1 𝜆1 1 − (1 + )2 − 4) = 7.74689. 𝑥 ̃︀4 = − (1 + 2 2 2 Таким образом, интерес для нас будет представлять задача поиска локально D-опти мального плана для модели 𝜂(𝑥, 𝜃) = 1 1 + , 𝑥 + 𝜃2 𝑥 + 𝜃4 на отрезке [0, 𝑑], где 𝑑 < 𝑥*4 = 7.746898. Возьмем значение 𝑑 = 7, запустим алгоритм, считая, что 𝑑 достаточно близко к 𝑑* и мы можем найти начальное приближение за один шаг (строим 𝜆0 по полученным в (2.40) значениям 𝑥 ̃︀2 , 𝑥 ̃︀3 ). Получим 𝜆01 = −39.07092, (2.41) 𝜆02 = −62.19083. Ограничения (2.39) будут выглядеть следующим образом: −40.2 < 𝜆1 < −26.2, (2.42) −2236 < 𝜆2 < 1630.8. Найдем точку минимума функции из теоремы 2, значение 𝜆: (𝜆1 , 𝜆2 ) = (−39.09320, −61.80159). По нему находим коэффициенты квадратного трехчлена 𝜓(𝑥): 𝜓0 = 1, 𝜓1 = −1.10679, 𝜓2 = 0.12536. 25 (2.43)

Находим корни 𝜓(𝑥) — это опорные точки плана 𝑥*2 , 𝑥*3 : 𝑥*2 = 0.12809, (2.44) 𝑥*3 = 0.97871. Таким образом, нашли решение, локально D-оптимальный план на отрезке X = [0, 𝑑] : ⎛ 𝜉1* = ⎝ 0 0.12809 0.97871 7 1 4 1 4 1 4 1 4 ⎞ ⎠. (2.45) Найденный план является оптимальным, поскольку если это не так, то должен ̂︀ Этому плану должен соот существовать другой, действительно оптимальный план 𝜉. ̂︀ подходящий под ограничения (2.39) и являющийся ветствовать некоторый вектор 𝜆, точкой минимума функции из теоремы 2. Это значит, что целевая функция имеет больше одного минимума, что противоречит единственности решения экстремальной задачи, о которой утверждается в теореме. Пример 2. Рассмотрим ещё один пример: 𝜃2 = 1 ,𝜃 10 4 = 10. Найдем локально D-опти мальный план по формулам (2.27) для случая достаточно большого промежутка: √︀ 𝜆1 = −(𝑎 + 3) − (𝑎 + 3)2 + 24 = −27.08606, √︂ 1 𝜆1 𝜆1 𝑥 ̃︀2 = − (1 + + (1 + )2 − 4) = 0.08029, 2 2 2 (2.46) 𝑥 ̃︀3 = 1, √︂ 1 𝜆1 𝜆1 𝑥 ̃︀4 = − (1 + − (1 + )2 − 4) = 12.46279. 2 2 2 Будем искать локально D-оптимальный план для модели 𝜂(𝑥, 𝜃) = 1 1 + , 𝑥 + 𝜃2 𝑥 + 𝜃4 на отрезке [0, 𝑑], 𝑑 = 10, который в данном случае будет являться малым промежутком, поскольку 𝑑 < 𝑥 ̃︀4 . Строим начальное приближение 𝜆0 по значениям 𝑥 ̃︀2 , 𝑥 ̃︀3 , получаем: 𝜆01 = −59.01976, (2.47) 𝜆02 = −137.58899. Ограничения (2.39) будут выглядеть следующим образом: −60.1 < 𝜆1 < −40.1, −4820 < 𝜆2 < 3392. 26 (2.48)

Найдем точку минимума функции из теоремы 2, значение 𝜆: (𝜆1 , 𝜆2 ) = (−59.06485, −136.31430). По нему находим коэффициенты квадратного трехчлена 𝜓(𝑥): 𝜓0 = 1, (2.49) 𝜓1 = −1.03515, 𝜓2 = 0.07594. Находим корни 𝜓(𝑥) — это опорные точки плана 𝑥*2 , 𝑥*3 : 𝑥*2 = 0.07946, 𝑥*3 (2.50) = 0.95569. Таким образом, нашли решение, локально D-оптимальный план на отрезке X = [0, 𝑑] : ⎛ 𝜉2* = ⎝ 0 0.07946 0.95569 10 1 4 1 4 1 4 27 1 4 ⎞ ⎠. (2.51)

Заключение Таким образом, мы построили L-оптимальный план для полиномиальной модели, т.е. для 𝑓 (𝑥) = (1, 𝑥, . . . , 𝑥𝑚−1 )T для 𝑚 = 2 : ⎞ ⎛ −1 0 1 ⎠, 𝜉2* = ⎝ 1 4 1 2 1 4 и для 𝑚 = 3 : ⎛ 𝜉3* = ⎝ −1 −0.439545 0.439545 0.145348 1 ⎞ ⎠. 0.354652 0.1453482 0.354652 Также мы рассмотрели нелинейную обобщенную модель Михаэлиса-Ментен с тре мя параметрами 𝜂(𝑥, 𝜃) = 𝜃1 + 𝜃2 𝑥 𝑥 + 𝜃3 и построили локально D-оптимальный план на отрезке X = [0, 𝑑] : ⎛ ⎞ 𝜃3 𝑑 0 𝑑 𝜉 = ⎝ 2𝜃3 +𝑑 ⎠ . 1 3 1 3 1 3 Далее, перейдя к рассмотрению дробно-рациональной модели с четырьмя пара метрами: 𝜂(𝑥, 𝜃) = 𝜃1 𝜃3 + , 𝑥 + 𝜃2 𝑥 + 𝜃4 мы подробно описали способ получения явной формулы для опорных точек в случае достаточно большого промежутка (когда d — достаточно большое число): √︂ 𝜆1 𝜆1 1 * 𝑥2,4 = − (1 + ± (1 + )2 − 4), 𝑥*3 = 1. 2 2 2 Наиболее интересным представляется новый результат относительно поиска ло кально D-оптимального плана для дробно-рациональной модели в случае малого про межутка (когда 𝑑 меньше некоторого критического значения 𝑑* ). Применив доказанную в работе [6, с. 41] теорему, мы получили метод нахождения численного решения для за дачи. Был написан алгоритм на языке 𝑅, реализующий данный метод. Частью алгоритма является минимизация специальной функции, описанной в утверждении теоремы 2. Для поиска минимума этой функции были получены необходи мые для минимизации ограничения (2.39) и схема построения начального приближения. 28

Результат работы алгоритма был показан на конкретных примерах случая малого промежутка: 1 𝜃2 = , 5 𝜃1 = 𝜃3 = 1, 𝜃4 = 5, X = [0, 7], для которого был найден локально D-оптимальный план: ⎛ ⎞ 0 0.12809 0.97871 7 ⎠, 𝜉1* = ⎝ 1 4 1 4 1 4 1 4 и 𝜃1 = 𝜃3 = 1, 𝜃2 = 1 , 10 𝜃4 = 10, X = [0, 10], для которого был найден локально D-оптимальный план: ⎞ ⎛ 0 0.07946 0.95569 10 ⎠. 𝜉2* = ⎝ 1 4 1 4 1 4 29 1 4

Литература 1. Fisher R. A. The Design of Experiments. –– New York : Hafner Pub. Co., 1966. –– 248 p. 2. Мелас В. Б., Шпилев П. В. Планирование и анализ для регрессионных моделей: Учеб. пособие. — Санкт-Петербург : Санкт-Петербургский государственный уни верситет, 2014. — 99 с. 3. Yu X.-Z., Gu J.-D. Differences in michaelis–menten kinetics for different cultivars of maize during cyanide removal // Ecotoxicology and environmental safety. –– 2007. –– Vol. 67. –– P. 254–259. 4. Clench H. How to make regional lists of butterflies: some thoughts // Journal of the Lepidopterists’ Society. –– 1979. –– Vol. 33. –– P. 216–231. 5. Rong C., Rappaport S. M. Relation between pulmonary clearance and particle burden: a michaelis-menten-like kinetic model // Occupational and Enviornmental Medicine. –– 1996. –– Vol. 53. –– P. 567–572. 6. Мелас В. Б. Локально оптимальные планы эксперимента. — Санкт-Петербург : Из дательство СПбГТУ, 1999. — 51 с. 7. Box G. E., Lucas H. Design of experiments in non-linear situations // Biometrika. –– 1959. –– Vol. 46. –– P. 77–90. 8. Katz D., Azen S., Schumitzky A. Bayesian approach to the analysis of nonlinear models: implementation and evaluation // Biometrics. –– 1981. –– Vol. 37. –– P. 137–142. 9. Sparrow P. E. Nitrogen response curves of spring barley // Journal of Agricultural Science. –– 1979. –– Vol. 92. –– P. 307–317. 10. Sparrow P. E. The comparison of five response curves for representing the relationship between the annual dry-matter yield of grass herbage and fertilizer nitrogen // Journal of Agricultural Science. –– 1979. –– Vol. 93. –– P. 513–520. 11. Ермаков С. М., Жиглявский А. А. Математическая теория оптимального экспери мента: Учеб. пособие. — Москва : Наука, 1987. — 320 с. 30

12. Сегё Г. Ортогональные многочлены. — Москва : Государственное издательство фи зико-математической литературы, 1962. — 500 с. 13. Chang F.-C. D-optimal designs for weighted polynomial regression – a functionalalgebraic approach // Statistica Sinica. –– 2005. –– Vol. 15. –– P. 153–163. 14. R: General-purpose optimization. –– https://stat.ethz.ch/R-manual/R-devel/ library/stats/html/optim.html. –– 2000. 31

Приложение А Реализация алгоритма 1 # количество строк и столбцов в матрице А 2 Nrow <- 5 3 Ncol <- 3 4 5 # получает матрицу А по значениям параметров a и d 6 getA <- function (a , d ) { 7 return ( matrix ( data = c (3 , 0 , 0 , 8 -a -5 *d , 2 , 0 , 9 -a *d -5 , -3 *d , 0 , 10 3 *d , -a *d -2 , 0 , 11 0 , d , 0) , nrow = Nrow , ncol = Ncol , byrow = TRUE )) 12 } 13 14 # получает матрицу E _ k 15 getE <- function ( k ) { 16 I <- matrix ( data = c (1 ,0 ,0 , 17 0 ,1 ,0 , 18 0 ,0 ,1) , nrow = Ncol , ncol = Ncol , byrow = TRUE ) 19 result <- I 20 if ( k > 0) { 21 for ( i in 1: k ) { 22 result <- rbind ( c (0 ,0 ,0) , result ) 23 24 } } 25 26 27 if ( Nrow - Ncol - k > 0) for ( i in 1:( Nrow - Ncol - k )) { 28 29 30 result <- rbind ( result , c (0 ,0 ,0)) } return ( result ) 31 } 32 33 # получает матрицу B по матрице A и значению lambda _ 0 34 # как функцию от lambda 32

35 getB <- function (A , l0 ) { 36 return ( function ( lambda ) { 37 38 ( A - l0 * getE (0) - lambda [1] * getE (1) - lambda [2] * getE (2)) }) 39 } 40 41 # получает T _ i ( определитель B _ i ) как функцию от i и lambda 42 getT <- function ( B ) { 43 return ( function (i , lambda ) { 44 45 det ( B ( lambda )[( i +1):( i +3) ,]) }) 46 } 47 48 # получает значение суммы из теоремы по данной функции T _ i 49 # как функцию от lambda 50 get _ sum _ T _ sqares <- function ( T ) { 51 return ( function ( lambda ) { 52 53 return (( T (1 , lambda )^2 + T (2 , lambda )^2)) }) 54 } 55 56 # получает psi по формулам из теоремы при известном 57 # значении psi _ 0 и матрице B 58 getPsi <- function ( psi0 , B ) { 59 psi <- NULL 60 psi [1] <- psi0 61 62 for ( s in 0:1) { 63 summa <- 0 64 for ( j in 0: s ) { 65 summa <- summa + B [ s +2 , j +1] * psi [ j +1] / B [ s +2 , s +2] 66 } 67 psi [ s +2] <- ( -1) * summa 68 } 69 70 return ( psi ) 71 } 72 73 # получает начальное приближение lambda по параметрам a и d 33

74 g et L a mb d aE s ti m a ti o n <- function (a , d ) { 75 lambda .1 <- -( a +3) - sqrt (( a +3)^2+24) 76 t2 <- ( -(1+ lambda .1 / 2) - sqrt ((1+ lambda .1 / 2)^2 -4)) / 2 77 t3 <- 1 78 t4 <- ( -(1+ lambda .1 / 2)+ sqrt ((1+ lambda .1 / 2)^2 -4)) / 2 79 psi .1 <- -t2 - t3 80 psi .2 <- t2 * t3 81 lambda _ 1 <- -a -5 * d - psi .1 82 lambda _ 2 <- 2 * a * d - 5 + 15 * d ^2 + lambda _ 1 * ( a +8 * d ) + lambda _ 1^2 - 3 * psi .2 83 lambda .0 <- c ( lambda _ 1 , lambda _ 2) 84 } 85 86 # решает задачу поиска опорных точек x _ 2 и x _ 3 используя теорему 87 solveByTheorem <- function (a , d ) { 88 lambda _ 0 <- 3 89 A <- getA (a , d ) 90 B <- getB (A , lambda _ 0) 91 T <- getT ( B ) 92 93 # получаем сумму из теоремы -- функцию от lambda для минимизации 94 sum _ T _ squares <- get _ sum _ T _ sqares ( T ) 95 96 # находим ограничения для lambda _ 1 и lambda _ 2 97 lower . lambda1 <- -a - 5 * d 98 upper . lambda1 <- -a - 3 * d 99 100 lower . lambda2 <- -10 * a * d - 38 * d ^2 -10 101 upper . lambda2 <- 2 * a * d + 32 * d ^2 - 10 102 103 # находим начальное приближение lambda 104 lambda . est <- g et L am b d aE s ti m a ti o n (a , d ) 105 106 # минимизируем целевую функцию 107 opt _ data <- optim ( lambda . est , sum _ T _ squares , method = "L - BFGS - B " , 108 lower = c ( lower . lambda1 , lower . lambda2 ) , 109 upper = c ( upper . lambda1 , upper . lambda2 )) 110 111 # получаем lambda , на котором достигается минимум 112 lambda <- opt _ data $ par 34

113 114 # считаем вектор psi по матрице B при найденном lambda 115 psi <- getPsi (1 , B ( lambda )) 116 117 # находим корни квадратного трехчлена psi ( x ) 118 roots <- polyroot ( rev ( psi )) 119 x2 <- as . numeric ( roots [1]) 120 x3 <- as . numeric ( roots [2]) 121 122 return ( c ( x2 , x3 )) 123 } 124 125 # пример 1 126 theta <- c (1 / 5 , 5) 127 a <- theta [1]+ theta [2] 128 d <- 7 129 130 # получаем решение для примера 1 131 solveByTheorem (a , d ) 132 133 # пример 2 134 theta <- c (1 / 10 , 10) 135 a <- theta [1]+ theta [2] 136 d <- 10 137 138 # получаем решение для примера 2 139 solveByTheorem (a , d ) 35

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв