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: MLE gradient

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))
gradient = np.dot(X.T,(y - y_pred))

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))

def gradient(y, y_pred, X):
    gradient = np.dot(X.T,(y - y_pred))
    return gradient

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: 
    # calculating gradient and updating weights
    calc_gradient = gradient(y, y_pred, X)
    weights -= lr * calc_gradient

Any help on what I am doing wrong here is highly appreciated!

  • $\begingroup$ 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? $\endgroup$ – microhaus Jun 10 at 17:35
  • $\begingroup$ @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? $\endgroup$ – Folo Molo Jun 10 at 18:33
  • $\begingroup$ 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. [...] $\endgroup$ – microhaus Jun 10 at 19:00
  • $\begingroup$ Coding algorithms from scratch in order to better understand them will definitely enhance your understanding. It might be worthwhile, before you do any coding, to sit down and derive the algorithm for parameter estimation, be it gradient descent or Newton-Raphson, for your particular problem case. I have found that this has never led me astray, in that it forces one to unambiguously and precisely define what will be implemented. Insufficient clarity at the pre-algorithmic stage will be revealed and more difficult to resolve when you additionally have to worry about coding bugs. [...] $\endgroup$ – microhaus Jun 10 at 19:09
  • $\begingroup$ And in my opinion, the process of working through the derivation allows one to naturally derive a checking facility to see if one has implemented their algorithm correctly. E.g. Monitoring log-likelihood for convergence in the case of maximum likelihood with gradient descent. There are other checks you can do if you have gradient expressions e,g. finite differences. If you are struggling with the derivation, consider ask another question. $\endgroup$ – microhaus Jun 10 at 19:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.