Подгонка линейной модели для части данных и экспоненциальной модели роста для остальных с помощью pymc3 - PullRequest
Новая
       61

Подгонка линейной модели для части данных и экспоненциальной модели роста для остальных с помощью pymc3

1 голос
/ 14 марта 2020

Давайте сгенерируем некоторые тестовые данные: 

import numpy as np

observed = np.hstack([np.arange(10)*0.1 + 3, np.exp(np.arange(15)*.2 + .2)])

Это выглядит так: 

import matplotlib.pyplot as plt

plt.plot(observed, marker='o', linestyle="None");

enter image description here

Я бы хотел для подгонки линейной модели (y = a + bx) к первой части данных и модели экспоненциального роста (y = exp(a + bx)) ко второй части данных. НО давайте притворимся, что я априори не знаю, где находится точка переключения.

Я пытался написать эту модель в pymc3: 

x = np.arange(len(observed))

with pm.Model() as model:
    sigma_1 = pm.HalfCauchy("sigma_1", beta=10)
    alpha_1 = pm.Normal("α_1", 0, sigma=20)
    beta_1 = pm.Normal("β_1", 0, sigma=20)
    sigma_0 = pm.HalfCauchy("sigma_0", beta=10)
    alpha_0 = pm.Normal("α_0", 0, sigma=20)
    beta_0 = pm.Normal("β_0", 0, sigma=20)

    switchpoint = pm.DiscreteUniform("switchpoint", lower=0, upper=len(x) - 1)

    exponential_growth = pm.Normal(
        "exponential_growth",
        mu=np.exp(alpha_1 + beta_1 * x),
        sigma=sigma_1,
        observed=observed,
    )

    linear_growth = pm.Normal(
        "linear_growth", mu=alpha_0 + beta_0 * x, sigma=sigma_0, observed=observed
    )

    likelihood = pm.math.switch(switchpoint >= x, linear_growth, exponential_growth)

    trace = sample(2000, cores=2)

хотя, конечно, это просто подходит две модели ко всем данным. Я неправильно комбинирую их.

Как правильно указать, что я хотел бы использовать linear_model до switchpoint и exponential_model после него?

1 Ответ

0 голосов
/ 14 марта 2020

Переключатель может использоваться для вычисления соответствующих mu и sigma для модели наблюдения с чем-то вроде: 

switchpoint = pm.DiscreteUniform("switchpoint", lower=0, upper=len(x) - 1)

mu, sigma = pm.math.switch(switchpoint >= x, 
                           (alpha_0 + beta_0 * x, sigma_0),
                           (np.exp(alpha_1 + beta_1 *x), sigma_1))


lik = pm.Normal('obs', mu=mu, sigma=sigma, observed=observed)
...