# How to do Maximum Likelihood Estimation (MLE) of a Poisson Regression using numpy

I am currently trying to learn how MLE in a poisson regression context works. As such I am trying to compute a poisson regression from scratch using numpy. Furthermore, I try to solve the MLE using gradient descent. However, the loss I am computing stays constant and my guess is that my representation (matrix representation) for the gradient is not correct.

Following an example I found here: https://web.stanford.edu/class/archive/stats/stats200/stats200.1172/Lecture27.pdf

Here it says the gradient should be represented like so: This is how I implemented it using numpy where X is the matrix of covariates and weights are the coefficients I am trying to learn:

y_pred = np.exp(np.dot(X, weights))


This is my full code:

n_samples, n_features = X.shape
weights = np.zeros(n_features)

def forward(X, weights):
return np.exp(np.dot(X, weights))

lr = 0.01
n_iter = 3000
for i in range(n_iter):
# predicting
y_pred = forward(X, weights)
# computing loss
loss = np.sqrt(np.mean((y-y_pred)**2))
if i % 250 == 0:
print(loss)

# calculating gradient and updating weights

• Having glanced at the source you have supplied, this looks to me like a statistical inference problem (i.e. estimating $\beta$ using maximum likelihood), using a numerical method (i.e. Newton-Raphson). On that basis, I am unable to understand why your code necessitates specifying a loss function; surely you only need to be monitoring the log-likelihood for convergence? Furthermore, I am confused by why you have chosen to endow this statistical inference problem with semantics from neural networks e.g. forward. Please may you clarify? – microhaus Jun 10 at 17:35
• @microhaus thanks for reaching out! The reason I use semantics from neural networks in my code is only because I was using pytorch previously which has a PoissonNLLLoss definition. Now I just transitioned my code to a more "from scratch" version using numpy, as I am trying to learn how this works under the hood. I think you are correct that I should be monitoring the log-likelihood instead of the RMSE, but so far I was not able to get this working :( In the example I provided they use Newton-Raphson but I thought it must also be possible to solve this using gradient descent or am I wrong? – Folo Molo Jun 10 at 18:33
• That you are adapting PyTorch clarifies things significantly. Here are some thoughts. 1. If developing understanding of how the statistical inference and numerical method works is your priority, then code it using Numpy. 2. If this Poisson regression wiki is what you have in mind, then yes, gradient descent and Newton-Raphson will work. 3. Depending on whether you wish to vectorise your code to do multivariate updating and the scope of your problem, Newton-Raphon might be computationally demanding, due to inversion of a Hessian. [...] – microhaus Jun 10 at 19:00