Commit 2ff43803 authored by Azat Garifullin's avatar Azat Garifullin
Browse files

metropolis-gibbs with noise

parent cf4e4c18
No preview for this file type
No preview for this file type
......@@ -10,36 +10,29 @@ from matplotlib import pyplot as plt
import random_fields as rf
def prior_sampler(_):
return [
np.random.uniform(0.5, 2.5),
np.random.uniform(0.5, 5),
np.random.uniform(0.5, 5)
]
def make_neg_log_posterior(observation):
def impl(state, proposal_field, field_cov):
return rf.matern_parameter_posterior_nl(observation, proposal_field, field_cov)
return impl
def main():
with open('./data/field.pickle', 'rb') as f:
xs, ys, field, nu, sigma, l = pickle.load(f)
print("Pre-evaluating distance matrix")
dist_mtx = rf.dist_matrix(xs, ys)
flat_gt_field = np.ravel(field)
noise_sigma = 1e-2
noise = np.random.normal(scale=noise_sigma, size=flat_gt_field.shape)
flat_field = np.ravel(field)
initial_field = flat_gt_field + noise
initial_theta = np.array([nu, sigma, l]) + 1e-2
initial = (prior_sampler(0), flat_field)
field_sampler = rf.make_matern_sampler(dist_mtx)
neg_log_posterior = make_neg_log_posterior(flat_field)
proposal_cov = 1e-4 * np.eye(3)
lower = np.array([0.5, 0.5, 0.5])
upper = np.array([3, 3, 3])
print("Starting MCMC")
chain_res = rf.gibbs_metropolis(initial, prior_sampler,
neg_log_posterior, field_sampler)
chain_res = rf.gibbs_within_metropolis(
initial_theta, flat_gt_field, noise_sigma,
proposal_cov, lower, upper,
dist_mtx, initial_field=initial_field)
with open('./data/matern_chain.pickle', 'wb') as f:
pickle.dump(chain_res, f)
......@@ -49,4 +42,4 @@ if __name__ == '__main__':
random.seed(42)
np.random.seed(42)
main()
\ No newline at end of file
main()
figs/dist_mtx.png

4.89 MB | W: | H:

figs/dist_mtx.png

477 KB | W: | H:

figs/dist_mtx.png
figs/dist_mtx.png
figs/dist_mtx.png
figs/dist_mtx.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_05_0.png

25.3 KB | W: | H:

figs/k_smooth_05_0.png

17.9 KB | W: | H:

figs/k_smooth_05_0.png
figs/k_smooth_05_0.png
figs/k_smooth_05_0.png
figs/k_smooth_05_0.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_05_1.png

25.1 KB | W: | H:

figs/k_smooth_05_1.png

17.2 KB | W: | H:

figs/k_smooth_05_1.png
figs/k_smooth_05_1.png
figs/k_smooth_05_1.png
figs/k_smooth_05_1.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_05_2.png

25.2 KB | W: | H:

figs/k_smooth_05_2.png

17.7 KB | W: | H:

figs/k_smooth_05_2.png
figs/k_smooth_05_2.png
figs/k_smooth_05_2.png
figs/k_smooth_05_2.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_05_cov_mtx.png

268 KB | W: | H:

figs/k_smooth_05_cov_mtx.png

183 KB | W: | H:

figs/k_smooth_05_cov_mtx.png
figs/k_smooth_05_cov_mtx.png
figs/k_smooth_05_cov_mtx.png
figs/k_smooth_05_cov_mtx.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_1_0.png

22.2 KB | W: | H:

figs/k_smooth_1_0.png

17.3 KB | W: | H:

figs/k_smooth_1_0.png
figs/k_smooth_1_0.png
figs/k_smooth_1_0.png
figs/k_smooth_1_0.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_1_1.png

21.6 KB | W: | H:

figs/k_smooth_1_1.png

19.1 KB | W: | H:

figs/k_smooth_1_1.png
figs/k_smooth_1_1.png
figs/k_smooth_1_1.png
figs/k_smooth_1_1.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_1_2.png

22 KB | W: | H:

figs/k_smooth_1_2.png

20.4 KB | W: | H:

figs/k_smooth_1_2.png
figs/k_smooth_1_2.png
figs/k_smooth_1_2.png
figs/k_smooth_1_2.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_1_cov_mtx.png

278 KB | W: | H:

figs/k_smooth_1_cov_mtx.png

191 KB | W: | H:

figs/k_smooth_1_cov_mtx.png
figs/k_smooth_1_cov_mtx.png
figs/k_smooth_1_cov_mtx.png
figs/k_smooth_1_cov_mtx.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_2_0.png

