# Simple Metropolis Hastings Model

I am trying to port some of my pymc models over to Edward to get a better feeling for what it can do, and am running into trouble that is probably simple but that I can’t figure out. Here is the code:

y = [0.6246, 0.5864, 0.5614, 0.5385, 0.5161, 0.4968, 0.4787, 0.4641, 0.4492, 0.4375, 0.427, 0.4188, 0.414, 0.4075, 0.4039, 0.4007, 0.3976, 0.3966, 0.3944, 0.393, 0.3919,
0.3907, 0.3895, 0.3887, 0.3874, 0.3873, 0.3864, 0.3869, 0.3857, 0.3849, 0.3844, 0.3844, 0.3833, 0.3829, 0.3837,
0.3829, 0.3824, 0.3814, 0.3816, 0.3819, 0.3795, 0.3801, 0.3806, 0.379, 0.3792, 0.3794, 0.3782, 0.3777, 0.3769,
0.3764, 0.3762, 0.3755, 0.3748, 0.3752, 0.3744, 0.374, 0.375, 0.3737, 0.3731, 0.3727, 0.3723, 0.3716, 0.3708,
0.3722, 0.371, 0.371, 0.3692, 0.3705, 0.3697, 0.369]

x = np.linspace(0, len(y)*16, len(y))

def rounded_hinge_edward(x, x0, a, b0, b1, delta):
# a + b0 * (x - x0) + (b1 - b0) * delta * tf.log(1.0 + tf.exp((x - x0) / delta))
# (x - x0)
x_minus_x0 = tf.subtract(x, x0)
# b0 * (x - x0)
second_term = tf.multiply(b0, x_minus_x0)
# a + b0 * (x - x0)

# (x - x0) / delta
inside_exp = tf.divide(x_minus_x0, delta)

# np.exp((x - x0) / delta)
exp_term = tf.exp(inside_exp)

# delta * tf.log(1.0 + tf.exp((x - x0) / delta)
delta_times_log = tf.multiply(delta, log_term)

# b - b0
b_minus_b0 = tf.subtract(b1, b0)
third_term = tf.multiply(b_minus_b0, delta_times_log)

def rounded_hinge_edward_model(x, f):
x0 = Uniform(name='x0', low=-1000., high=1000.)
a = Uniform(name='a', low=1E-9, high=1000.)
b1 = Uniform(name='b1', low=-1000., high=1000.)
b0 = Uniform(name='b0', low=-1000., high=1000.)
delta = Uniform(name='delta', low=1E-9, high=1000.)
sigma = Uniform(name='sigma', low=1E-9, high=1.0)

X = tf.placeholder(tf.float32, len(x), name='X')
one = tf.Variable(1.0, tf.float32)
model = rounded_hinge_edward(X, x0, a, b0, b1, delta)
y = Normal(loc=model, scale=1.0)
print(model, sigma, y)

# INFERENCE
T = 75000  # number of posterior samples

q_x0 = Empirical(params=tf.Variable(tf.random_uniform([T], minval=-1000., maxval=1000.)))
q_a = Empirical(params=tf.Variable(tf.random_uniform([T], minval=1E-9, maxval=1000.)))
q_b1 = Empirical(params=tf.Variable(tf.random_uniform([T], minval=-1000., maxval=1000.)))
q_b0 = Empirical(params=tf.Variable(tf.random_uniform([T], minval=-1000., maxval=1000.)))
q_delta = Empirical(params=tf.Variable(tf.random_uniform([T], minval=1E-9, maxval=1000.)))
q_sigma = Empirical(params=tf.Variable(tf.random_uniform([T], minval=1E-9, maxval=1.0)))

# Inference arguments
latent_vars = {x0: q_x0, a: q_a, b1: q_b1, b0: q_b0, delta: q_delta, sigma: q_sigma}
data={X: x, y: f}

proposal_x0 = Uniform(name='proposal_x0', low=1E-9, high=1000.)
proposal_a = Uniform(name='proposal_a', low=1E-9, high=1000.)
proposal_b1 = Uniform(name='proposal_b1', low=-1000., high=1000.)
proposal_b0 = Uniform(name='proposal_b0', low=-1000., high=1000.)
proposal_delta = Uniform(name='proposal_delta', low=1E-9, high=1000.)
proposal_sigma = Uniform(name='proposal_sigma', low=1E-9, high=1.0)
proposals = {x0: proposal_x0, a: proposal_a, b1: proposal_b1, b0: proposal_b0, delta: proposal_delta, sigma: proposal_sigma}

inference = ed.MetropolisHastings(latent_vars, proposals, data=data, )
inference.run()

mean_x0 = tf.reduce_mean(q_x0.params[60000:]).eval()
mean_sigma = tf.reduce_mean(q_sigma.params[60000:]).eval()
mean_a = tf.reduce_mean(q_a.params[60000:]).eval()
mean_b1 = tf.reduce_mean(q_b1.params[60000:]).eval()
mean_b0 = tf.reduce_mean(q_b0.params[60000:]).eval()
mean_delta = tf.reduce_mean(q_delta.params[60000:]).eval()
print(mean_x0, mean_a, mean_b1, mean_b0, mean_delta, mean_sigma)
plt.plot(x, x*mean_b0+mean_a)
plt.plot(x, f)

rounded_hinge_edward_model(x, np.asarray(y))

I don’t get a good fit and don’t get a very high acceptance probability. I am sure I am doing something simple wrong as I am still learning both Tensorflow and Edward, but some help would be greatly appreciated!

James

Can you post you original (or corresponding) pymc code? Besides perhaps the very wide ranges in the priors and the proposals, I don’t see anything that stands out. Also, I would first try with simpler model – sub some of the variables with constants and see how the inference behaves.
TF tip: you can actually code rounded_hinge_edward like this:

def rounded_hinge_edward(x, x0, a, b0, b1, delta):
return a + b0 * (x - x0) + (b1 - b0) * delta * tf.log(1.0 + tf.exp((x - x0) / delta))

It will produce tensor.