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)
first_two_terms = tf.add(a, second_term)
# (x - x0) / delta
inside_exp = tf.divide(x_minus_x0, delta)
# np.exp((x - x0) / delta)
exp_term = tf.exp(inside_exp)
log_term = tf.log(tf.add(1.0, exp_term))
# 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)
return tf.add(first_two_terms, third_term)
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