Poisson-gamma regression


#1

Hi all,

I am new to Bayesian inference via Edward. I am implementing a Poisson-gamma regression in Edward. I compiled my Edward’s model. However, I have been unsuccessful with running Variational Bayes inference in Edward. Here is a snippet of code with the error that I’m getting:

N = 120
M = 5

x_train = np.reshape(z[:N], (N, 1))
y_train = np.reshape(w[:N], (N, 1))
t_train = np.reshape(range(N), (N, 1))
s_train = np.array([[12.0], [36.0], [60.0], [84.0], [108.0]])

X = tf.placeholder(tf.float32, [N, 1])
t = tf.placeholder(tf.float32, [N, 1])
s = tf.placeholder(tf.float32, [M, 1])

# Define priors
theta = Normal(mu=tf.zeros(1), sigma=tf.ones(1))
theta_star = Normal(mu=tf.zeros(1), sigma=tf.ones(1))
nu0 = Normal(mu=tf.ones(1)*1.5, sigma=tf.ones(1))
nu1 = Normal(mu=tf.ones(1)*2.0, sigma=tf.ones(1))

alpha_ = tf.exp(nu0)
beta_ = tf.exp(nu1)
gamma = Gamma(alpha=tf.ones([M, 1])*alpha_, beta=tf.ones([M, 1])*beta_)


def Kernel(t, s):
    mat = []
    for i in range(N):
        mat += [[]]
        for j in range(M):
            mat[i] += [1.2*tf.exp(-0.25*tf.pow(t[i] - s[j], 2))]
            
        mat[i] = tf.stack(mat[i])
    return tf.reshape(tf.stack(mat), [N, M])

# Poisson intensity
mu = theta*X + theta_star*tf.matmul(Kernel(t, s), gamma)

# Define likelihood
y = Poisson(lam=mu, value=tf.zeros_like(mu))


qtheta = Normal(mu=tf.Variable(tf.random_normal([1])),
            sigma=tf.nn.softplus(tf.Variable(tf.random_normal([1]))))
qtheta_star = Normal(mu=tf.Variable(tf.random_normal([1])),
            sigma=tf.nn.softplus(tf.Variable(tf.random_normal([1]))))

# qnu0 = Normal(mu=tf.Variable(tf.random_normal([1])), sigma=tf.nn.softplus(tf.Variable(tf.random_normal([1]))))
# qnu1 = Normal(mu=tf.Variable(tf.random_normal([1])), sigma=tf.nn.softplus(tf.Variable(tf.random_normal([1]))))

# qgamma = Gamma(alpha=tf.nn.softplus(tf.Variable(tf.random_normal([M, 1]))), beta=tf.nn.softplus(tf.Variable(tf.random_normal([M, 1]))))


inference = ed.KLqp({theta: qtheta, theta_star: qtheta_star}, data={X: x_train, t: t_train, s: s_train, y: y_train})
---------------------------------------------------------------------------
InvalidArgumentError                      Traceback (most recent call last)
<ipython-input-48-5b104cd13dec> in <module>()
----> 1 inference = ed.KLqp({theta: qtheta, theta_star: qtheta_star}, data={X: x_train, t: t_train, s: s_train, y: y_train})

/usr/local/lib/python2.7/site-packages/edward/inferences/klqp.pyc in __init__(self, *args, **kwargs)
     52   """
     53   def __init__(self, *args, **kwargs):
---> 54     super(KLqp, self).__init__(*args, **kwargs)
     55 
     56   def initialize(self, n_samples=1, kl_scaling=None, *args, **kwargs):

/usr/local/lib/python2.7/site-packages/edward/inferences/variational_inference.pyc in __init__(self, *args, **kwargs)
     21   """
     22   def __init__(self, *args, **kwargs):
