Реализация алгоритма Гиббса в Python ⇐ Python

[Anonymous]

Я пытаюсь реализовать подход, использованный в статье Ф. Руджери и С. Сиваганесан «О моделировании точек изменения в неоднородных пуассоновских процессах». В своей статье они определяют полные условия


В статье говорится, что значения бета можно определить с помощью выборки Гиббса. n — количество точек данных, где точки данных — это значения t_i. t_0 = 0, а t_n+1 — это значение, немного большее, чем t_n. Я выбрал это максимальное значение t_n + 1. Другие константы имеют (начальные) значения:

Код: Выделить всё

alpha, delta, rho, tau, mu, phi = 0.01, 0.01, 0.01, 0.01, 0, 0

Обратите внимание, что бета-версии пропорциональны, а не равны, что означает, что мне нужно определить константу нормализации. Вот в этом для меня и заключается проблема.
Я попробовал следующую реализацию на Python.

Код: Выделить всё

def M_sample(alpha, delta, n, t, beta):
shape = alpha + n
first_val = np.power(t[0], beta[0]) - np.power(0, beta[0])
last_val = np.power(np.ceil(t[-1]), beta[-1]) - np.power(t[-1], beta[-1])
scale = 1/(delta + sum(np.power(t[i], beta[i]) - np.power(t[i-1], beta[i]) for i in range(1, n-1)) + first_val + last_val)
return gamma.rvs(shape, scale=scale)

def sample_sigma2(rho, tau, n, beta, phi):
shape = rho + n/2
scale = tau + 0.5*sum(np.power(np.log(beta[i]) - phi, 2) for i in range(n+1))
#print("scale",scale)
return invgamma.rvs(a=shape, scale=scale)

def sample_phi(mu, tau, n, beta, sigma2):
tau2 = np.power(tau, 2)
mu1 = (mu*sigma2 + tau2*sum(np.log(beta[i]) for i in range(n+1)))/(sigma2 + (n+1)*tau2)
denom = (tau2*(n+1) + sigma2)
tau1_sq = (tau2*sigma2)/denom
return norm.rvs(mu1, np.sqrt(tau1_sq))

def unnormalized_posterior_beta_i(beta, M, t, phi, sigma2, i):
return np.power(t[i], beta) * np.exp(-M*(np.power(t[i], beta) - np.power(t[i-1], beta)) - (np.power(np.log(beta)- phi, 2)/(2*sigma2)))

def unnormalized_posterior_beta_0(beta, M, t, phi, sigma2):
return np.power(t[0], beta) * np.exp(-M*(np.power(t[0], beta) - np.power(0, beta)) - (np.power(np.log(beta)- phi, 2)/(2*sigma2)))

def unnormalized_posterior_beta_n(beta, M, t, phi, sigma2):
y = np.ceil(t[-1]+1)
return (beta**(-1)) * np.exp(-M*(np.power(y, beta) - np.power(t[-1], beta)) - (np.power(np.log(beta)- phi, 2)/(2*sigma2)))

def sample_beta(n, M, t, phi, sigma2, beta):
for i in range(1, n):
integral_i, _ = quad(unnormalized_posterior_beta_i, 0, 10, args=(M, t, phi, sigma2, i))
beta_i = unnormalized_posterior_beta_i(beta[i], M, t, phi, sigma2, i)
beta[i] = beta_i/integral_i

integral_0, _ = quad(unnormalized_posterior_beta_0, 0, 10, args=(M, t, phi, sigma2))
beta_0 = unnormalized_posterior_beta_0(beta[0], M, t, phi, sigma2)
beta[0] = beta_0/integral_0

integral_n, _ = quad(unnormalized_posterior_beta_n, 0, 10, args=(M, t, phi, sigma2))
beta_n = unnormalized_posterior_beta_n(beta[-1], M, t, phi, sigma2)
print(beta_n, integral_n)
beta[-1] = beta_n/integral_n

return beta

def gibbs_sampler(n, t, n_samples=5):
alpha, delta, rho, tau, mu, phi = 0.01, 0.01, 0.01, 0.01, 0, 0
M_samples = []
sigma2_samples = []
phi_samples = []
beta_samples = []

beta = np.ones(n+1)

for _ in range(n_samples):
M = M_sample(alpha, delta, n, t, beta)
sigma2 = sample_sigma2(rho, tau, n, beta, phi)
phi = sample_phi(mu, tau, n, beta, sigma2)
beta = sample_beta(n, M, t, phi, sigma2, beta)
print(beta)

M_samples.append(M)
sigma2_samples.append(sigma2)
phi_samples.append(phi)
beta_samples.append(beta)

return M_samples, sigma2_samples, phi_samples, beta_samples

Основная проблема заключается в том, что значения бета становятся слишком большими, особенно beta[-1]. С помощью линий интегрирования, которые я делаю, я должен вычислять нормализующую константу, которая должна привести к разумным значениям бета. Поэтому я предполагаю, что я ошибаюсь с константой нормализации. Где я допускаю ошибку при реализации алгоритма?

Подробнее здесь: https://stackoverflow.com/questions/790 ... -in-python