Edward modeling to an artificial data

Overview

On the article below, I switched method on Edward model from variational method to Hamiltonian Monte Carlo.

• Hamiltonian Monte Carlo on TensorFlow and Edward

As an another example, I'll try same thing to the model of the following article.

• Edward modeling to artificial data with random effects

In a nutshell, I'll make model for an artificial data and get sample by Hamiltonian Monte Carlo.

Data

Here, I’ll use the following artificial data from the article, Edward modeling to artificial data with random effects.

import numpy as np
np.random.seed(49)

x_1 = list(range(1, 100))
x_2 = list(range(1, 100))
x_3 = list(range(1, 100))

y_1 = [np.random.normal(3.0 * i + 4, 4.0) for i in x_1]
y_2 = [np.random.normal(4.0 * i + 5, 4.0) for i in x_2]
y_3 = [np.random.normal(2.0 * i + 2, 20.0) for i in x_3]

x_1_sign = [1 for i in x_1]
x_2_sign = [2 for i in x_2]
x_3_sign = [3 for i in x_3]

X = x_1 + x_2 + x_3
Y = y_1 + y_2 + y_3
GID = x_1_sign + x_2_sign + x_3_sign

By visualization, we can check the data behavior.

It contains three types of data dependent on GID.

Edward Modeling

The original code is on the article, Edward modeling to artificial data with random effects. By changing some points, we can switch the method.
On the importing part, you need to import Empirical. And on the inference part, we also need to set the Empirical responding to the parameters. Be careful about the dimension and the values of Empirical.

import tensorflow as tf
import edward as ed
from edward.models import Normal, Empirical

# data
X_data = np.array(X).reshape([len(X), 1])
g_n = len(set(GID))
GID_data = np.array([x-1 for x in GID])

# placeholder
x_ph = tf.placeholder(tf.float32, [None, 1])
gid_ph = tf.placeholder(tf.int32, None)

# fixed effect
a_common = Normal(loc=tf.zeros(1), scale=tf.ones(1))
b_common = Normal(loc=tf.zeros(1), scale=tf.ones(1))

# random effect
random_effect = Normal(loc=tf.zeros([g_n]), scale=tf.ones([g_n]))

# model
y_hat = ed.dot(x_ph, b_common) + a_common + tf.gather(random_effect, gid_ph)
output = Normal(loc=y_hat, scale=tf.ones(1))

# inference
T = 10000
q_a_common = Empirical(params=tf.Variable(tf.random_normal([T, 1])))
q_b_common = Empirical(params=tf.Variable(tf.random_normal([T, 1])))
q_random_effect = Empirical(params=tf.Variable(tf.random_normal([T, g_n])))

inference = ed.HMC({a_common: q_a_common, b_common: q_b_common, random_effect: q_random_effect},
                   {x_ph: X_data, output: np.array(Y), gid_ph: GID_data})
inference.run(step_size=1e-3)

10000/10000 [100%] ██████████████████████████████ Elapsed: 24s | Acceptance Rate: 0.920

We can get sampled points and check those on histogram easily.

import matplotlib.pyplot as plt
plt.subplot(3,1,1)
plt.title("b")
plt.hist(q_b_common.sample(1000).eval())
plt.subplot(3,1,2)
plt.title("a")
plt.hist(q_a_common.sample(1000).eval())
plt.subplot(3,1,3)
plt.title("random_effect")
plt.hist(q_random_effect.sample(1000).eval())
plt.show()

On the same manner as the original article, we can generated the predicted points.

n_samples = 100
prob_lst = []
samples = []
b_samples = []
a_samples = []
random_effect_samples = []
for _ in range(n_samples):
    b_samp = q_b_common.sample()
    a_samp = q_a_common.sample()
    random_samp = q_random_effect.sample()
    b_samples.append(b_samp)
    a_samples.append(a_samp)
    random_effect_samples.append(random_samp)
    prob = ed.dot(tf.cast(X_data, tf.float32), b_samp) + a_samp + tf.gather(random_samp, GID_data)

    prob_lst.append(prob.eval())

for i,val in enumerate(prob_lst):
    if i == 0:
        added = val
    else:
        added += val

added_mean = added / len(prob_lst)

By visualizing it, we can see the three lines. On this model, I reflected the GID effect only on the interception. For appropriate modeling, we basically need to think about the slop on this case.

plt.scatter(X_data, added_mean)
plt.show()