Simple Metropolis Hastings Model


#1

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


#2

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.