How to code regression using maximum likelihood in Python using only numpy?

Question

I am trying to do Poisson regression from scratch (only using numpy). I have the following code, but the output makes no sense to me.

import numpy as np

First I generate some data:

# generate data
x1 = np.random.normal(10, 5, 1000)
x2 = np.random.normal(20, 5, 1000)
error = np.random.normal(0, 1, 1000)
y = 3*x1 + 4*x2 + error
X = np.stack((x1,x2)).T

From this data generating process it is clear that the correct values are close to 3 and 4. Next I define my class to perform the regression:

class GLM:
def __init__(self):
    self.weights = np.ones(np.shape(X)[1])

def logLikelihoodDerivative(self, X, y):
    derivative = np.dot(y - np.exp(np.dot(X, self.weights)), X)

    print(self.weights) # to observe how weights are updated
    return derivative

def fit(self, X, y, numberOfIterations=10000):
    for _ in range(0, numberOfIterations):
        self.gradientAscent(X, y)

def gradientAscent(self, X, y, learningRate=0.0001):
    gradient = self.logLikelihoodDerivative(X, y)
    self.weights = self.weights + learningRate * gradient
    return self.weights

Now I apply my model:

test = GLM()
test.fit(X, y)
print("My output:     ", test.weights[0], " and ", test.weights[1])

Finally I try to check my outcome with the Python library statsmodels:

import statsmodels.api as sm
check = sm.GLM(y, X, family=sm.families.Poisson()).fit()
print("Correct output:", check.params[0], " and ", check.params[1])

The results are:

My output:      nan  and  nan
Correct output: 2.9999999999999982  and  4.000000000000005

(Note: I am assuming a Poisson distribution).

In the code above the weights explode and become inf/-inf very quickly. I do not understand why.

What is going wrong in my code?

Answer 1

You are doing multiple things wrong.

First, as noticed in the comments, your data-generating process is incorrect. If you want to do Poisson regression, the conditional distribution of your data should be Poisson, not Gaussian. The code should be something like

x1 = np.random.normal(10, 5, 1000)
x2 = np.random.normal(20, 5, 1000)
y = np.random.poisson(np.exp(3*x1 + 4*x2), 1000)
X = np.stack((x1,x2)).T

But this would throw

ValueError: lam value too large

Because np.exp(3*x1 + 4*x2) returns huge numbers.

The second problem is your choice of the optimization algorithm. Gradient descent is a rather poor algorithm. Both R and Python's statmodels use the iteratively reweighted least squares algorithm. There are also other alternatives like the Newton-Raphson algorithm. You could also use a general-purpose optimizer, but if you use a different algorithm than implemented in the software you are comparing to, you should not expect to see the same results.

Finally, gradient descent is sensitive to the choice of the learning rate. In the code you showed us, you are using a hard-coded learning rate. With different learning rates, you could see better results.