построение модели PyMC3, включающей различные измерения - PullRequest
Новая
       5

построение модели PyMC3, включающей различные измерения

2 голосов
/ 17 мая 2019

Я пытаюсь включить различные типы и копии измерений в одну модель в PyMC3.

Рассмотрим следующую модель: P (t) = P0 * exp (-k B t) где P (t), P0 и B - концентрации.к это скорость.Мы измеряем P (t) в разное время и B один раз, все путем подсчета частиц.k является интересующим параметром, который мы пытаемся вывести.

Мой вопрос состоит из двух частей: (1) Как объединить измерения P (t) и B в одну модель?(2) Как использовать переменное число повторяющихся экспериментов для информирования о значении k?

Я думаю, что могу ответить на часть (1), но не уверен в том, правильно ли это или сделано в правильном виде,Мне не удалось обобщить код, чтобы включить переменное число повторений.

Для одного эксперимента (один повтор): 

ts=np.asarray([time0,time1,...])
counts=np.asarray([countforB,countforPattime0,countforPattime1,...])
basic_model = pm.Model()
with basic_model:
    k=pm.Uniform('k',0,20)
    B=pm.Uniform('B',0,1000)
    P=pm.Uniform('P',0,1000)
    exprate=pm.Deterministic('exprate',k*B)
    modelmu=pm.math.concatenate(B*(np.asarray([1.0]),P*pm.math.exp(-exprate*ts)))
    Y_obs=pm.Poisson('Y_obs',mu=modelmu,observed=counts))

Я попытался включить разные реплики в соответствии с приведенным выше описанием,но безрезультатно: 

...
    k=pm.Uniform('k',0,20) # same within replicates
    B=pm.Uniform('B',0,1000,shape=numrepl) # can vary between expts.
    P=pm.Uniform('P',0,1000,shape=numrepl) # can vary between expts.
    exprate=???
    modelmu=???

1 Ответ

2 голосов
/ 17 мая 2019

Несколько наблюдаемых

PyMC3 поддерживает несколько наблюдаемых, то есть вы можете добавить несколько объектов RandomVariable к графику с установленным аргументом observed.

Single Trial

В вашем первом случае это дало бы некоторую ясность модели: 

counts=[countforPattime0, countforPattime1, ...]

with pm.Model() as single_trial:
    # priors
    k = pm.Uniform('k', 0, 20)
    B = pm.Uniform('B', 0, 1000)
    P = pm.Uniform('P', 0, 1000)

    # transformed RVs
    rate = pm.Deterministic('exprate', k*B)
    mu = P*pm.math.exp(-rate*ts)

    # observations
    B_obs = pm.Poisson('B_obs', mu=B, observed=countforB)
    Y_obs = pm.Poisson('Y_obs', mu=mu, observed=counts)

Multiple Trials

С этой дополнительной гибкостью,надеюсь, это делает переход к нескольким испытаниям более очевидным.Он должен выглядеть примерно так: 

B_cts = np.array(...) # shape (N, 1)
Y_cts = np.array(...) # shape (N, M)
ts = np.array(...)    # shape (1, M)

with pm.Model() as multi_trial:
    # priors
    k = pm.Uniform('k', 0, 20)
    B = pm.Uniform('B', 0, 1000, shape=B_cts.shape)
    P = pm.Uniform('P', 0, 1000, shape=B_cts.shape)

    # transformed RVs
    rate = pm.Deterministic('exprate', k*B)
    mu = P*pm.math.exp(-rate*ts)

    # observations
    B_obs = pm.Poisson('B_obs', mu=B, observed=B_cts)
    Y_obs = pm.Poisson('Y_obs', mu=mu, observed=Y_cts)

Может быть несколько дополнительных синтаксических элементов для правильного умножения матриц, но это по крайней мере включает в себя правильные формы.

Приоры

Как только вы запустите эту настройку, вам будет интересно пересмотреть приоры.Я подозреваю, что у вас есть больше информации о типичных значениях для тех, которые включены в настоящее время, тем более что это похоже на химическую или физическую модель.

Например, сейчас модель говорит:

Мы полагаем, что истинное значение B остается неизменным на протяжении всего испытания, но между испытаниями является совершенно произвольным значением в диапазоне от 0 до 1000 , и его повторное измерение в ходе испытания будет распределено по Пуассону..

Как правило, следует избегать усечений, если они не исключают бессмысленные значения.Следовательно, нижняя граница 0 хороша, но верхние оценки произвольны.Я бы порекомендовал взглянуть на Stan Wiki по выбору приоров .

...