Обновление модели на PyMC3 новыми данными наблюдений - PullRequest
Новая
       14

Обновление модели на PyMC3 новыми данными наблюдений

0 голосов
/ 08 ноября 2018

В прошлом году я измерил диаметр 80 плодов, и после проверки наилучшего распределения значений я создал модель PyMC3 

with Model() as diam_model:
    mu = Normal('mu',mu=57,sd=5.42)
    sigma = Uniform('sigma',0,10)

после того, насколько я понимаю, я "натренировал" модель с моими предыдущими данными (80 значений) 

with diam_model:
    dist = Normal('dist',mu=mu,sd=sigma, observed=prior_data.values)

with diam_model:
    samples=fit().sample(1000)

затем я использовал plot_posterior из samples, возвращая также среднее значение и HPD.

Моя идея состоит в том, чтобы снова провести измерения в этом году с использованием байесовского обновления, чтобы уменьшить размер выборки. Как я могу добавить отдельные значения и обновить апостериор, ожидая, что HPD становится все меньше и меньше?

1 Ответ

0 голосов
/ 19 ноября 2018

Предполагаемая плотность ядра обновлена. Приоры

Используя другой ответ, предложенный как дубликат, можно извлечь приблизительные версии приоров, используя код из этого ноутбука Jupyter .

Первый раунд

Я предполагаю, что у нас есть данные из первого раунда выборки, которые мы можем наложить среднее значение 57,0 и стандартное отклонение 5,42. 

import numpy as np
import pymc3 as pm
from sklearn.preprocessing import scale
from scipy import stats

# generate data forced to match distribution indicated
Y0 = 57.0 + scale(np.random.normal(size=80))*5.42

with pm.Model() as m0:
    # let's place an informed, but broad prior on the mean
    mu = pm.Normal('mu', mu=50, sd=10)
    sigma = pm.Uniform('sigma', 0, 10)

    y = pm.Normal('y', mu=mu, sd=sigma, observed=Y0)

    trace0 = pm.sample(5000, tune=5000)

Извлечение новых приоров из задней части

Затем мы можем использовать результаты этой модели для извлечения постеров KDE по параметрам со следующим кодом из записной книжки : 

def from_posterior(param, samples, k=100):
    smin, smax = np.min(samples), np.max(samples)
    width = smax - smin
    x = np.linspace(smin, smax, k)
    y = stats.gaussian_kde(samples)(x)

    # what was never sampled should have a small probability but not 0,
    # so we'll extend the domain and use linear approximation of density on it
    x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
    y = np.concatenate([[0], y, [0]])
    return pm.Interpolated(param, x, y)

Второй раунд

Теперь, если у нас будет больше данных, мы сможем запустить новую модель, используя обновленные априоры KDE: 

Y1 = np.random.normal(loc=57, scale=5.42, size=100)

with pm.Model() as m1:
    mu = from_posterior('mu', trace0['mu'])
    sigma = from_posterior('sigma', trace0['sigma'])

    y = pm.Normal('y', mu=mu, sd=sigma, observed=Y1)

    trace1 = pm.sample(5000, tune=5000)

И аналогичным образом можно использовать эту трассировку для извлечения обновленных апостериорных оценок для будущих циклов поступающих данных.

Конъюгат Модель

Приведенная выше методология дает приближения к истинно обновленным априорам и будет наиболее полезна в тех случаях, когда сопряженные априоры невозможны Следует также отметить, что я не уверен, в какой степени такие приближения на основе KDE вносят ошибки и как они распространяются в модели при многократном использовании. Это хитрый трюк, но следует осторожно использовать его в производстве без дальнейшей проверки его надежности.

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

Нормально-обратная гамма-модель

Нормально-обратная гамма-модель оценивает как среднее значение, так и дисперсию наблюдаемой нормальной случайной величины. Среднее смоделировано с нормальным априорным значением; дисперсия с обратной гаммой. Эта модель использует четыре параметра для: 

mu_0  = prior mean
nu    = number of observations used to estimate the mean
alpha = half the number of obs used to estimate variance
beta  = half the sum of squared deviations

Учитывая вашу исходную модель, мы могли бы использовать значения 

mu_0  = 57.0
nu    = 80
alpha = 40
beta  = alpha*5.42**2

Затем вы можете построить логарифмическую вероятность предыдущего, как показано ниже: 

# points to compute likelihood at
mu_grid, sd_grid = np.meshgrid(np.linspace(47, 67, 101), 
                               np.linspace(4, 8, 101))

# normal ~ N(X | mu_0, sigma/sqrt(nu))
logN = stats.norm.logpdf(x=mu_grid, loc=mu_0, scale=sd_grid/np.sqrt(nu))

# inv-gamma ~ IG(sigma^2 | alpha, beta)
logIG = stats.invgamma.logpdf(x=sd_grid**2, a=alpha, scale=beta)

# full log-likelihood
logNIG = logN + logIG

# actually, we'll plot the -log(-log(likelihood)) to get nicer contour
plt.figure(figsize=(8,8))
plt.contourf(mu_grid, sd_grid, -np.log(-logNIG))
plt.xlabel("$\mu$")
plt.ylabel("$\sigma$")
plt.show()

enter image description here

Обновление параметров

Учитывая новые данные, Y1, каждый обновляет параметры следующим образом: 

# precompute some helpful values
n = Y1.shape[0]
mu_y = Y1.mean()

# updated NIG parameters
mu_n = (nu*mu_0 + n*mu_y)/(nu + n)
nu_n = nu + n
alpha_n = alpha + n/2
beta_n = beta + 0.5*np.square(Y1 - mu_y).sum() + 0.5*(n*nu/nu_n)*(mu_y - mu_0)**2

Для иллюстрации изменений в модели давайте сгенерируем некоторые данные из немного другого распределения и затем построим результирующую заднюю логарифмическую вероятность: 

np.random.seed(53211277)
Y1 = np.random.normal(loc=62, scale=7.0, size=20)

, что дает

enter image description here

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

...