CC BY
16ISSN 0868-5886 НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 3, c. 3-11
= ТЕОРЕТИЧЕСКИЕ ИССЛЕДОВАНИЯ
УДК 57.087+519.254+519.213.2
© А. Л. Буляница
МЕТОДЫ ОЦЕНИВАНИЯ ПАРАМЕТРОВ КРИВОЙ ЛОГИСТИЧЕСКОГО РОСТА. Ч. 1. ОПТИМИЗАЦИЯ УСЛОВИЙ ОЦЕНИВАНИЯ ПРИ НАЛИЧИИ АДДИТИВНОЙ СЛУЧАЙНОЙ ПОМЕХИ
Логистический рост является одной из традиционных математических моделей, применяемых для описания многих популяционных процессов. В частности, динамика количества ДНК при проведении полимеразной цепной реакции в реальном времени (ПЦР-РВ) адекватно описывается подобной моделью. Параметры кривой логистического роста служат основой для выявления необходимой аналитической информации об исходном количестве генетического материала. Объектом исследования является последовательность измерений, каждое из которых содержит в общем случае помимо полезной информативной составляющей случайную центрированную аддитивную помеху и систематическую составляющую в форме линейного тренда первого порядка с априорно неизвестными параметрами. Задачами данного раздела исследования является разработка алгоритмов эффективного оценивания параметров кривой логистического роста при наличии аддитивной центрированной ограниченной случайной помехи.
Кл. сл.: полимеразная цепная реакция (ПЦР), логистический рост, оценивание параметров, случайная помеха, стохастическая аппроксимация
ВВЕДЕНИЕ
Основам осуществления полимеразной цепной реакции (ПЦР) посвящено множество работ, в частности [1-6]. Ранее было доказано, что информативный сигнал в приборах, реализующих полиме-разную цепную реакцию в реальном масштабе времени, в первом приближении описывается кривой логистического роста. Даже предложенные альтернативные модели трехстадийной динамики (практическое отсутствие роста количества ДНК, экспоненциальный (геометрический) рост и выход на стадию насыщения) [7, 8] не противоречат на качественном уровне гипотезе о логистическом росте.
Одной из принятых форм этой кривой является
X = -
X *
1 + (B - 1К
(1)
Здесь Хп — величина сигнала на п-м шаге измерения (п = 0, 1,..., N где N обычно 30-40).
Подобная форма зависимости получается с помощью квантования по времени (дискретизации) решения задачи Коши
d x /d t = Xx (1 - x / x*); x(0) = x0.
(2)
Соответствующее дифференциальное уравнение первого порядка есть уравнение Бернулли (см., например, [9]) Оно решается в квадратурах отно-
сительно у = х ) при любых зависимостях от времени параметров Л = х* = х*(^). В част-
ности, в работе [10] рассмотрены формы решений уравнения (2) при линейных зависимостях параметров уравнения от времени.
Определяющими параметрами являются следующие.
1) Коэффициент а0, обычно находящийся в пределах 1.78-1.98 (при теоретическом пределе, равном 2), характеризующий эффективность проведения реакции. Этот коэффициент определяет средний прирост количества ДНК за 1 цикл.
2) Величина Х характеризует максимальный информативный сигнал, т. е. либо величину сигнала при максимальном количестве ДНК (при достижении режима насыщения), либо максимально допустимую величину сигнала, определяемую чувствительностью измерения, разрядностью устройств (в частности, АЦП) и другими факторами, связанными не с ПЦР, а с выбранным режимом измерения.
3) Параметр В численно равен отношению Х и Х0 — сигнала, соответствующего начальному (достаточно малому) количеству ДНК. Соответствующее отношение В обычно не менее 10 , но может достигать 10-10 .
Основные особенности рассматриваемых информативных сигналов.
п
1) Наличие аддитивной помехи, которую будем полагать ограниченной в диапазоне ±3 (возможно, от единиц до 10-20 отсчетов).
2) Максимальный сигнал Х полагаем превосходящим 1000 (возможно, до 3000-4000 отсчетов). Эта величина, в частности, определяется разрядностью используемого АЦП.
3) Сигнал, соответствующий начальному количеству ДНК Х0, как правило, не превосходит 1 отсчета (возможно, 10-3 и менее).
Главные методические проблемы анализа информативных сигналов состоят в следующем:
а) значимые параметры Х и а0 априорно неизвестны, а величина Х0 по существу и является искомой аналитической информацией;
б) при определенных условиях анализа кривая может не достигать стадии насыщения;
в) калибровка прибора обычно осуществляется на основе наиболее устойчивой комплексной характеристики сигнала — так называемого порогового цикла, т. е. стадии анализа, когда информативный сигнал обладает наибольшей линейностью. Аналитически положение порогового цикла может быть вычислено из соотношения п = 1п(В -1) / 1п(а0). При этом следует учесть, что
п необходим для определения В. Например, если при а0 = 2 пороговый цикл равен 20.00, то этому значению В при а0 = 1.8 будет соответствовать сильное смещение положения порогового цикла до 23.58. Таким образом, при осуществлении калибровки по положению порогового цикла п необходимо учитывать и эффективность ПЦР (т. е. величину а0).
Дальнейшие рассуждения касаются реализации алгоритмов последовательного оценивания базовых параметров информативного сигнала (а0, Х , а также Х0 или/и п).
НАЧАЛЬНЫЕ УСТАНОВКИ АЛГОРИТМА.
ОСНОВНЫЕ ПРИНЦИПЫ ОЦЕНИВАНИЯ ИНФОРМАТИВНЫХ ПАРАМЕТРОВ —
ПОСЛЕДОВАТЕЛЬНОСТЬ, РЕАЛИЗАЦИЯ И КАЧЕСТВО ОЦЕНИВАНИЯ
Полагаем, что информативный сигнал Хп содержит аддитивную ограниченную помеху С,п. (Для конкретизации считаем помеху равномерной и независимой).
Исследуем вторую разность информативного сигнала:
2Хп - Хп-1 - Хп+1 = 2С - С-1 - С+1 +
+X*(а0 - 2 +1/ а0)
* (1 - *)
(1 + * )(1 + * / а0)(1 + а0*)
(3)
Здесь * = (В -1) / ап . Исследование экстремумов
второй разности (и связанной с ним кривизны) зависимости (3) дает следующие результаты: * = 1 — точка нулевой кривизны (соответствует наиболее
линейному участку кривой и положению порого-
* \
вого цикла п = п ); в зависимости от а0 максимальная кривизна кривой реализуется при *, близких к 4 и % (точнее: при а0 = 2 эти значения 4.013 и 0.249, при а0 = 1.9 и менее — 3.87 и 0.26). Практически эти точки отстоят примерно на 2 отсчета от точки перегиба (или порогового цикла).
Ориентировочно оценим максимальную вторую разность в этих и соседних точках. Зададим а0 = 2, т. к. при достаточно близких значениях параметра результаты будут количественно схожи. Следует отметить, что замена * на 1/* означает движение по логистической кривой в другую сторону от точки перегиба. Соответствующие величины второй разности в этом случае будут совпадать по модулю и иметь противоположный знак. При * >1 эти отсчеты предшествуют точке перегиба, и вторая разность отрицательна; при * < 1 рассматриваются отсчеты после точки перегиба или порогового цикла.
Оценим величину второй разности при отсчетах, соответствующих * = 8, 4, 2, т. е. в трех предшествующих пороговому циклу отсчетах (полагаем при этом, что п — целая величина). При дробных п характерные величины второй разности в трех ближайших к пороговому циклу отсчетах будут соответствовать названным величинам *, деленным на а0 в степени {п} — дробная часть от п — (т. е.
на величину от 1 до 72).
В отсчетах * = 8, * = 4, * = 2 вторая разность — по абсолютной величине Х7 27.32, X/ 22.5 и Х7 30.0 соответственно. Суммарная аддитивная помеха при формировании второй разности (3) достоверно не превзойдет по величине 43. Тем самым гарантировано, что при условии 3 < XV 120 по крайней мере в 3 последовательных отсчетах, предшествующих точке перегиба, будет сохраняться знак второй разности. Далее при переходе через точку перегиба произойдет смена знака второй разности с отрицательного на положительный. В указанной ситуации с точностью примерно 0.5 отсчета может быть установлено положение точки перегиба (по интервалу смены знака второй разности).
Более подробно исследуем аддитивную помеху для второй разности. В первом приближении закон распределения помехи каждого измерения Хп полагаем равномерным. В принципе допустимо рассматривать любую гипотезу о распределении случайной составляющей помех при наложенных ограничениях о симметричности (и как следствие центрированности) и ограниченности. При независимости помех на (п - 1), п и (п + 1) отсчетах плотность распределения вероятностей (ПРВ) результирующей помехи получается методом сверт-
ки. Удобнее использовать аппарат характеристических функций (см., например, [11]). Переход к характеристическим функциям аналогичен преобразованию Лапласа с заменой аргумента р на мнимый аргумент (/^).
Для равномерной в диапазоне [-5, 5] помехи характеристическая функция имеет вид
©*( jw) =
1 ejwS _ q-jwS
25
jw
(4)
Результирующей помехе соответствует произведение характеристических функций вида (4), что приводит к выражению:
©( jw) =
1 (ejw 25 _ e_jw25 )(qjw5 _ Q-jw5 )2
1653
( jw)3
(5)
Здесь первый множитель соответствует слагаемому 2С,п, два других множителя в числителе соответствуют помехам С,п-\ и £П+1.
Обращение характеристической функции (5) аналогично переходу к оригиналу для преобразования Лапласа позволяет получить выражения для ПРВ в форме
Р( х) =
-[(х + 45)2^(х + 45) - 2(х + 25)2^(х + 25) +
+ 2( х - 25)2^( х - 25) + (х - 45)2^( х - 45)].
Здесь г)(х) — функция Хевисайда (функция единичного скачка). Обращение характеристической функции можно выполнить либо с использованием таблиц преобразования Лапласа [12], либо на основе теории вычетов. Соответственно для характерных интервалов явное выражение ПРВ примет вид:
325
(6)
x < _45 : p(x) = 0,
x е [_45,_25) : p(x) = (x + 45)2 / (3253), x е [_25;0] : p(x) = (852 _ x2)/(3253); Vx ^ p (_ x) = p(x).
График ПРВ представлен на рис. 1, а соответствующая функция распределения иллюстрируется рис. 2.
На основе (6) легко рассчитать вероятности попадания суммарной помехи в заданный диапазон значений: вероятность попадания в интервал [0, 15] равна 23/96, [15, 25] — 17/96, [25, 35] — 7/96 и [35, 45] — 1/96. Математическое ожидание равно нулю, дисперсия как сумма дисперсий 3 независимых слагаемых будет (4+1+1) 52 / 3. Таким образом видно, что подобная суммарная помеха "фокусируется" в окрестности x = 0.
0.25 Н
0.20-
0.15-
0.10-
0.05
0.00-
2 ! Я Л
x / 5
Рис. 1. ПРВ случайной составляющей второго разностного сигнала
1.0-1 фр 0.80.60.40.2-0.0-»—»—•
~I-'-г
2 / X 4
x/5
Рис. 2. Функция распределения (ФР) случайной составляющей второго разностного сигнала
Алгоритмы оценивания информативных параметров сигнала а0, X и В базируются на принципе построения постоянного сигнала (тренда нулевого порядка) на основе имеющихся измерений Хп при аддитивной случайной помехе С,п.
На первом этапе строится последовательность оценок величины а0 в форме
X ,(X . _X )
. n+ 2 v n+1_n '
Xn (Xn+ 2 _ Xn+1)
(7)
Указанная схема, учитывающая три последовательных измерения, начиная с отсчета п (далее трехточечная схема), обладает большей алгоритмической простотой, но меньшей точностью. В принципе возможно решение задачи оптимизации выбора базовых отсчетов п, п + 1 и п + 2 по крите-
рию минимума дисперсии ошибки определения а0. Соответствующее положение отсчетов должно удовлетворять условию * — 1.7. (Это положение точки оптимума достаточно слабо зависит от измеряемой величины а0). На практике следует учесть, что этот оптимум может реализоваться приближенно, поскольку п должно быть целым.
Более точной схемой для построения оценки величины а0 и как следствие а0 является 4-точечная схема с использованием сигнала в отсчетах п, п + 1, п + к и п + к + 1:
^.к _ n+k n+k+\
tfA —
■^n^n+1 (Xn + k + 1 - Xn
)
(8)
Таким образом, используются разделенные на к отсчетов две пары последовательных измерений.
Дисперсия ошибки оценивания (в линейном приближении) зависит от п и к (задача оптимизации по критерию минимума дисперсии решается поиском нетривиального минимума функции 2 переменных * и к). Положение оптимума также слабо зависит от величины а0: при а0 достаточно больших (1.90 и более) положение оптимума соответствует * - 4.5, к = 4; при менее эффективной ПЦР (а0 - 1.80) оптимальными будут условия * - 5 и к = 5. В принципе незначительное отклонение базовых отсчетов (п, п + 1, п + к и п + к + 1) от оптимальных слабо влияет на величину дисперсии ошибки оценивания а0.
Вторым этапом является определение максимального сигнала Х в соответствии с выражением
X* =-
ак -1
ак / X.^ -1/ X.
(9)
После вычисления Х* имеется возможность оценить 1п(В -1), далее вычислить В или/и положение порогового цикла п . Используются уже выбранные отсчеты со значениями величин сигнала Х„:
ln(B -1) = ln(X* / Xn -1) + nln(a0).
(10)
ИМИТАЦИОННОЕ МОДЕЛИРОВАНИЕ И ПРАКТИЧЕСКАЯ РЕАЛИЗАЦИЯ ПРОЦЕДУРЫ НАХОЖДЕНИЯ ПАРАМЕТРОВ КРИВЫХ ЛОГИСТИЧЕСКОГО РОСТА
В качестве генератора равномерно распределенной случайной величины использована программа RANDOM на языке С++. Генерация соответствующей случайной величины с диапазоном [-5, 5] выполнена как
28 (RANDOM(1001) /1000 -1/2) .
В круглых скобках показан способ генерации равномерно распределенной центрированной случайной величины в диапазоне [-0.5, 0.5].
Далее представлены результаты имитационного моделирования расчета параметров информативного сигнала. Заданы следующие параметры информативного сигнала и случайной помехи: а0 = = 1.920; X* = 3500; Х0 = 0.02; 3= 10. Этим параметрам соответствует п = 18.507; В = 174999.
По критерию смены знака второй разности выявлено положение порогового цикла на интервале 18-19 отсчетов. Выбрана 4-точечная схема (8) для оценивания а0: к = 4 и п = 13-18. Значению п = 13 соответствует * = 41.5 (исходя из предположения, что а0 = 1.9), т. е. этот отсчет примерно на 3-4 такта предшествует оптимальному режиму оценивания а0.
6 оценок а0 равны: 1.918, 1.896, 1.918, 1.906, 1.917, 1.927. Результирующая оценка 1.917. На основе этих же измерений получены, согласно (9), оценки X*: 3395, 3385, 3491, 3490, 3504, 3505, и итоговая оценка 3490.2, а также оценки (10) величины 1п(В -1): 12.029, 12.031, 12.042, 12.045, 12.038, 12.051. На базе результирующей оценки 1п(В -1) получено В = 169388 и как следствие — X0 = 0.0206; п = 18.497. Все оценки значения достаточно близки к заданным значениям.
Оценки параметров вычислялись на основе различных базовых измерений вблизи области оптимальных значений (п, к). В качестве результирующей оценки была использована медиана выборки 6 измерений, т. е. полусумма 2 центральных порядковых статистик. Как известно, она обладает робастностью и одновременно с этим учитывает и различие оценок, полученных на основе разных отсчетов. Далее будет рассмотрена альтернативная процедура получения результирующей оценки.
В целом наблюдается хорошая сходимость оценок параметров к истинным значениям. Условие приближенного "отслеживания" точки перегиба в форме 3 < X* /120 было выполнено во всех случаях более чем с двукратным запасом.
АНАЛИЗ ПОГРЕШНОСТЕЙ ОЦЕНИВАНИЯ ПАРАМЕТРОВ КРИВОЙ ЛОГИСТИЧЕСКОГО РОСТА
Общий вид оценки — аа0 = а0 +Да = а0(1 + еа) при Да и £а — абсолютной и относительной погрешностях. В первом приближении (приближении дифференциала первого порядка) для 3-точечной схемы (7) абсолютная погрешность имеет вид:
а
Да - (Кп^п + Кп+£п+1 + К+2С2),
где
К = -
-(* + 2 +1/*):
К = а0 +1
Кп+1 = —
п+2
-(* / а0 + 2 + а0 / *) .
(11)
—(* / а02 + 2 + а02 / *) 1
( X * )'
-3Ч2( К1 + кп+1 + к+2).
(12)
4 • 242.268 3
(3 / X * )2 - 323(3 / X * )2.
= ао
2 32&1-1-(К + К+1 + К2+к + К2+к+1), (13)
где
(X *)
К =—+ 2 +1/*) ,
п а0-Л '' К„+1 =—Ц (*/ а0 + 2 + а0/ *)
п 1 а0 -1 0 0
К„+к = (* / а0 + 2 + а0к / *) , 1
К
-(* / а0к+1 + 2 + а0к+1/ *) .
Правые множители формируют функцию 2 пе-
ременных (*, к), и возможен поиск минимума указанной функции. Точка минимума при заданном а0 искалась приближенным методом спирального (координатного) спуска. Как говорилось ранее, координаты этой точки слабо зависят от а0. При эффективно реализованной ПЦР (т. е. при а0 - 2) точка оптимума лежит вблизи * = 4.7, к = 4. В этом случае минимальная дисперсия ошибки оценива-
4 • 26.500,
Учитывая характер аддитивных случайных погрешностей, в частности их центрированность и независимость, можно сделать выводы:
а) оценка — несмещенная;
б) ее дисперсия равна
ния определяется как аа
3
(3 / X* )2
Здесь (г02 — дисперсия исходной помехи, приведенная к диапазону [-1, 1]. В случае равномерной помехи = 1/ 3 . Выражение в круглых скобках является функцией *, и возможны постановка и решение задачи поиска минимума соответствующей функции. Как говорилось ранее, имеется нетривиальное решение * - 1.7, слабо зависящее от а0 (* - 1.73 при а0 = 2, * - 1.68 при а0 = 1.9, * - 1.61 при а0 = 1.8). В частности, при а0 = 2 минимальная дисперсия ошибки оценивания
В общем случае при произвольных истинных значениях параметров кривой оптимуму может соответствовать нецелое значение п, и потребуется выбрать ближайший к точке оптимума отсчет п.
Схожий вид (и аналогичные свойства) имеет дисперсия ошибки определения а0 при 4-точечной схеме (8):
- 35.3(3 /X*) . Тем самым стандартное отклонение оценки (13) в случае 4-точечной оптимальной схемы будет приблизительно в 3 раза меньшим по сравнению с (12), т. е. аналогичной 3-точечной схемой.
Проведенное имитационное моделирование, включающее генерацию идеальной кривой логистического роста и аддитивной помехи с заданным законом распределения (100-250 различных реализаций), подтверждает указанные выше закономерности оценки а0. В частности, близкая к оптимальной 4-точечная схема при заданных условиях (а0 = 2.000; Х* = 3000; Х0 = 10; 3 = 5) дает следующие характеристики оценки а0: математическое ожидание 1.99939, среднеквадратичное отклонение 0.01023 (в соответствии с расчетной формулой (13), среднеквадратичное отклонение в указанных условиях должно составить 0.009907). Некоторое расхождение может быть объяснено неидеальностью характеристик датчика равномерной случайной величины. В частности, не в полной мере обеспечивается несмещенность случайной величины, имеются незначительные асимметрия и превышение коэффициента эксцесса.
При анализе ошибки оценивания Х*, связанной как с наличием аддитивной погрешности в базовых измерениях Xп и Xn+k, так и с погрешностью оценивания а0, следует учитывать, что погрешности могут не быть независимыми, если оценивание а0 и Х базируется на одних и тех же измерениях (на отсчетах п и п + к). Представляет интерес решение задачи о выборе точек для оценивания Х: что является более выгодным по точности — использование общих точек или различных? Решение этой задачи будет приведено далее.
В первом приближении приращение функции нескольких переменных можно вычислить на основе аналога дифференциала первого порядка:
8Х' А у
-Да ^ при
п+т г
ах* ах*
ДХ =-Да +-дх +-
аа0 0 ахп п ахп+т
0 п п+т
пользовании оценки X* =-
ат -1
/ Хп+т - 1/Хп
ис-
ана-
логичной (9). Здесь приращения ДХп и ДХп+т равны соответственно £п и £п+т. (Номера п и п + т могут не совпадать с аналогичными номерами отсчетов,
а
0
1
а
по которым осуществлялось оценивание а0). Задавая X* = (В -1) / аЩ, получим, что
ах*
дап
т х*г* ах*
ат -1 а0 аХп ат -1
1 „ ах*
(1 + о и
ах„
= —0— (1 + X* / ат ) . В случае несовпадения отсче-
а0т -1 '
тов, по которым осуществлялось оценивание а0 и X*, все три слагаемых в АХ* будут независимы, и соответствующие дисперсии оценок сложатся:
°2Х- =
1
(ат -1)
^ (х * )2 ()2 +
+а02т(1 +X* /а0т)252ст02 + (1 + X*)252ст0
У этой функции имеется тривиальный асимптотический минимум 0. Практически это означает, что наиболее точные измерения х выполняются при максимальных номерах отсчетов п и п + т. Таким образом, для нескольких последних измерений требуется вычислить X*; в случае его малости дисперсия ошибки оценивания Х" будет
2 0
примерно равна <уг~ —-
а 2т +1
(ат -1)2
5гст}, , т. е. лишь не-
значительно превзойдет дисперсию аддитивной помехи, равную 52ст2 .
Другим подходом является оценивание х на основе тех же измерений, что и при вычислении оценки а0, т. е. X = X, т = k для 4-точечной схемы (8). При этом для оценивания а0 и х (по (9)) использованы одни и те же отсчеты п, п + 1, п + k и п + k + 1 и
ах *
а0 -
1 (КпСп + Кп+1^п+1 + Кп+кСп+к + Кп+к+1^п+к+1 ),
где
а X
К* + 2 +1/X) - (1 +X),
п а0 - 1
Кп*+1 =--~л (X/ а0 + 2 + а0/X),
а0 -1
а X
К;+к =--0—(Г / а'к + 2 + а'к / X) + а0к + X,
а -1
а X
К+к+1 = -°- (X / а0к+1 + 2 + ак+1/ X). а0 - 1
Поскольку все четыре слагаемых являются независимыми, то их дисперсии складываются при вычислении дисперсии оценивания х. Наиболь-
шая по абсолютной величине ошибка оценивания х определяется полуразмахом помехи 5, а также абсолютными значениями коэффициентов К . Если значение X находится достаточно близко к 4 (т. е. вблизи оптимума для вычисления а0), то максимальная погрешность оценивания х теоретически может составить около 125, т. е. огромную величину (хотя и вероятность этого ничтожно мала).
Вывод. Многократно более эффективной является следующая реализация алгоритма:
- на первом этапе при режимах, близких к оптимальным, производится оценивание параметра
а0;
- на втором этапе по набору измерений, достаточно близких к участку насыщения (т. е. по заключительным имеющимся отсчетам кривой логистического роста, соответствующим достаточно малым X*), производится, согласно (9), оценивание другого параметра х*. Эффективность подобного подхода при оценивании х иллюстрируется далее.
Пусть, для примера, сгенерирована последовательность 40 отсчетов с параметрами а0 = 1.910; х = 3000; х0 = 0.50 и добавлена равномерная центрированная помеха с 5 = 10. Далее по (8) с к = 4 и п = 10-15 проведено оценивание а0. После этого выполнен расчет х по тем же отсчетам. Получены следующие оценки: 3063, 2997, 2982, 3002, 3001, 3009 (стандартное отклонение ст» 29.6). Для сравнения по измерениям п = 29-34 последовательность таких же оценок (9) для х имеет вид: 2990, 3010, 2994, 3006, 3004, 3004 (ст » 7.8). Таким образом, оценивание х по более поздним отсчетам является существенно более эффективным, исходя из величины среднеквадратичного отклонения ошибки оценивания. Аналогичная тенденция наблюдается и при иных параметрах а0, х , х0 и диапазоне помехи.
Заключительным этапом является оценивание параметра В и на его основе величины порогового цикла п как точки перегиба кривой логистического роста. Используется формула (10) и строятся выборки значений F = 1п(В -1) = 1п(X * / хп -1) + +п 1п(а0). Тогда погрешность оценивания величины F в приближении первого дифференциала бу-
дF дF дF дет АF =-Аа0 +--Ах * +--Ах . Полагаем,
да0 0 дх* дх п
0 п
что предшествующие параметры были оценены ранее в оптимальных условиях и, следовательно, по различным отсчетам. Тем самым соответствующие погрешности оценивания были независимы, и дисперсии могут складываться.
Частные производные равны соответственно
да
= п / ап
дF
1
дх * х * - х.
_дF дх"
2
1
и
X *
X •/ X -1 XI
X (X * - X,,)
. Полагая, что
измерение Xn = X* /(1 +1'), где t' не совпадает ни со значениями t (при условиях оптимального вычисления а0), ни со значениями t (при оптимальном вычислении X) получим следующее.
а) Дисперсии а0 и X при оптимальных режимах
оценивания равны соответственно 1065 2о"0 / (X* )2
и д2&02 (т. к. дисперсия ошибки оценивания X* в оптимальных условиях оценивания лишь незначительно превосходит дисперсию аддитивной помехи любого из измерений).
б) В приведенных выше обозначениях
1 1+t' дF
дF
дX * X* t' ' - ln(t')) /1п(а0).
дX,,
X *t'
и п = (1п(В -1) -
мерно равна
(X * Г
t'
(1 +1 ')2 п2 + ^^ +106 п_
t '2 а2
212
ао21п2(ао)
1n(t') - 4/1 '2 - 6/1'+ 4t'+ 2t '2 =
2121п(В -1)
ао21п2(ао) '
Наиболее удобной представляется робастная модификация Фабиана—Цыпкина. Основанием для использования указанного алгоритма служит модель сигнала (выборка оценок, представляющая собой постоянную величину — истинное значение с аддитивной помехой с априорно неизвестным законом распределения).
Общие принципы организации алгоритма оценивания достаточно детально рассмотрены в работах [13, 14]. Выполняется рекурсивное построение оценок по правилу:
= сп - 3п -¥(сп - хп+:).
(15)
При этом \у( г) =
-1, для 2 < -А,
в) Тогда дисперсия ошибки оценивания F при" (1 +1 ')4
В принципе возможен поиск минимума дисперсии по параметру t', который, безусловно, зависит от истинного значения В. Необходимое условие экстремума при этих данных сведется к поиску положительного решения трансцендентного уравнения
(14)
В случае нахождения такого корня t' > 1, соответствующая вторая производная, безусловно, положительна, что свидетельствует об экстремуме типа минимума.
В частности, при В « 1000 и а0 = 2, оптимум (14) находится близко от t' = 15 ; соответствующее стандартное отклонение оценки логарифма (т. е. функции Р) примерно равно 20.55/X* (при 5 = 10, X* = 3000 отклонение логарифма на соответствующую величину 0.0683 приведет к отклонению оценки В в пределах 936-1071).
АЛГОРИТМ ПОСТРОЕНИЯ РЕЗУЛЬТИРУЮЩЕЙ ОЦЕНКИ ИНФОРМАТИВНЫХ ПАРАМЕТРОВ НА ОСНОВЕ КОРОТКОЙ ВЫБОРКИ
Ранее в качестве такого алгоритма предлагалось использование медианы по четному числу отсчетов. Практически аналогичного результата можно достичь при применении алгоритма стохастической аппроксимации Роббинса—Монро.
0, < А,
1, 2 > А;
хп+=с*+Е>п+1 — измерение; Е>п+\ — помеха; сп, сп+\ — оценки величины с на п-м и п+1-м шагах оценивания; / — параметр алгоритма; 2А — величина зоны нечувствительности. В предельном случае (15) при А = 0 получается сигнатурный алгоритм Я.З. Цыпкина, т. к. ¥(г) = sign(г). Начальное приближение с1 оценок а0, X и X0 может совпадать с первой из полученных оценок. При оценивании а0 в качестве первого приближения можно выбрать и априорную оценку 1.90, соответствующую средней эффективности ПЦР.
Степень близости оценки к истинному значению можно определить в соответствии с [13, 14] на основе анализа знака "поправок", т. е. при наличии немонотонности последовательности оценок. В противном случае при монотонном увеличении (уменьшении) оценок достаточно вероятно расхождение оценки и истинного значения параметра. При оценивании X и X0 априорная информации о возможном значении этих параметров отсутствует, и в качестве первого приближения следует брать первую оценку величины.
Параметр / в (15) соизмерим с дисперсией соответствующей оценки: при оценивании X сопоставим с величиной 5 (т. е. 10-15 квантов), при оценивании а0 дисперсия оценки близка к 65 / X (т. е. около 0.03). Соответствующий выбор / требует использования значений 0.05-0.10. Наконец, при оценивании величины Р = 1п(В - 1) стандартное отклонение примерно равно 205 / X, т. е. около 0.1. Тем самым следует выбрать /3 порядка 0.2. Выбранная величина зоны нечувствительности 2А должна быть сопоставима с требуемой точностью определения параметра. С практической точки зрения достаточным является оценивание а0 с точностью до 0.01, X с точностью до 1 или даже нескольких квантов, точность при оценивании Р определяет с масштабным множителем 1/ 1п(а0) « 1.4^1.7 погрешность нахождения положения порогового цикла п.
с
Выборки измерений параметров и оценок, полученных по алгоритму (15)
p A Параметр * Измерения / оценки Итоговая оценка
0.03 0.001 а0 1.918 1.896 1.918 1.906 1.917 1.927 1.917
1.918 1.888 1.903 1.903 1.911 1.917
10 0.5 Х* 3492 3511 3494 3506 3504 3504 3504.8
3492 3502 3497 3500.3 3502.8 3504.8
0.5 0.01 F 12.029 12.031 12.042 12.045 12.038 12.051 12.05
12.03 12.03 12.08 12.05 12.03 12.05
Измерения в верхней строке, оценки — в нижней.
Представляется достаточным нахождение п* с точностью до 0.2-0.3 цикла (отсчета). Тогда получаем результирующие оценки параметров на базе процедуры стохастической аппроксимации (в таблице). По правилу знаков "поправок", эти оценки близки к истинному значению. Исходные оценки, приведенные в таблице, получены при условиях, близких к оптимальным. Результирующие оценки: ао = 1.917, = 3504.8, F = 12.05, из чего следует, что В » 171100, Х0 » 0.0205 и п » 18.52. Они весьма близки к исходным данным, приведенным в разделе "Имитационное моделирование...".
ЗАКЛЮЧЕНИЕ
В работе исследованы различные алгоритмы оценивания параметров кривой логистического роста. Объединяющим принципом является переход к построению оценок в форме постоянного сигнала (тренда нулевого порядка) с аддитивной случайной помехой. Были проанализированы возможные алгоритмы построения подобных оценок, описаны законы формирования случайной погрешности оценки с доказательством центрированности соответствующих плотностей распределения вероятностей, а также решены задачи оптимизации выбора базовых отсчетов по критерию минимизации дисперсии ошибок оценивания. Учет влияния систематической составляющей сигнала планируется осуществить в ч. 2 статьи.
Данная работа выполнялась при финансовой поддержке в рамках программ и проектов: Аналитической ведомственной целевой программы "Развитие научного потенциала высшей школы" — проект "Исследования и диагностика клеточных структур: новые методические подходы и инстру-
ментальные решения на основе сканирующей зон-довой микроскопии и микрочиповых технологий" (№ 4247); Программы фундаментальных исследований Президиума РАН № 20 "Создание и совершенствование методов химического анализа и исследования структуры веществ и материалов" — проект "Микрочиповые аналитические системы для метода молекулярных колоний"; Программы фундаментальных исследований Президиума РАН на 2009 г. № 27 "Основы фундаментальных исследований нанотехнологий и наноматериалов".
СПИСОК ЛИТЕРАТУРЫ
1. Arnheim N. Polymerase Chain Reaction Strategy // Annual review of biochemistry. 1992. V. 61. P. 131-156.
2. Erlich H.A., Gelfand D., Sninsky J.J. Recent Advances in the Polymerase Chain Reaction // Sciences. 1991. V. 252, N 5013. P. 1643-1651.
3. Higuchi R., Fockler C., Dollinger C., Watson R. Kinetic PCR Analysis: Real-Time Monitoring of DNA Amplification Reactions // Biotechnology. N. Y., Sept. 1993. V. 11, N 9. P. 1026-1030.
4. Mullis K.B., Faloona F.A. Specific Synthesis of DNA in Vitrovia a Polymerase-Catalized Chain Reaction // Method in Enzymology. 1987. V. 155. P. 335-350.
5. White T.J., Madej R., Persing D.H. The Polyme-rase Chain Reaction: Clinical Applications // Advances in Clinical Chemistry. 1992. V. 29. P.161-196.
6. Livak K.J., Flood S.A., Marmaro J. et al. Oligonucleotides with Fluorescent Dyes at Opposite Ends Provide a Quenched Probe System Useful for Detecting PCR Product and Nucleic Acid Hybridization // PCR Methods and Applic. 1995. N 4. P.357-362.
7. Чернышев А.В., Друца В.Л. Проблемы создания оборудования для медицинской ПЦР-диагностики // Биомедицинские технологии и радиоэлектроника. 2004. № 12. С. 18.
8. Чернышев А.В. Создание теории рабочих процессов, методов расчета и разработка оборудования для ПЦР-диагностики. Автореф. дисс. ... доктора техн. наук. М.: МГТУ (МВТУ им. Н.Э. Баумана), 2006. 38 с.
9. Камке Э. Справочник по обыкновенным дифференциальным уравнениям. Пер. с нем. М.: Наука, 1976. С. 38-39.
10. Архипов Д.Б., Буляница А.Л., Щербаков А.П. Моделирование тенденций аналитического приборостроения логистическим уравнением с коэффициентами, зависящими от времени // Сб. "Тенденции развития аналитического приборостроения". СПб.: Русская классика, 2008. С. 20-21.
11. Горяинов В.Т., Журавлев А.Г., Тихонов В.И. Статистическая радиотехника: Примеры и задачи. Учеб. пособие для вузов. М.: Сов. радио, 1980. С. 50.
12. Первозванский А.А. Курс теории автоматического управления. Учеб. пособие. М.: Наука, 1986. 616 с.
13. Буляница А.Л., Курочкин В.Е. Исследование свойств и программно-аппаратная реализация алгоритма стохастической аппроксимации в модификации Я.З. Цыпкина // Научное приборостроение. 2002. Т. 12, № 2. С. 30-49.
14. Буляница А.Л., Курочкин В.Е., Бурылов Д.А. Реализация процедуры оценивания постоянного сигнала на основе метода стохастической аппроксимации в модификации Я.З. Цыпкина // Радиотехника и электроника. 2002. Т. 47, № 3. С. 343-346.
Институт аналитического приборостроения РАН, Санкт-Петербург
Материал поступил в редакцию 24.04.2009.
METHODS FOR LOGISTIC GROWTH CURVE PARAMETERS ESTIMATION.
PART. 1. OPTIMIZATION OF ESTIMATION CONDITIONS AT THE PRESENCE OF ADDITIVE RANDOM ERROR
A. L. Bulyanitsa
Institute for Analytical Instrumentation RAS, Saint-Petersburg
Logistic growth is one of the traditional mathematical models, used for describing many population processes in particular, the dynamics of DNA quantity in polymerase chain reaction in real time (PCR-RT) is adequately described by a similar model. The parameters of the logistic growth curve serve the basis for the development of the necessary analytical information about the initial quantity of genetic material. The subject of the study is the sequence of measurements, each of which, in a general case, besides useful informative component contains random centered additive error and systematic component as a linear trend of the first order with a priori unknown parameters. The aim of the present part of the study is the development of the algorithms of the effective parameters estimation of logistic growth in the presence of the additive centered limited random error.
Keywords: polymerase chain reaction (PCR), logistic growth, estimation of parameters, random error, stochastic approximation
