Tensorflow for Statisticians (2)

Generalized Linear Regression

Last updated on Apr 17, 2022 2 min read

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()
>>> model.add(layers.Dense(1, activation = "sigmoid"))
>>> model.compile(optimizer = tf.optimizers.Adam(0.05),
        loss = tf.keras.losses.BinaryCrossentropy())
>>> model.fit(data, Y, epochs = 10, batch_size = 1024)
>>> model.weights[0].numpy().flatten()
# array([0.5853649 , 0.12772334, 0.50310177, 0.40288347, 0.36750564,
#        0.45062354, 0.53615963, 0.31267476, 0.38474655, 0.8759412 ,
#        0.24433382, 0.31499726, 0.48656067, 0.91565496, 0.04293687,
#        0.4549262 , 0.10319919, 0.22904651, 0.19099398, 0.22391798],
#       dtype=float32)

technical tensorflow

Tensorflow for Statisticians (2)

Yishu Xue

Data Scientist / Coder / Novice Sprinter / Gym Enthusiast