---> 23     super(VariationalInference, self).__init__(*args, **kwargs)
     24 
     25   def initialize(self, optimizer=None, var_list=None, use_prettytensor=False,

/usr/local/lib/python2.7/site-packages/edward/inferences/inference.pyc in __init__(self, latent_vars, data)
    121           var = tf.Variable(ph, trainable=False, collections=[])
    122           self.data[key] = var
--> 123           sess.run(var.initializer, {ph: value})
    124         elif isinstance(value, np.number):
    125           if np.issubdtype(value.dtype, np.float):

/usr/local/lib/python2.7/site-packages/tensorflow/python/client/session.pyc in run(self, fetches, feed_dict, options, run_metadata)
    765     try:
    766       result = self._run(None, fetches, feed_dict, options_ptr,
--> 767                          run_metadata_ptr)
    768       if run_metadata:
    769         proto_data = tf_session.TF_GetBuffer(run_metadata_ptr)

/usr/local/lib/python2.7/site-packages/tensorflow/python/client/session.pyc in _run(self, handle, fetches, feed_dict, options, run_metadata)
    963     if final_fetches or final_targets:
    964       results = self._do_run(handle, final_targets, final_fetches,
--> 965                              feed_dict_string, options, run_metadata)
    966     else:
    967       results = []

/usr/local/lib/python2.7/site-packages/tensorflow/python/client/session.pyc in _do_run(self, handle, target_list, fetch_list, feed_dict, options, run_metadata)
   1013     if handle is None:
   1014       return self._do_call(_run_fn, self._session, feed_dict, fetch_list,
-> 1015                            target_list, options, run_metadata)
   1016     else:
   1017       return self._do_call(_prun_fn, self._session, handle, feed_dict,

/usr/local/lib/python2.7/site-packages/tensorflow/python/client/session.pyc in _do_call(self, fn, *args)
   1033         except KeyError:
   1034           pass
-> 1035       raise type(e)(node_def, op, message)
   1036 
   1037   def _extend_graph(self):

InvalidArgumentError: Node 'inference_15065146640/0/Variable_24/read': Unknown input node 'Variable_24'

#2

I swapped your x_train, y_train, t_train lines with the following in order to run the script

x_train = np.zeros((N, 1))
y_train = np.zeros((N, 1))
t_train = np.zeros((N, 1))

The script terminated successfully.

Could you try it again? My guess is that the iPython kernel was interacting strangely with TensorFlow sessions, so you need only restart the kernel.


#3

Hi Dustin,

I restarted the kernel, and I swapped the few lines of code you suggested. The script ran successfully.


#4

Hi Dustin,

I realize the problem I am having is with the input tensors shapes and types. x_train, y_train, and t_train are integer matrices. The input data looks like this:

u = np.reshape(z, (120, 1))

u = u.astype(float)

w = np.reshape(w, (120, 1))

z = np.reshape(range(121)[1:], (120, 1))

z =z.astype(float)

u[:5],w[:5], z[:5]

(array([[ 1.],
        [ 2.],
        [ 3.],
        [ 4.],
        [ 5.]]), array([[4],
        [5],
        [9],
        [6],
        [8]]), array([[ 1.],
        [ 2.],
        [ 3.],
        [ 4.],
        [ 5.]]))

Here is the code:

N = 120
M = 5

x_train = u
y_train = w
t_train = z
s_train = np.array([[12.00], [36.00], [60.00], [84.00], [108.00]])

X = tf.placeholder(tf.float32, [N, 1])
t = tf.placeholder(tf.float32, [N, 1])
s = tf.placeholder(tf.float32, [M, 1])

# Define priors
theta = Normal(mu=tf.zeros(1), sigma=tf.ones(1))
theta_star = Normal(mu=tf.zeros(1), sigma=tf.ones(1))
nu0 = Normal(mu=tf.ones(1)*1.5, sigma=tf.ones(1))
nu1 = Normal(mu=tf.ones(1)*2.0, sigma=tf.ones(1))

alpha_ = tf.exp(nu0)
beta_ = tf.exp(nu1)
gamma = Gamma(alpha=tf.ones([M, 1])*alpha_, beta=tf.ones([M, 1])*beta_)


def Kernel(t, s):
    mat = []
    for i in range(N):
        mat += [[]]
        for j in range(M):
            mat[i] += [1.2*tf.exp(-0.25*tf.pow(t[i] - s[j], 2))]
            
        mat[i] = tf.stack(mat[i])
    return tf.reshape(tf.stack(mat), [N, M])

# Poisson intensity
mu = theta*X + theta_star*tf.matmul(Kernel(t, s), gamma)

# Define likelihood
y = Poisson(lam=mu, value=tf.zeros_like(mu))


qtheta = Normal(mu=tf.Variable(tf.random_normal([1])),
            sigma=tf.nn.softplus(tf.Variable(tf.random_normal([1]))))
qtheta_star = Normal(mu=tf.Variable(tf.random_normal([1])),
            sigma=tf.nn.softplus(tf.Variable(tf.random_normal([1]))))

qnu0 = Normal(mu=tf.Variable(tf.random_normal([1])), sigma=tf.nn.softplus(tf.Variable(tf.random_normal([1]))))
qnu1 = Normal(mu=tf.Variable(tf.random_normal([1])), sigma=tf.nn.softplus(tf.Variable(tf.random_normal([1]))))
qgamma = Gamma(alpha=tf.nn.softplus(tf.Variable(tf.random_normal([M, 1]))), beta=tf.nn.softplus(tf.Variable(tf.random_normal([M, 1]))))

inference = ed.KLqp({theta: qtheta, theta_star: qtheta_star, nu0: qnu0, nu1: qnu1, gamma: qgamma}, data={X: x_train, t: t_train, s: s_train, y: y_train})

inference.run(n_samples=5, n_iter=250)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-34-ec1ea2c18385> in <module>()
----> 1 inference.run(n_samples=5, n_iter=250)

/usr/local/lib/python2.7/site-packages/edward/inferences/inference.pyc in run(self, variables, use_coordinator, *args, **kwargs)
    188       Passed into ``initialize``.
    189     """
--> 190     self.initialize(*args, **kwargs)
    191 
    192     if variables is None:

/usr/local/lib/python2.7/site-packages/edward/inferences/klqp.pyc in initialize(self, n_samples, kl_scaling, *args, **kwargs)
     76     self.n_samples = n_samples
     77     self.kl_scaling = kl_scaling
---> 78     return super(KLqp, self).initialize(*args, **kwargs)
     79 
     80   def build_loss_and_gradients(self, var_list):

/usr/local/lib/python2.7/site-packages/edward/inferences/variational_inference.pyc in initialize(self, optimizer, var_list, use_prettytensor, *args, **kwargs)
     62       var_list = list(var_list)
     63 
---> 64     self.loss, grads_and_vars = self.build_loss_and_gradients(var_list)
     65 
     66     if optimizer is None:

/usr/local/lib/python2.7/site-packages/edward/inferences/klqp.pyc in build_loss_and_gradients(self, var_list)
    124       #    return build_score_entropy_loss_and_gradients(self, var_list)
    125       else:
--> 126         return build_score_loss_and_gradients(self, var_list)
    127 
    128 

/usr/local/lib/python2.7/site-packages/edward/inferences/klqp.pyc in build_score_loss_and_gradients(inference, var_list)
    531     for x in six.iterkeys(inference.data):
    532       if isinstance(x, RandomVariable):
--> 533         x_copy = copy(x, dict_swap, scope=scope)
    534         p_log_prob[s] += tf.reduce_sum(
    535             inference.scale.get(x, 1.0) * x_copy.log_prob(dict_swap[x]))

/usr/local/lib/python2.7/site-packages/edward/util/random_variables.pyc in copy(org_instance, dict_swap, scope, replace_itself, copy_q)
    174     kwargs['name'] = new_name
    175     # Create new random variable with copied arguments.
--> 176     new_rv = rv.__class__(*args, **kwargs)
    177     return new_rv
    178   elif isinstance(org_instance, tf.Tensor):

/usr/local/lib/python2.7/site-packages/edward/models/random_variable.pyc in __init__(self, *args, **kwargs)
     70 
     71     if value is not None:
---> 72       t_value = tf.convert_to_tensor(value, self.dtype)
     73       expected_shape = (self.get_batch_shape().as_list() +
     74                         self.get_event_shape().as_list())

/usr/local/lib/python2.7/site-packages/tensorflow/python/framework/ops.pyc in convert_to_tensor(value, dtype, name, preferred_dtype)
    635       name=name,
    636       preferred_dtype=preferred_dtype,
--> 637       as_ref=False)
    638 
    639 

/usr/local/lib/python2.7/site-packages/tensorflow/python/framework/ops.pyc in internal_convert_to_tensor(value, dtype, name, as_ref, preferred_dtype)
    700 
    701         if ret is None:
--> 702           ret = conversion_func(value, dtype=dtype, name=name, as_ref=as_ref)
    703 
    704         if ret is NotImplemented:

/usr/local/lib/python2.7/site-packages/tensorflow/python/ops/variables.pyc in _TensorConversionFunction(v, dtype, name, as_ref)
    647       raise ValueError(
    648           "Incompatible type conversion requested to type '%s' for variable "
--> 649           "of type '%s'" % (dtype.name, v.dtype.name))
    650     if as_ref:
    651       return v._ref()  # pylint: disable=protected-access

ValueError: Incompatible type conversion requested to type 'float32' for variable of type 'int32_ref'

#5

Good catch. This is because the Poisson random variable uses a float data type and not int:

>>> from edward.models import Poisson
>>> x = Poisson(lam=5.0, value=0)
>>> print(x)
RandomVariable("Poisson_1/", shape=(), dtype=float32)

This means any data binded to Poisson random variables must be float and not int. Similarly, you defined your placeholders for t and s to be of type float, although you’re binding it to data of type int.

I added a Github issue to raise a more informative error message if there’s a type mismatch (https://github.com/blei-lab/edward/issues/560).

As an aside, I’m not actually sure why Poisson uses float and not int. It is a result of this [line]
(https://github.com/tensorflow/tensorflow/blob/master/tensorflow/contrib/distributions/python/ops/poisson.py#L86) in TensorFlow; maybe because of unbounded support?


#6

Hi Dustin,

I’m a postdoc in Dr. Barbara Engelhardt’s lab at Princeton. I’m working on a spatial regression for a marked point Poisson process with the Poisson mean that is a sum of fixed effects part and a doubly stochastic process represented by a convolution of a Gaussian kernel and gamma process. The placeholders t and s are inputs to the kernel and gamma is gamma process in the code. I am getting nan’s for the loss when I ran the variational inference. I’m going to try to run a simple Poisson regression to understand how Poisson distribution works in Edward. You have helped me a lot.


#7

Great! Another dimension worth experimenting is to not simplify the model but to simplify the inference. That is, I recommend just trying MAP. Suchi Saria’s lab has produce some excellent healthcare applications using Bayesian nonparametric models and MAP.