# Tensorflow for Statisticians (2)

Generalized Linear Regression

Following my last post on using tensorflow for linear regression, in this post I am going to extend the scope to generalized linear models.

>>> import tensorflow as tf
>>> import tensorflow_probability as tfp
>>> tfd = tfp.distributions
>>> from tensorflow.keras import layers


Let us simulate some data first.

>>> tf.compat.v1.random.set_random_seed(1) # for reproduce
>>> data = tfd.Normal(loc = 0., scale = 1.).sample([10000, 20])
>>> data
# <tf.Tensor: id=22, shape=(1000000, 20), dtype=float32, numpy=
# array([[-1.1012203 ,  1.5457517 ,  0.383644  , ..., -0.02147584,
#         -0.31968883,  0.37332553],
#        [ 0.25279108,  0.6437664 ,  2.146308  , ..., -0.15791255,
#         -0.3566943 , -0.20044029],
#        [ 1.613107  ,  0.6796728 ,  0.08133233, ...,  1.9076501 ,
#         -0.7215866 ,  1.2887349 ],
#        ...,
#        [-0.34464952,  0.33008808,  1.6495917 , ...,  0.35497805,
#          0.8145206 ,  1.4482048 ],
#        [-0.32808402, -0.20314096, -2.0368671 , ..., -0.50538176,
#         -0.04394155,  0.03310239],
#        [-0.7274147 , -0.35093465,  0.8078732 , ..., -0.10988711,
#         -1.7783372 ,  0.15070562]], dtype=float32)>
>>> tf.compat.v1.random.set_random_seed(3)
>>> pars = tfd.Beta(2, 3).sample([20, 1])
>>> lp = tf.matmul(data, pars)
>>> Y = tfd.Bernoulli(lp).sample(1)
>>> Y
# <tf.Tensor: id=1579, shape=(1000000,), dtype=float32, numpy=array([0., 0., 1., ..., 1., 0., 0.], dtype=float32)>


Tensorflow-probability provides a glm module which enables easy fitting of a collection of generalized linear models such as logistic regression (tfp.glm.Bernoulli), probit regression (tfp.glm.BernoulliNormalCDF), Poisson regression (tfp.glm.Poisson), and many more. An one-line implementation of logistic regression is:

>>> w, linear_response, is_converged, num_iter = tfp.glm.fit(
model_matrix = data,
response = Y,
model = tfp.glm.Bernoulli())
>>> w
>>> w
# <tf.Tensor: id=554, shape=(20,), dtype=float32, numpy=
# array([0.58070004, 0.08244729, 0.48235574, 0.39969572, 0.37823585,
#        0.39500824, 0.5433686 , 0.27827388, 0.38219643, 0.8608786 ,
#        0.22414763, 0.34881216, 0.4582049 , 0.9232365 , 0.02907711,
#        0.47394714, 0.16831397, 0.26486722, 0.20210025, 0.17360945],
#       dtype=float32)>
>>> w.numpy() - pars.numpy().flatten()
# array([ 2.2941828e-03,  9.4745308e-04, -2.6564598e-03,  3.9328039e-03,
#         2.1443665e-03,  4.0888786e-04,  1.4651418e-03, -1.3917387e-03,
#        -2.0698011e-03,  1.9533038e-03,  7.5770915e-04, -2.2533536e-04,
#         3.8230419e-04,  4.8930049e-03, -1.8272493e-03, -6.6980720e-04,
#         2.6237667e-03,  1.7036200e-03,  2.6491731e-03, -3.5420060e-05],
#       dtype=float32)


Another approach is to use keras. We build a model with one dense layer and sigmoid activation function. Then we can use model.fit() to train the model on the dataset.

>>> model = tf.keras.Sequential()