CC BY
50
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 - 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Параметрическая идентификация модели взаимодействующих
клеточных популяций на основе байсовского подхода
# 11, ноябрь 2012
DOI: 10.7463/1112.0490900
Виноградова М. С.
УДК 51.76: 517.9: 57.085.23
Россия, МГТУ им. Н.Э. Баумана mathmod@bmstu.ru
1. Введение
Клеточная трансплантология — новое перспективное и быстро развивающееся направление в современной лечебной практике. Предназначенные для пересадки клетки получают в нужном количестве путем их размножения в лабораторных условиях (in vitro), поскольку исходного количества клеток, взятых из организма пациента, недостаточно.
В процессе размножения в силу естественной изменчивости в культуре клеток могут появиться клетки с хромосомными мутациями (аномальные клетки), причем клоны аномальных клеток могут обладать селективным преимуществом по сравнению с нормальными клетками.
При трансплантации такого материала имеется риск возникновения канцерогенеза (рака) у пациента. Выращивание трансплантата является дорогой и длительной процедурой. При этом выявление факта преобладания аномальных клеток на поздних стадиях выращивания ведет к необходимости повторения процесса, что что связано со значительным увеличением материальных и временных затрат и, что самое важное, к увеличению времени ожидания операции, а задержка операции трансплантации может нести опасность для жизни пациента.
Непосредственное слежение за составом клеточной популяции лабораторными методами трудоемко, а его результаты не всегда достоверны. Для отбра-
ковки культуры клеток на ранних стадиях развития предлагается использовать современные информационные технологии, основанные на методах математического моделирования, методах параметрической идентификации моделей на ограниченных выборках экспериментальных данных о значениях вектора состояния изучаемой системы и статистических методах.
В работах [11, 12] предложена математическая модель развития клеточной популяции, состоящую из двух видов стволовых клеток человека: нормальных и аномальных клеток.
Основной проблемой при параметрической идентификации математической модели является малость массивов измерений ее переменных состояния в различные моменты времени.
В работе обоснован алгоритм решения задачи параметрической идентификации, и на основе байесовского подхода получены плотности распределения вероятности параметров модели.
2. Задача оценивания параметров модели
Рассмотрим задачу параметрической идентификации математической модели, предложенной в работах [11, 12]. Напомним, что в указанной математической модели динамика сумммарной численности популяции нормальных клеток полностью определяется соотношениями
X((п0 + 1)т0) = а0Х(п0т0) , п0 е {0,1, 2,...};
X (0) = Х0 ,
где
а0 = (1 - М0)(1 - А0) + 2М0(1 - А0)(1 - 70), (2)
а динамика сумммарной численности популяции аномальных клеток определяется соотношениями
У((п1 + 1)т= а1У(п1т+ Ь1Х(п1т1) , п1 е {0,1, 2,...};
(3)
у (0) = о,
где
а1 = (1 - М 1)(1 - А1) + 2М 1(1 - А1), Ь1 = 2М0(1 - А0)70.
Здесь Х0 — начальная численность нормальных клеток; А0, А1 — доля нормальных и аномальных клеток, погибающих на временных интервалах длительности т0 и т0 соответственно; 70 — доля нормальных клеток, каждая из которых в процессе деления на временном интервале длительности т0 становится аномальной; М0, М0 — доля нормальных и аномальных клеток, разделившихся на временных интервалах длительности т0 и т0 соответ-
Заметим, что непосредственно оценить значения параметров указанной модели не удается, поскольку по результатам экспериментов можно оценить только количество нормальных и аномальных клеток по различным пассажам [1, 2, 3]. Соответствующие измерения количества нормальных и аномальных клеток можно считать проведенными через равные промежутки времени.
Для решения задачи параметрической идентификации указанной математической модели на ограниченных выборках экспериментальных данных воспользуемся методами, предложенными в работе [8].
Приведем общую схему решения задачи оценивания параметров исходной математической модели, для чего рассмотрим дискретную динамическую систему
где V ) — п-мерный вектор состояния системы, А — матрица типа п х п параметров системы, которая полагается неизвестной. Пусть в дискретные моменты времени tk, к = 1,..., 2Х, известны значения {, к = 1, 2Х} наблюдений вектора состояния системы (5).
Сделаем следующие допущения:
1) последовательность tk, к = 1, 2Х, возрастающая и
ственно.
V (^+1) = AV (^),
(5)
= t2k - t2k-l = еопв1;
2) наблюдения являются прямыми и содержат случайные ошибки наблюдения {ек, к = 1, }:
т = V (1к) + ек, (6)
т.е. тк = тк(и) и ек = ек(и) — случайные величины;
3) ошибки наблюдения ек представляют собой независимые случайные векторы, распределенные по п-мерному нормальному закону с нулевом математическим ожиданием и неизвестной положительно определенной ковариационной матрицей
Под задачей оценивания будем понимать задачу нахождение оценки матрицы А неизвестных параметров математической модели (5) по имеющимся данным наблюдений {тк, к = 1, }.
При дальнейших исследованиях будем опираться на следующий результат.
Теорема 1 (см. [8]). Задача линейного оценивания неизвестных параметров исходной дискретной линейной модели (5) эволюционного процесса на ограниченных выборках экспериментальных данных о значениях векторов состояния (6) и внешних воздействий в смысле совпадения функций правдоподобия эквивалентна соответствующей ей задаче линейного оценивания с ошибками в инструментальных переменных, записанной на максимально возможной системе непересекающихся шаблонов разностного уравнения математической модели изучаемого эволюционного процесса.
Для построения непересекающихся шаблонов сформируем из множества {тк, к = 1, } множество четных и множество нечетных узлов. Тогда задача линейного оценивания матрицы А неизвестных параметров модели (5) по известным значениям наблюдений вектора состояния (6) в смысле совпадения функций правдоподобия будет эквивалентна задаче линейного оценивания с ошибками в инструментальных переменных:
т2к = AV(¿2к-1) + е2к, к = ТТ^,
__(7)
т2к-1 = V(¿2к—1) + б2к-1, к = 1, N
где тк — п-мерный случайный вектор.
В соответствии с введенными шаблонами представим значения вектора V состояния изучаемой системы и значения наблюдений вектора т состояния
системы (вектора экспериментальных данных) в к-е моменты времени в виде матриц размерностью п х N:
V 4 [V (*э),...,^ (¿2Ж-1)], Н = [Н (¿2), Н (¿4),..., Н (¿2Ж )],
Нч = [Н (¿1), Н (¿з),..., Н (¿2Ж-1)], (8)
бч = [б(^2),б(^),... , )],
бнч = б(£з), . . . , б(^2Ж-1)].
С учетом введенных матриц задача линейного оценивания (7) с ошибками в инструментальных переменных примет вид
Н = ^ + бч, к = 1, N (9)
Ннч = V + бнч, к = 17^. ( )
При известных экспериментальных данных объемом 2nN, представленных матрицами Нч и Н^, необходимо оценить nN элементов матрицы V неизвестных значений вектора состояния изучаемой системы, п2 элементов матрицы А ее неизвестных параметров и п(п + 1)/2 элементов неизвестной положительно определенной ковариционной матрицы 2 размерности п х п случайных ошибок измерений {бк, к =1, 2N}, имеющих нормальный закон распределения N(0пЬ 2).
Для статистической разрешимости задачи (9) должно выполняться условие [17]
{2nN > nN + п2 + п(п + 1/2)} ^ ^ > (3п + 1)/2}. (10)
При выполнении условия (10) задача линейного оценивания с ошибками в инструментальных переменных (9) может быть представлена как задача регрессионного анализа — получение оценки неизвестной матрицы А, определяемой соотношениями
Нч = АНнч + С,
( 4 бч - Абнч, к = 1, N (11)
где каждый столбец случайной матрицы ( является независимым п-мерным случайным вектором, имеющими закон распределения N(#п1, 2^), где 2^ — ковариционная матрица размерностью п х п. Эта матрица имеет вид
2С = 2 + А2АТ (12)
и зависит от матрицы неизвестных параметров A исходной математической модели (5). С позиций классического регрессионного анализа [5] это означает, что задача оценивания (11) лишь внешне похожа на задачу классического регрессионного анализа, но таковой не является [5]. Ее разрешимость весьма проблематична. Однако можно показать [8], что при условии коммутативности матриц A и (AE)T имеют место равенства
|Е + AEAT | = |E||1n + AAT |, (Е + AEAT )-1 = (In + AAT )-1(Е)-1.
В этом случае функция правдоподобия задачи регрессионного анализа (11), (12) может быть представлена в следующем виде:
J(A, Ес|ЖчЖнч) -
- |ЕС|-N exp{- 1tr[(W4 - AW^XW - AW^fЕ--1]} =
N m N
= |Е|-N |/n + AAT|-n X X exp{- 1tr[(W4 - AWh4)(W4 - AW^f (4 + AAT)-1E-1]}, (13)
где tr[M] — след матрицы M, а W4 и Wнч — реализации случайных матриц W4 и W^ соответственно.
В работах [7, 8] показано, что реализация A оценки матрицы неизвестных параметров A определена равенством
A = W4 W+, (14)
где Wj+j — матрица, псевдообратная по отношению к матрице W^ [10]. Таким образом, реализация оценки максимального правдоподобия [5,10] для матрицы неизвестных параметров A в рассматриваемом случае совпадает с ее оценкой, полученной методом наименьших квадратов. Сама оценка максимального правдоподобия для матрицы неизвестных параметров A в рассматриваемом случае тоже совпадает с ее оценкой наименьших квадратов и
ah 4 wjwj+j.
Поскольку в дальнейшем необходимо различать вектор состояния системы и его наблюдаемые значения, то для обозначения вектора состояния воспользуемся символами Х, У и представим математическую модель (1), (3) в следующем виде:
Х((п0 + 1)т0) = а0Х(п0т0), п0 е {0, 1, 2,...};
X (0) = Х0,
^((п1 + 1)т1) = а11>(п1т1) + Ь1 Х^т1), п1 е {0, 1, 2,...}; У(0) = 0,
где Х((п0 + 1)т0) и ^((п1 + 1)т1) — численности нормальных и аномальных клеток в моменты времени £ = (п0 + 1)т0 и £ = (п1 + 1)т1 соответственно.
Напомним, что согласно (2), параметр а0, определяет характер динамики сумммарной численности популяции нормальных клеток, а параметры а1 и Ь1, согласно (4) определяют характер динамики сумммарной численности популяции аномальных клеток, при этом а1 характеризует поведение аномальных клеток, а Ь1 определяет перерождение нормальных клеток в аномальные.
Поскольку в результате эксперимента можно оценить только количество нормальных и аномальных клеток в фиксированнные моменты времени [1,2,3], соответствующие измерения (фиксацию) проводим через равные промежутки времени, равные £к. Будем предполагать, что т1 = т2 = т и положим £к = кт. Для упрощения записи будем использовать следующие обозначения: Х(£к) = = Хк, У(£к) = Ук. Через Хк и Ук обозначим наблюдаемые значения величин Хк и Ук. С учетом введенных обозначений рассматриваемая математическая модель динамики численности изучаемой популяции примет вид
а0 0 Х
* аЧ \ftj- (15)
Заметим, матрица параметров системы имеет нижнетреугольный вид, т.е. динамика развития популяции аномальных клеток не оказывает влияния на динамику развития популяции нормальных клеток. Поэтому задачу идентификации можно разбить на две независимые подзадачи: нахождение оценки параметра а0, определяющего динамику развития популяции нормальных клеток, и
нахождение оценок параметров а1 и Ь1, определяющих динамику развития популяции аномальных клеток. В общем случае решение задачи оценивания типа (11) с частично определенной матрицей А связано с преодолен ием трудностей принципиального характера [8].
Решение задачи параметрической идентификации математической модели (15) удобно провести в три этапа:
1) айти точеч ые оце ки параметров модели;
2) проверить статистическую гипотезу о виде закона распределения случайных возмущений;
3) найти законы распределения параметров модели.
Отметим, что на основе полученных законов распределения параметров модели мож о построить доверитель ые и тервалы для этих параметров, од ако этой задаче будет посвящено отдельное исследование.
В случае, если проверка статистической гипотезы о виде закона распреде-ле ия случай ых возмуще ий показала, что гипотеза о ормаль ом распределении случайных возмущений отклонена, решить задачу параметрической идентификации математической модели (15) предложенным ниже методом не удается.
Перейдем к решению задачи нахождения точечных оценок параметров модели и задачи построения плотности распределения вероятности параметров модели.
3. Основные принципы байесовского подхода
В теории статистических выводов используются два ос ов ых аправле-ния: «классический» подход, развитый в работах Дж. Неймана, Е.С. Пирсона и др., используемый при выборках большого объема, и байесовский подход, основанный на применении теоремы Байеса, часто используемый при небольших объемах данных [9] и при неизвестной ковариционной матрице Заметим, что байесовский подход приме яется в ос ов ом для оце ива ия биологических и экон омических параметров. Поскольку объем экспериментальных дан ных ограничен [1, 2, 3], то исходную задачу оценивания будем решать на основе байесовского подхода [9, 16].
Для построения функций плотности распределения вероятности параметров а0, а1, Ь1 модели (15) воспользуемся байесовским подходом и теорией инвариантности Джеффриса [9, 15].
Следуя [9, 16], приведем основные сведения о байесовском подходе для оценивания параметров. Этот подход к интерпретации истинных параметров модели ос ова а использова ии имеющейся априор ой и формации об изучаемом объекте и пересмотре этой и формации с учетом получаемых выборочных данных об исследуемом объекте.
При таком подходе полагают, что точное значение параметров неизвестно, и они являются случайными величинами с какими-то заранее заданными распределениями. Случайность рассматривается как неотъемлемое свойство объекта, присущее ему по самой природе. При этом полагают, что сам физический объект подвержен непрерывным случайным физическим изменениям. Случайность здесь можно интерпретировать как меру незнания, что позволяет проводить рассуждения в условиях неопределенности.
Идея байесовского подхода заключается в том, что объединяя априорную функцию плотности распределения вероятности вектора параметров с информацией выборки при помощи теоремы Байеса, получают апостериорную функцию плотности распределения вероятности вектора параметров.
Формула Байеса для непрерывных случайных величин может быть пред-ставле а в виде
/ (%) = (16)
Здесь /(9|у) — апостериорная функция плотности распределения вероятности вектора 9 при данной эмпирической информации относительно вектора у. Апостериорная функция плотности распределения вероятности, содержит как априор ую, так и выбороч ую и формацию и используется в байесовском подходе для построения вероятностных утверждений о параметре р(9), например для построения доверительных интервалов для параметров.
Функция /(у|9) в формуле (16)есть условная плотность распределения вероятности новых наблюдений у при определенных значениях параметра 9, рассматривается как функция от второго аргумента 9 при фиксированн ом первом аргументе, и является функцией правдоподобия для параметра р(9).
Функция ](0) в формуле (16) есть априорная функция плотности распределения вероятности для вектора 0. Она базируется на первоначальной информации о векторе параметров, полученной из предыдущих наблюдений или моделирования.
Априорное распределение р(0) вектора параметров 0 необязательно задавать точно, необходимо только определить семейство априорных распределений, к которому принадлежит р(0).
Функция f (у) в формуле (16) есть плотность распределения вероятности наблюдений у для всех возможных значений параметра 0.
Учитывая, что плотность распределения вероятности наблюдений у для всех возможных значений параметра 0 имеет вид
и является константой (f (y) = const = 0), формула Байеса может быть записана в виде
В случае, если априорная информация базируется на данных, полученных ранее в результате эксперимента, то априорную функцию плотности распределения вероятности называют априорной функцией плотности распределения вероятности, базирующейся на данных.
Если же таких данных нет или их очень мало и они носят случайный характер, то априорная информация может получена из теоретических соображений или моделированием. В подобных случаях априорную функцию плотности распределения вероятности называют априорной функцией плотности распределения вероятности, не базирующейся на данных.
При отсутствии ранее полученных данных и крайне ограниченной информации об изучаемом явлении или объекте, говорят, что априорная информация является «неясной» или «расплывчатой». Если априорная информация относится к параметрам некоторой модели и является «неясной» или «расплывчатой», то для анализа данных используют расплывчатую априорную функцию плотности распределения вероятности.
f(%) - fW(y|0).
В случае, когда значение параметра неизвестно, согласно [15] предлагаются два правила выбора априорного распределения.
1. Если параметр существует в интервале от —то до то, то его априорная вероятность должно считаться равномерно распределенной.
Данное правило для представления незнания значения предлагает принять
р(в)—в — —в, —то < в < то,
т.е. р(в) — const. Это плотность распределения вероятности является несобственной, поскольку
/то
р(в) —в = то.
то
2. Если параметр в существует в интервале от 0 до то, то априорная вероятность его логарифма должно считаться равномерно распределенной, т.е. полагают $ = log(e).
Это правило для представления незнания значения предлагает принять
p($) d$ — d$, —то <$< то,
где $ = log в. Поскольку d$ = ^, можно записать
в
—в
р(в) —в - —, 0 <в< то. в
Следуя данному подходу перейдем к рассмотрению параметров модели (15) и законов их распределения.
4. Нахождение точечных оценок параметров линейной модели
Для нахождения точечных оценок а0, а1, b1 элементов матрицы параметров системы а0, а1, b1 математической модели (15), описывающей динамику суммарных численностей субпопуляций нормальных и аномальных клеток, воспользуемся ее разностным аналогом, записанным на двухточечном шаблоне:
£)=(a:)(£)+(-
где матрицы-строки £ = [£Ь £2, ..., £м] е М1Хж(Я) и п = [П1, П2, ... , Пж] е е М1хж (Я) являются реализациями случайных возмущений из уравнений состояния в математической модели (15).
Учитывая, что задачу идентификации матрицы параметров системы можно разбить на две независимые подзадачи, сначала сформируем массивы экспериментальных данных для нахождения оценки параметра а0, затем массивы экспериментальных данных для нахождения оценок параметров а1 и Ь1.
Пусть {£2п-1 ; £2п} — п-й двухточечный шаблон уравнения состояния (1), где £2п-1 = (2п — 1)т, £2п = 2пт, в узлах которого известны экспериментальные значения Х2п—1,Х2п суммарной численностей популяций нормальных клеток в моменты фиксации £2п—1 и £2п соответственно.
Формируем массивы экспериментальных данных
Хч 4 [Х2,Х4,...,Х2ж ] е М1хж (Я),
(18)
Хнч = [X1, X3, . . . , Х2Ж —1] е М1хЖ(Я).
Согласно (1) и (18)
Хч = а0Хнч + £, (19)
где матрица-строка £ = [£1, £2,..., £Ж] е М1хЖ(Я) является реализацией случайного возмущения первого из уравнений состояния в исходной математической модели (1), обусловленного как случайными возмущениями самой математической модели, так и так и погрешностями в экспериментальных данных.
Таким образом, на основе теории псевдообратных матриц [10] в соответствии с равенством (19) можно утверждать, что реализация оценки параметра а0 методом наименьших квадратов может быть представлена в явном виде:
а° = Хч Х+, (20)
где Х+ — матрица, псевдообратная по отношению к матрице Хнч, определенной в (18).
Пусть далее {£2п—1; £2п} — п-й двухточечный шаблон уравнения состояния (3), где £2п—1 = (2п — 1)т, £2п = 2пт, в узлах которого известны экспериментальные значения У2п—1, У2п суммарной численности популяций аномальных
клеток в моменты фиксации £2п-1 и £2п соответственно. Формируем массивы экспериментальных данных:
Уч 4 [У2,У4,...,^2Ж] е М1Хж(Я), Унч 4 [У1,Уз,..., У2Ж-1] е М1ХЖ(Я),
(21)
4 У1,Уз,.. . , У2Ж-1 _ У Унч е М2хж(Я)
ХЪ Х3, . . . , Х2Ж-1 Хнч
и, согласно (3), (21), получаем
Уч = [а1, Ь1]^ + п, (22)
где матрица-строка
П = [П1,П2,... ,пж] е М1ХЖ(Я)
является реализацией случайного возмущения из уравнений состояния в математической модели (3). Таким образом, на основе теории псевдообратных матриц [10] и равенства (22) методом наименьших квадратов находим реализацию оценки параметров а1, Ь1:
[а1, Ь1] = Уч ^ +, (23)
где Z + — матрица, псевдообратная по отношению к матрице Z, определенной в (21).
5. Законы распределения параметров
Второй этап решения задачи параметрической идентификации математической модели (15) подробно рассмотрен в работе [12].
На третьем этапе решения указанной задачи необходимо найти законы распределения параметров математической модели (15) и параметры соответствующих законов распределения.
Если с заданным уровнем значимости нет оснований для отклонения основ ной гипотезы о том, что матрицы-строки невязок £,п, определенные равен ствами
£ = [£1, £2, . . . , ] = Хч — а0Хнч,
ч а Х ^ нч ,
(24)
П = [П1,П2,... ,Пм] — Гч — [а1, Ь1 ]£,
являются матрицами реализаций нормальных случайных величин, то можно доказать [6], что байесовские апостериорные законы распределения параметров а0, а1, Ь1, построенные с учетом теории инвариантности Джеффриса [9, 15], являются распределением Стьюдента.
Найдем параметры соответствующих распределений.
Построение функции плотности распределения вероятности параметра а0. Согласно (1) и (19) характер динамики сумммарной численности популяции нормальных клеток определяется выражением
Хч = а0Хнч + £ ,
где а0 — параметр, определенный в (2), матрица-строка £ = [£1, £2,... ] Е Е М1хж(Я) является реализацией случайного возмущения из уравнения состояния в исходной математической модели (1), матрицы-строки
Хч — [Х2,Х4,...,Х2Ж] Е М1хЖ(Я),
Хнч — [Х1 , . . . ,Х2Ж —1] Е М1хЖ(Я)
представляют используемую выборку.
Элементы матриц-строк Хнч и Хч будем обозначать и соответственно (к = 1, N). Эти элементы являются реализациями случайных независимых величин [6].
Примем допущение [9], что все элементы матрицы Хнч являются фиксированными нестохастическими переменными.
В случае, если гипотеза о том, что матрица-строка невязок £, определенная первым равенством (24), является матрицей реализаций нормальных случайных величин с нулевым математическим ожиданием, не была отклонена, воспользуемся основными положениями регрессионого анализа [13, 14].
Будем считать, что независимые случайные величины £к, к = 1, N, имеют распределение N(0,а2), т.е. для любого к = 1, N математическое ожида-
ние случайной величины равно нулю (М(ек) = 0), дисперсия этой величины постоянна ) = а2), случайные величины и е не коррелированы
М(е*е?) = 0 при к _ ]).
Воздействие неучтенных случайных факторов и ошибок в модели (1) опре-
2
деляется с помощью дисперсии случайных возмущений а , несмещенной оценкой которой является выборочная остаточная дисперсия
о2 _ (Хч — а°Хнч)(Хч — а°Хнч)Т _ _ 1 ^^ а2
О _ Х-1 _ (X - 1) _ X-а,
где _ ж) - а°жнч — выбочная оценка реализации случайного возмущения .
При построении функции плотности распределения вероятности параметра модели а° для оценки этого параметра а° применим метод наименьших квадратов, не прибегая к использованию псевдообратной матрицы.
Обозначим через О2 величину, пропорциональную оценке дисперсии а2 и представляющую собой сумму квадратов отклонений расчетных значений состояния а°Хнч от их экспериментальных значений Хч:
О2 _ (Хч - о°Хнч)(Хч - а°Хнч)Т _ . (25)
Оценка параметра а° методом наименьших квадратов предполагает, что вели-
N
чина е_ минимизируется. Имеем
:=1
О2 _ еет _ (Хч - а°Хнч)(Хч - а°Хнч)Т _
_ ХчХ^ — а° ХнчХч^ — ХчХта° + й°ХнчХнчй° _
_ ХчХчТ - 2а°ХнчХ + (а°)2ХнчХТ. (26)
Последнее равенство в (26) выполняется, поскольку
а°Хнч Х^ _ Хч а°,
N
так как Хнч Х^ _ ^ (ж^ ж)) и а° являются скалярными величинами.
:=1
Для нахождения экстремума продифференцируем функцию б"2 (а0) и при-
равняем нулю:
^ 2
— — 2ХнчХт + 2ХнчХТга0 — 0.
-нч^нч^
Отсюда = Х^Х^а0. Решая это уравнение относительно а0, находим
'Т а0
а = (Хнч ХНЧ) Хнч Хч = Хч ^
(27)
что полностью соответствует (20).
Поскольку Хнч Хч > 0, то найденное значение а0 есть минимум функции £2(а0).
Совместная апостериорная функция плотности распределения вероятности параметра а0 и дисперсии а2 согласно теореме Байеса [9] может быть представлена в виде
/(а0,а2|Хч,Хнч) - /(а0, а2) f (Хч|Хнч, а0, а2), (28)
где /(Хч|Хнч, а0, а2) — условная плотность распределения вероятности новых наблюдений Хч при определенных значениях параметра а0, а2 и Хнч, рассматриваемая как функция правдоподобия, а /(а0, а2) — совместная априорная функция плотности распределения вероятности величин а0 и а2. Здесь знак — обозначает пропорциональность, коэффициент пропорциональности может быть найден из условий нормировки после определения остальных параметров совместной апостериорной функции плотности распределения вероятности.
Запишем апостериорную функцию плотности распределения вероятности /(Хч|Хнч, а0, а2) при нормальном распределении ошибок {£«}, заданных значениях а0, Хнч и выборочной дисперсии а
2
N
/ (Хч|Хнч,а0,а2) — (а2)—-ехр
1
(Хч — Хнча0)(Хч — Хнча0)Т . (29)
2а2
Учитывая (25) и (27), приходим к следующим эквивалентным представлениям для произведения (Хч — Хнч а0)(Хч — Хнч а0)
0Т
(Хч — Хнч а )(Хч — Хнч а ) =
= {(Хч — Хнч а0) — (а0 — а0)Хнч}{ (Хч — Хнч а0) — (
— (а —а
ОХнч}
Т
= (Хч — Хнч а0)(Хч — Хнч а0)Т + (Х^Х^Ха0 — а0)2 =
-т
02
= £2 + (ХнчХТ,)(а0 — а0)2. (30)
02
Запишем совместную функцию правдоподобия для параметра а0 и дисперсии а2. Она задается как плотность вероятности совместного появления результа-
тов выборки жГ, хк , к = 1, N:
N
/(а0, а2^, Хч) = Ц / (жГ|а0, а2)
к=1
Рассматривая плотность /(Хч|Хнч, а 0, а2 ), определенную соотношением (29) как функцию от а 0 и а2 при фиксированных Хч и Хнч и учитывая (30), получим функцию правдоподобия
/(а0,а2|Хнч,Хч) - (а2) 2 ехр
2ч-К
2\ 2
52 - (ХчХТХа0 -а0)2 2а2
(31)
Совместную априорную функцию плотности распределения вероятности /(а 0, а2) согласно [15] будем считать расплывчатой. Также примем допущение, что параметр а0 и выборочная дисперсия а2 независимо распределены, т.е. /(а0, а2) = /(а°)/(а2).
Априорную функцию плотности распределения вероятности параметра а0 согласно [15] будем считать константой, поскольку о значении а 0 ничего сказать нельзя и формально а0 может принимать значения от —то до то. При этом
/ (а 0) = сопб1;.
(32)
Заметим, что данное предположение не ограничивает общности, так как совместная апостериорная функция плотности распределения вероятности будет найдена с точностью до константы. Поскольку а2 может принимать значения от 0 до то, априорная функция плотности распределения вероятности для а, будет определяться известным соотношением [15]:
/(а2) - 1.
а
С учетом (31), (32) и (33) выражение (28) примет вид
(33)
/(а0, а2|Хч,Хнч) - (а2)-^ ехр
2\-
(52 + (ХнЧХТ)(а0 — а0)2) 2а 2
(34)
Для дальнейших исследований положим х = а 2, тогда выражение (34) можно переписать следующим образом:
N +1
/(а0,х|Хч,Хнч) — (х) 2 ехр
(£2 + (Хнч Х1Т1)(а0 — а0)2 )х 2
(35)
Для получения апостериорной функции плотности распределения вероятности параметра а0 необходимо проинтегрировать (35) по х . В результате получим
/(а0|Хч,Хнч) — (£2 + (ХнчХнч)(а0 — а0)2)—* х
N-1
^ X 2 ехр
х
—1(52 + (ХнчХТ,)(а0 — а0)2)
(£2 + (ХнчХнч)(а0 — а0)2)—я/2
¿х- (36)
Из (36) следует, что плотность распределения вероятности параметра а0 может быть представлена в следующем виде:
/ (а0|Хч,Хнч) — {^2 + (ХнчХнчТ )(а0 — а0)2}—*,
(37)
где £2 4 (Хч — а0Хнч)(Хч — а0Хнч)Т = иТ
Построение функций плотности распределения вероятности параметров модели а1 и Ь1. Для построения совместной функции плотности распределения вероятности параметров модели а1 и Ь1 воспользуемся предложенным выше способом.
Напомним, что согласно (22) вектор экспериментальных данных (численность аномальных клеток) может быть представлен в виде
Уч = [а1, Ь1] £ + п,
где матрица-строка п = [пъ П2, - - -, Пж] е (Л) является реализацией случайного возмущения из уравнения состояния в математической модели (3), а матрица-строка Уч и матрица £ согласно (21) имеют вид
Уч 4 [У2,У4,...,У2Ж ] е М1хж (Л),
£ 4
У У нч
Хнч
е М2хм(Л),
где
Унч 4 [Yi, Ys,..., Y2N-1] g Mixn(R),
Хнч — [X1, Xs, . . . , X2N-1] G M1xN(R).
Далее элементы матриц Z и Y4 будем обозначать zk и уЧ соответственно
(k = TTM).
Зависимая переменная Y4 есть линейная функция от Z с точностью до реализации п случайного возмущения.
Как и при построении функции плотности распределения вероятности параметра модели а0, для построения совместной функции плотности распределения вероятности параметров модели а1 и b1 примем допущение [9], что все элементы матрицы Z являются фиксированными нестохастическими переменными.
Если гипотеза о том, что матрица-строка невязок п, определенная вторым равенством (24), является матрицей реализаций нормальных случайных величин с нулевым математическим ожиданием и неизвестной дисперсией а , не была отклонена, то считаем, что случайные величины Пк, k = T, N (или зависимые переменные уч с учетом допущения об их нестохастической природе), являются независимыми нормально распределенными случайными величины с математическим ожиданием, равным нулю (M(пк) = 0, k = T, N) и неизвестной дисперсией а2 [9].
Применим метод наименьших квадратов для оценки вектора параметров модели [а1,Ь1]. Обозначим через S2 величину, равную сумму квадратов отклонений расчетных значений состояния [a1,b1] Z от их экспериментальных значений Уч, пропорциональную оценке дисперсии а2:
N
S2 — (Уч - [а1, b1]Z)(Уч - [а1, b1]Z)T = ппТ = £ П2 (38)
k=1
где пк — выбочная оценка случайной величины пк, П = Уч - [а1, b1]Z.
Оценка вектора параметров [а1, b1 ] методом наименьших квадратов предполагает, что величина S2 минимизируется:
S2 = ппТ = (Уч - [а1, b1]Z)(Уч - [а1, b1]Z)T ^ min. (39)
Раскрыв скобки в (39) и учитывая, что
Уч([а 1,Ь1]^)Т = У^Т [а 1,Ь1]Т = [а 1,61]^Учт,
Тг ¿1 а 1Т г-1 а 1 ^ллТ
поскольку матрицы [а1, Ь1]^УчТ и ^ [а1, Ь1^ имеют размерности 1 х 1, полу им
Б2 = УчУчТ — 2[ а1, Ь^УчТ + [а 1,Ь1]^^Т[а 1,Ь1]Т. (40)
Далее, дифференцируя функцию 52(а1, Ь1), определенную соотношением (40), и используя необходимое условие экстремума функции нескольких переменных полу им
—2^УчТ )Т + 2([а 1,Ь1])(^^Т) = 0.
Из последнего соотношения находим оценку вектора [а1, Ь1 ] в виде:
уТ г-1 ги Т
[а 1,ь1] = (Уч^Т )—1 =
Т\( пггуТл-1
(41)
Поскольку матрица ZZТ является положительно определенной, найденное значение [а1, Ь1] соответствует минимуму функции 52. Отметим, что полученное соотношение совпадает с (23).
Совместная апостериорная функция плотности распределения вероятности вектора параметров [а1, Ь1 ] и дисперсии а2 2 согласно формуле Байеса [9] может быть представлена в виде
/(а 1,Ь1,а2|Уч , Z) - /(а1, Ь1, а2)/(У,|Z, а\Ь\а2),
1 1 2
1 1 2
(42)
где /(У!| Z, а1, Ь1 , а2) — условная плотность распределения вероятности новых наблюдений У! при определенных значениях вектора параметров [а1, Ь1], £ и Z, рассматриваемая как функция правдоподобия,
/ (а1,Ь1,а2) — совместная априорная функция плотности распределения вероятности для вектора параметров [ а1, Ь1] и а2.
Запишем апостериорную функцию плотности распределения вероятности /(У!а1, Ь1, а2) при нормальном распределении ошибок п», заданных знач е-ниях Z, а1, Ь1 и дисперсии а2:
/(Уч|Z,а1,Ь1,а2) - (а2Г4/2
ехр
1
(Уч — [а )(Уч — [а )
11
2а2
Т
. (43)
Выражение
(Уч - [а1, 6^2)(УЧ - [а1, 6^2)т с учетом (41) может быть представлено в следующем виде:
(Уч - [а1, 6^2)(Уч - [а1, Ь^2)т =
= {(Уч - [а1, ) - [а1 - а1, б1 - 6^2} х
т
х {(Уч - [а1, ) - [а1 - а1, б1 - Ь1^} = = (Уч - [а1, )(Уч - [а1 ,6^2)т -
- (Уч - [а1, ь1^)([а1 - а1, ь1 - Ь1^)т -
- ([а1 - а1, ь1 - Ь1]^)(Уч - [а1, Ь1]^)т +
+ [а1 - а1, Ь1 - Ь1 ]22т [а1 - а1, Ь1 - Ь1]7 =
I1 - а1,Ь1 - Ь1 ]22т [а1 = 52 - (Уч - Уч2т(22т)-12)([а1 - а1, Ь1 - 6^2)т -- ([а1 - а1, Ь1 - Ь1]^)(Уч - Уч2т(22т)-12)т +
+ [а1 - а1, Ь1 - ^22т[а1 - а1, Ь1 - Ь1 ]т.
Учитывая, что матрица 2 т (22т )-12 равна единичнои матрице, получаем
52 + [а1 - а1, Ь1 - Ь1]22т[а1 - а1, Ь1 - Ь1]т. (44)
Функция правдоподобия для параметров а1,Ь1 и а2 задается как плотность распределения вероятности совместного появления результатов выборки ,
Уч, к = 1,Х:
N
/(а1, Ь1, а212, Уч) = П /(**, УЧ 1а1, Ь1, а2).
ч=1
Получить функцию правдоподобия для параметров а1, Ь1 и а2 можно, рассматривая (43) как функцию от а1, Ь1 и а2 при фиксированных Уч и 2. С учетом (44) получим:
/(а1, Ь1, а2|2, Уч) - (а2)-^ ехр
1 2 -
2а2
(а1 - а1, Ь1 - Ь1 )22т (а1 - а1, Ь1 - Ь1)т
. (45)
Априорную функцию плотности распределения вероятности параметров а1, b1
Совместную априорную функцию плотности распределения вероятности f (а1, b1, а2) согласно [15, 9] будем считать расплывчатой. Также примем допущение, что элементы вектора [ а1, b1] и дисперсия а2 независимо распределены, т.е.
f (а1, b1, а2) = f (а 1,b1) f (а2).
1,
согласно [15] положим константой, поскольку о значении параметров а1 и b1 ничего сказать нельзя, и формально а1 и b1 могут принимать значения от —то до +то. При этом
f (а 1,b1) = const. (46)
Поскольку дисперсия а2 может принимать значения из интервала (0, +то), ее априорная функция плотности распределения вероятности, согласно [15], будет иметь вид
f (а2) - 1. (47)
а
С учетом (45), (46) и (47) соотношение (42), определяющее совместную апосте-риорною функцию плотности распределения вероятности неизвестных параметров а1, b1 и выборочной дисперсии а2, может быть представлено следующим образом:
f (а1, b1, а2|Уч, Z) - (а2)-^ х
1
х exp
(S2 + [а1 — a1,b1 — b1]ZZT[а1 — а 1,b1 — b1]7) . (48)
2а2
Для дальнейших исследований положим
X = а—2. (49)
Тогда выражение (48) примет вид
N +1
f(а 1,b1,x|Y4,Z) - х* х
1
х exp
(S2 + [а1 — а1, b1 — b1]ZZT[а1 — a1,b1 — b1]7)х
. (50)
Для полу ения апостериорной функции плотности распределения вероятности параметров модели а1 и Ь1 необходимо проинтегрировать (50) по х.
В результате получим
_ N
/(а1,Ь1|К1,^) - (52 + [ЯаЛ]ЯЯт[¿аЛ]Т)-2 х
ехр
х / -
■ N
(52 + [£аА]^т [¿«А]т) -2
¿X, (51)
где ¿а = а1 — а 1, = Ь1 — Ь1.
Интеграл в (51) сходится и его значение равно константе. Поэтому из (51) следует, что плотность совместного распределения вероятности параметров модели а1 и Ь1 может быть представлена как
Л Л Л _ N
f (а1, Ь1|К1, ^) - {52 + [а1 — а1, б1 — Ь1]^[а1 — а 1, Ь1 — Ь1 ]т} 2 , (52)
где
52 4 (Уч — [а1,ь1]^)(Уч — [а1,Ь1]^)т = .
6. Заключение
Разработана методика получения точечных оценок параметров линейной модели, описывающей динамику селективного размножения клонообразую-щей популяции аномальных клеток в культуре стволовых клеток ч еловека при стандартных лабораторных условиях культивирования без у ета ограни ений питания.
На основе баесовского подхода получены функции плотности распределения вероятности параметров а0, а1 и Ь1 математич еской модели.
На основе найденных функций плотности распределения вероятностей можно найти интервальные оценки для параметров модели и использовать их для анализа возможных сценариев развития популяции.
Работа выполнена при поддержке грантов РФФИ № 12-07-00329 и № 1207-00267.
Список литературы
1. Бочков Н.П., Никитина В.А. Цитогенетика стволовых клеток человека // Молекулярная медицина. 2008. № 3. С. 40-47.
2. Бочков Н.П., Никитина В.А., Рослова Т.А., Чаушев И.Н., Якушина И.И. Клеточная терапия наследственных болезней // Вестник РАМН. 2008. № 10. С. 20-28.
3. Бочков Н.П., Никитина В.А., Воронина Е.С., Кулешов Н.П. Методическое пособие по тестированию клеточных трансплантатов на генетическую безопасность // Клеточные технологии в биологии и медицине. 2009. № 4. C.183-189.
4. Бочков Н.П., Веденков В.Г., Волков И.К. и др. Построение математической модели для оценки соотношения клеток, прошедших разное число делений в культуре // Доклады АН СССР. 1984. Т. 274, № 1. C. 186-189.
5. Демиденко Е.З. Линейная и нелинейная регрессии. М.: Финансы и статистика, 1981. 302 с.
6. Волков И.К. Условия идентифицируемости математических моделей эволюционных процессов по дискретным косвенным измерениям вектора состояния // Известия РАН. Теория и системы управления. 1994. № 6. C. 55-72.
7. Волков И.К. Идентифицируемость математических моделей эволюционных процессов // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2005. № 3. C. 64-73.
8. Волков И.К. Математическое моделирование эволюционных процессов в системе генетического мониторинга по экспериментальным данным: дисс. ... докт. физ.-мат. наук. М., 1992. 281 с.
9. Зельнер А. Байесовские методы в эконометрии. - М.: Статистика, 1980. -440 с. (Zellner A. An introduction to Bayesian inference in econometrics. John Wiley and Sons, New York, 1971, 448 p.)
10. Алберт А. Регрессия, псевдоинверсия и рекурентное оценивание: пер. с англ. М.: Наука, 1977. 224 с.
11. Бочков Н.П., Виноградова М.С., Волков И.К., Кулешов И.К. Математическая модель суммарных численостей взаимодействующих клеточных популяций // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2011. № 1. С. 18-24.
12. Бочков Н.П., Виноградова М.С., Волков И.К. Оценка вероятности реализации вариантов развития взаимодействующих клеточных популяций // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2011. № 3. С. 31-43.
13. Кремер Н. Ш. Теория вероятностей и математическая статистика. М.: Юнити, 2006. 574 с.
14. Шметтерер Л. Введение в математическую статистику. М.: Наука, 1976. 520 с. (Schmetterer L. Einführung in die mathematische Statistik. SpringerVerlag, Wien-New York, 1966, 567 p.)
15. Jeffreys H. Theory of probability. Oxford: Clarendon Press, 1983. 459 с.
16. Боровков А. А. Математическая статистика. М.: Наука, 1984. 472 с.
17. Закс Ш. Теория статистических выводов: пер. с англ. М.: Мир, 1975. 776 с.
18. Прудников А.П., Брычков Ю.А., Маричев О.И. Интегралы и ряды. М.: Наука, 1981. 800 с.
SCIENTIFIC PERIODICAL OF THE BAUMAN MSTU
SCIENCE and EDUCATION
EL № FS77 - 48211. №0421200025. ISSN 1994-0408
electronic scientific and technical journal
Parametrical identification of cooperating cellular populations model on the basis of Bayesian inference # 11, November 2012 DOI: 10.7463/1112.0490900 Vinogradova M. S.
Russia, Bauman Moscow State Technical University
mathmod@bmstu.ru
The algorithm of receiving point estimations of model parameters is developed for the mathematical model describing dynamics of selective reproduction of klone-forming population of abnormal cells in culture of human stem cells under standard laboratory conditions of cultivation without restrictions of a delivery on limited sampling of experimental data about condition vector values. Overviews on Bayesian inference and Jeffreys's invariancy theory are given. Probability density functions of distribution of the specified mathematical model parameters are received on the basis of these theories.
References
1. Bochkov N.P., Nikitina V.A. Tsitogenetika stvolovykh kletok cheloveka [Cytogenetics of human stem cells]. Molekuliarnaia meditsina [Molecular medicine], 2008, no. 3, pp. 40-47.
2. Bochkov N.P., Nikitina V.A., Roslova T.A., Chaushev I.N., Iakushina I.I. Kle-tochnaia terapiia nasledstvennykh boleznei [Cell therapy of inherited diseases]. Vestnik RAMN [Annals of the Russian Academy of Medical Sciences], 2008, no. 10, pp. 20-28.
3. Bochkov N.P., Nikitina V.A., Voronina E.S., Kuleshov N.P. Metodicheskoe posobie po testirovaniiu kletochnykh transplantatov na geneticheskuiu bezopas-nost' [Methodical manual for testing cell transplants for genetic safety]. Kle-
tochnye tekhnologii v biologii i meditsine [Cell technologies in biology and medicine], 2009, no. 4, pp. 183-189.
4. Bochkov N.P., Vedenkov V.G., Volkov I.K. i dr. Postroenie matematicheskoi modeli dlia otsenki sootnosheniia kletok, proshedshikh raznoe chislo delenii v kul'ture [Constructing a mathematical model for the evaluation of the ratio of the cells that have passed different number of divisions in the culture]. Doklady ANSSSR [Reports of Academy of Sciences of the USSR], 1984, vol. 274, no. 1, pp. 186-189.
5. Demidenko E.Z. Lineinaia i nelineinaia regressii [Linear and nonlinear regression]. Moscow, Finansy i statistika, 1981. 302 p.
6. Volkov I.K. Usloviia identifitsiruemosti matematicheskikh modelei evoliut-sionnykh protsessov po diskretnym kosvennym izmereniiam vektora sos-toianiia [Identifiability conditions of mathematical models of evolutionary processes over discrete indirect measurements of the state vector]. Izvestiia Rossiiskoi akademii nauk. Teoriia i sistemy upravleniia [Bulletin of the Russian Academy of Sciences. Theory and systems of control], 1994, no. 6, pp. 55-72.
7. Volkov I.K. Identifitsiruemost' matematicheskikh modelei evoliutsionnykh protsessov [Feasibility of identification of mathematical models of evolutionary processes]. VestnikMGTUim. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2005. no. 3. pp. 64-73.
8. Volkov I.K. Matematicheskoe modelirovanie evoliutsionnykh protsessov v sisteme geneticheskogo monitoringa po eksperimental'nym dannym. Dokt. diss. [Mathematical modeling of evolutionary processes in the system of genetic monitoring from experimental data. Dr. diss.]. Moscow, 1992. 281 p.
9. Zellner A. An introduction to Bayesian inference in econometrics. John Wiley & Sons, New York, 1971,448 p. (Zel'ner A. Baiesovskie metody v ekonometrii. Moscow, Statistika, 1980. 440 p.)
10. Albert A. Regression, And The Moor-Penrose Pseudoinverse. New York and London, Academic Press, 1972. 180 p. (Mathematics in Science and Engineering, vol. 94). (Russ. ed.: Albert A. Regressiia, psevdoinversiia i rekurentnoe otsenivanie. Moscow, Nauka, 1977. 224 p.)
11. Bochkov N.P., Vinogradova M.S., Volkov I.K., Kuleshov I.K. Matematich-
eskaia model' summarnykh chislenostei vzaimodeistvuiushchikh kletochnykh populiatsii [Mathematical model of dynamics of total quantities of interacting cell's populations]. VestnikMGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2011, no. 1, pp. 18-24.
12. Bochkov N.P., Vinogradova M.S., Volkov I.K. Otsenka veroiatnosti realizatsii variantov razvitiia vzaimodeistvuiushchikh kletochnykh populiatsii [Estimating the probability of implementation of variants of development of interacting cell's populations]. VestnikMGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2011, no. 3, pp. 31-43.
13. Kremer N.Sh. Teoriia veroiatnostei i matematicheskaia statistika [Theory of probability and mathematical statistics]. Moscow, Iuniti, 2006. 574 p.
14. Schmetterer L. Einführung in die mathematische Statistik. Springer-Verlag, Wien-New York, 1966, 567 p. (in German). (Russ. ed.: Shmetterer L. Vve-denie v matematicheskuiu statistiku. [Introduction to mathematical statistics]. Moscow, Nauka, 1976. 520 p.)
15. Jeffreys H. Theory of probability. Oxford Clarendon Press. 1983. 459 p.
16. Borovkov A.A. Matematicheskaia statistika [Mathematical statistics]. Moscow, Nauka, 1984.472 p.
17. Zacks S. The Theory of Statistical Inference John Wiley & Sons, 1971. 626 p. (Russ. ed.: Zaks Sh. Teoriia statisticheskikh vyvodov. Moscow, Mir, 1975. 776 p.)
18. Prudnikov A.P., Brychkov Iu.A., Marichev O.I. NmjIntegraly i riady [Integrals and series]. Moscow, Nauka, 1981. 800 p.
