Надежная нелинейная регрессия с использованием PyMC(2)

Этот вопрос аналогичен Подходит для нелинейного к данным/наблюдениям с помощью pyMCMC/pyMC, поскольку я пытаюсь выполнить нелинейную регрессию с помощью PyMC.

Однако мне было интересно, знает ли кто-нибудь, как заставить мою наблюдаемую переменную следовать ненормальному распределению (т. Е. T-распределению) с использованием PyMC. Я знаю, что они включают T-распределения, но я не был уверен, как включить их в качестве наблюдаемой переменной.

Вот краткая демонстрация с использованием некоторых поддельных данных, где я сталкиваюсь с проблемами: я хотел бы использовать выходное распределение, которое защищает от некоторых из явно выпадающих точек данных.

import numpy as np
import pymc as pm
import matplotlib.pyplot as plt

# For reproducibility
np.random.seed(1234)

x = np.linspace(0, 10*np.pi, num=150)

# Set real parameters for the sinusoid
true_freq = 0.9
true_logamp = 1.2
true_decay = 0.12
true_phase = np.pi/4


# Simulate the true trajectory
y_real = (np.exp(true_logamp - true_decay*x) *
          np.cos(true_freq*x + true_phase))

# Add some noise
y_err = y_real + 0.05*y_real.max()*np.random.randn(len(x))


# Add some outliers
num_outliers = 10
outlier_locs = np.random.randint(0, len(x), num_outliers)
y_err[outlier_locs] += (10 * y_real.max() *
                        (np.random.rand(num_outliers)))



# Bayesian Regression

def model(x, y, p0):

    log_amp = pm.Normal('log_amp', mu=np.log(p0['amplitude']),
                        tau=1/(np.log(p0['amplitude'])))

    decay = pm.Normal('decay', mu=p0['decay'],
                        tau=1/(p0['decay']))

    period = pm.TruncatedNormal('period', mu=p0['period'],
                                tau=1/(p0['period']),
                                a=1/(0.5/(np.median(np.diff(x)))),
                                b=x.max() - x.min())

    phase = pm.VonMises('phase', mu=p0['phase'], kappa=1.)

    obs_tau = pm.Gamma('obs_tau', 0.1, 0.1)

    @pm.deterministic(plot=False)
    def decaying_sinusoid(x=x, log_amp=log_amp, decay=decay,
                          period=period, phase=phase):

        return (np.exp(log_amp - decay*x) *
                np.cos((2*np.pi/period)*x + phase))

    obs = pm.Normal('obs', mu=decaying_sinusoid, tau=obs_tau, value=y,
                    observed=True)

    return locals()

p0 = {
    'amplitude' : 2.30185,
    'decay'     : 0.06697,
    'period'    : 7.11672,
    'phase'     : 0.93055,
}


MDL = pm.MCMC(model(x, y_err, p0))
MDL.sample(20000, 10000, 1)

# Plot fit
y_min = MDL.stats()['decaying_sinusoid']['quantiles'][2.5]
y_max = MDL.stats()['decaying_sinusoid']['quantiles'][97.5]
y_fit = MDL.stats()['decaying_sinusoid']['mean']
plt.plot(x, y_err, '.', label='Data')
plt.plot(x, y_fit, label='Fit')
plt.plot(x, y_real, label='True')

plt.fill_between(x, y_min, y_max, color='0.5', alpha=0.5)
plt.legend()

ненадежная посадка :(

Спасибо!!


python regression robust pymc
person pstjohn    schedule 24.10.2014    source источник

Ответы (1)


arrow_upward
2
arrow_downward

PyMC2 имеет встроенное распределение t, pm.T, но оно центрировано на нуле, поэтому вы не можете использовать его непосредственно в этом приложении. Вместо этого вы можете использовать функцию pm.t_like(x, nu), которая вычисляет логарифмическое правдоподобие по значению x и параметру глубины резкости nu, чтобы определить пользовательский наблюдаемый стохастик. Сделать такое пользовательское распределение для наблюдаемой переменной очень просто: измените строки 59-60 на:

@pm.observed
def obs(mu=decaying_sinusoid, tau=obs_tau, value=y):
    return pm.t_like(value-mu, tau)
person Abraham D Flaxman    schedule 27.10.2014
comment
Спасибо! Это работает, но, к сожалению, занимает вечность. Первоначальный сценарий занимает 8,2 секунды, чтобы выполнить 20 000 шагов, а замена его пользовательским логарифмическим правдоподобием увеличивает это время до 6284,6 секунды. Мне нужно просто попытаться удалить выбросы раньше времени. - person pstjohn; 27.10.2014
comment
Для меня это увеличивает время работы только на 30%, что того стоит для 10-кратного увеличения точности. Вот блокнот, в котором они соревнуются. - person Abraham D Flaxman; 27.10.2014