20.2 KB | W: | H:

figs/k_smooth_2_0.png

20.6 KB | W: | H:

figs/k_smooth_2_0.png
figs/k_smooth_2_0.png
figs/k_smooth_2_0.png
figs/k_smooth_2_0.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_2_1.png

20 KB | W: | H:

figs/k_smooth_2_1.png

18 KB | W: | H:

figs/k_smooth_2_1.png
figs/k_smooth_2_1.png
figs/k_smooth_2_1.png
figs/k_smooth_2_1.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_2_2.png

20.5 KB | W: | H:

figs/k_smooth_2_2.png

18.6 KB | W: | H:

figs/k_smooth_2_2.png
figs/k_smooth_2_2.png
figs/k_smooth_2_2.png
figs/k_smooth_2_2.png
  • 2-up
  • Swipe
  • Onion skin
figs/k_smooth_2_cov_mtx.png

290 KB | W: | H:

figs/k_smooth_2_cov_mtx.png

191 KB | W: | H:

figs/k_smooth_2_cov_mtx.png
figs/k_smooth_2_cov_mtx.png
figs/k_smooth_2_cov_mtx.png
figs/k_smooth_2_cov_mtx.png
  • 2-up
  • Swipe
  • Onion skin
......@@ -41,77 +41,118 @@ def dist_matrix(xs, ys):
return dist_mtx
def make_matern_sampler(dist_mtx, n_sigma=eps()):
def sampler_impl(state, mean):
nu, sigma, l = state
def sample_matern_cond_posterior(mean_field, cov_mtx, noise_sigma):
i_mtx = np.eye(cov_mtx.shape[0])
nc_inv = np.linalg.inv(cov_mtx + (noise_sigma**2) * i_mtx)
cnc = np.matmul(cov_mtx, nc_inv)
cov_mtx = matern_cov_dm(dist_mtx, nu, sigma, l)
mu = np.matmul(cnc, mean_field)
cov = cov_mtx - np.matmul(cnc, cov_mtx)
chol_cov = np.linalg.cholesky(cov)
n = np.random.normal(size=(cov_mtx.shape[0],))
nc_inv = np.linalg.inv(cov_mtx + n_sigma * np.eye(cov_mtx.shape[0]))
cnc = np.matmul(cov_mtx, nc_inv)
return mu + np.matmul(chol_cov, n)
mu = np.matmul(cnc, mean)
cov = cov_mtx - np.matmul(cnc, cov_mtx)
chol_cov = np.linalg.cholesky(cov)
n = np.random.normal(size=(cov_mtx.shape[0],))
return mu + np.matmul(chol_cov, n), cov
return sampler_impl
def matern_parameter_posterior_nl(observed, field, field_cov, n_sigma=0):
def matern_parameter_nlpost(observed, field, field_cov, noise_sigma):
inv_noise_cov = np.eye(np.size(field)) # * (1.0 / noise_sigma**2)
field_diff = observed - field
nll = np.dot(field_diff, field_diff)
nll = np.matmul(np.matmul(field_diff, inv_noise_cov), field_diff)
field_inv_cov = np.linalg.inv(field_cov)
field_nl_prior = np.matmul(np.matmul(field.T, field_inv_cov), field)
chol_cov = np.linalg.cholesky(field_cov)
tril = chol_cov[np.tril_indices(chol_cov.shape[0])]
log_det_prior = 2 * np.sum(tril)
cc_diag = np.diag(chol_cov)
log_det_prior = 2 * np.sum(np.log(cc_diag))
return nll + field_nl_prior + 0.5 * (log_det_prior + chol_cov.shape[0] * np.log(2 * np.pi))
return nll + field_nl_prior + log_det_prior
def gibbs_metropolis(initial, proposal_sampler, f_nlp, field_sampler,
n_samples=1000, burn_in=100, n_save=500):
def gibbs_within_metropolis(initial_theta, observed_field, noise_sigma,
proposal_cov, lower, upper,
dist_mtx, initial_field=None,
adaptation_start=1,
n_samples=10000, n_save=1000):
chain = []
nlls = []
p_accs = []
state, observed_field = initial
field, field_cov = field_sampler(state, observed_field)
state_nll = f_nlp(state, field, field_cov)
for i in range(n_samples):
proposal_field, field_cov = field_sampler(state, observed_field)
proposal_state = proposal_sampler(state)
proposal_nll = f_nlp(proposal_state, proposal_field, field_cov)
acceptance = np.exp(state_nll - proposal_nll)
if np.random.random() < acceptance:
state = np.log(initial_theta)
field = initial_field
if field is None:
field = np.copy(observed_field)
nu, sigma, l = initial_theta
field_cov = matern_cov_dm(dist_mtx, nu, sigma, l)
state_nlp = matern_parameter_nlpost(observed_field, field,
field_cov, noise_sigma)
p_mean_intext = np.copy(state)
sd = 5.76 / state.size
eps = np.finfo(np.float32).eps
eps_eye = eps * np.eye(state.size)
proposal_cov += eps_eye
if adaptation_start is None:
adaptation_start = n_samples
for step in range(n_samples):
pr_chol = np.linalg.cholesky(proposal_cov)
dstate = np.matmul(pr_chol, np.random.normal(size=state.shape))
proposal_state = state + dstate
proposal_theta = np.exp(proposal_state)
up_lower = np.all(proposal_theta >= lower)
lo_upper = np.all(proposal_theta <= upper)
if not (lo_upper and up_lower):
print(
"Gibbs-Metropolis [{}]: rejected with out of bounds".format(step))
continue
nu, sigma, l = proposal_theta
proposal_field_cov = matern_cov_dm(dist_mtx, nu, sigma, l)
proposal_field = sample_matern_cond_posterior(
observed_field, proposal_field_cov, noise_sigma)
proposal_nlp = matern_parameter_nlpost(
observed_field, proposal_field, proposal_field_cov, noise_sigma)
# print("Proposal theta: {}".format(proposal_theta))
# print("Proposal field diff: {}".format(
# np.mean(np.abs(proposal_field - observed_field))))
# print("State neg-log posterior: {}".format(state_nlp))
# print("Proposal neg-log posterior: {}".format(proposal_nlp))
acceptance = np.exp(state_nlp - proposal_nlp)
accepted = np.random.random() < acceptance
if accepted:
state = proposal_state
field = proposal_field
state_nll = proposal_nll
print("Gibbs-Metropolis [{}]: acceptance with p={}".format(i, acceptance))
state_nlp = proposal_nlp
print(
"Gibbs-Metropolis [{}]: acceptance with p={}".format(step, acceptance))
else:
print("Gibbs-Metropolis [{}]: rejected with p={}".format(i, acceptance))
print(
"Gibbs-Metropolis [{}]: rejected with p={}".format(step, acceptance))
if i > (n_samples - n_save):
if step > (n_samples - n_save):
chain.append((state, field))
p_accs.append(acceptance)
nlls.append(state_nll)
return chain, p_accs, nlls
n_mean_intext = p_mean_intext + (state - p_mean_intext) / (step + 1)
if step >= adaptation_start:
base = (step - 1) / step * proposal_cov
pmmt = np.outer(p_mean_intext, p_mean_intext) * step
nmmt = np.outer(n_mean_intext, n_mean_intext) * (step + 1)
sst = np.outer(state, state)
upd_cov = pmmt - nmmt + sst + eps_eye
proposal_cov = base + sd / step * upd_cov
# print("[{0}] Proposal cov trace: {1}".format(
# step, np.trace(proposal_cov)))
p_accs.append((acceptance, accepted))
nlls.append(state_nlp)
return chain, p_accs, nlls
......@@ -12,9 +12,9 @@ import random_fields as rf
def main():
w, h = 64, 64
xs = np.linspace(0, 5, w)
ys = np.linspace(0, 5, h)
w, h = 32, 32
xs = np.linspace(0, 2, w)
ys = np.linspace(0, 2, h)
n_pts = np.size(xs) * np.size(ys)
......@@ -22,7 +22,8 @@ def main():
dist_mtx = rf.dist_matrix(xs, ys)
print(dist_mtx.min())
dist_mtx_img = (dist_mtx - dist_mtx.min()) / (dist_mtx.max() - dist_mtx.min()) * 255
dist_mtx_img = (dist_mtx - dist_mtx.min()) / \
(dist_mtx.max() - dist_mtx.min()) * 255
cv2.imwrite('./figs/dist_mtx.png', dist_mtx_img)
ns = [
......@@ -62,6 +63,7 @@ def sample_and_save_fig(k_mtx, ns, w, h, prefix):
plt.figure()
plt.imshow(field_sample)
plt.colorbar()
plt.savefig('./figs/{}_{}'.format(prefix, i))
plt.close()
......@@ -72,4 +74,4 @@ if __name__ == '__main__':
random.seed(42)
np.random.seed(42)
main()
\ No newline at end of file
main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment