I'm looking for someone to give me a hand with a calculus / max. likelihood problem & Python. I can't seem to get a MLE to converge, and I suspect it's because my gradient calculations are wrong (or at least, I want to rule that out).

I'm using scipy's minimize and BFGS to fit a lognormal survival function. Here are some notes on the lognormal distribution you may find useful: 1 and 2

The log likelihood for a survival function is

$$ ll(\theta) = \sum_{\text{only observed events}} \log{h(\theta, t_i)} + \sum_{\text{all events}} \log{S(\theta, t_i)} $$

where $h(\theta, t)$ is the hazard ratio, and $S(\theta, t)$ is the survival function.

Here's my implementation in my Python library: 1. This code is what I want to confirm is correct.

Okay so what do I know: For testing, I generate some fake log-normal($\mu, \sigma$) data and censor some of it (script below).

  1. I get a valid convergence only when ~0.06 < $\sigma$ < ~4. All other values will fail convergence with 'Desired error not necessarily achieved due to precision loss.' $\mu$ does not have this behaviour.

  2. When I remove the jac part from scipy's minimize, convergence seems to work. So this leads to suspect it's my computed gradients.

  3. When all values are observed (no censorship), the model always converges.

    from lifelines import LogNormalFitter
    import numpy as np

    N = 160000
    mu = 3 * np.random.randn()
    sigma = 3 # try 0.03, and 12.0 to see it fail

    X, C = np.exp(sigma * np.random.randn(N) + mu), np.exp(np.random.randn(N) + mu)
    E = X <= C
    T = np.minimum(X, C)

    lnf = LogNormalFitter().fit(T, E)

    assert abs(mu - lnf.mu_) < 0.01
    assert abs(sigma - lnf.sigma_) < 0.01
  • $\begingroup$ Just a little comment before I check out this in detail but if you are trying to rule out gradient computations, you should numerically compute the gradient of your objective function using a Finite Difference technique and make sure your gradient computations are close to the Finite Difference result. $\endgroup$ – spektr Jan 31 at 20:37
  • $\begingroup$ Good suggestion. Scipy has a feature to do this, check_grad, and thing looks okay there - at least, the value outputted is very small. Ex: check_grad(_negative_log_likelihood, gradient_function, [0, 0], log(T), E) $\endgroup$ – Cam.Davidson.Pilon Jan 31 at 20:40
  • $\begingroup$ Does the removal of jac make the optimizer tool in scipy default to a Finite Difference approximation for the gradient? I checked out the documentation on the function and it appears it might. If so, it does seem your gradient computation might be off. Can you try doing the gradient check with random values for your parameters and tell me what you get? $\endgroup$ – spektr Jan 31 at 20:51
  • 1
    $\begingroup$ Ah good hunch: gist.github.com/CamDavidsonPilon/… some clear problems there. $\endgroup$ – Cam.Davidson.Pilon Jan 31 at 20:58
  • 1
    $\begingroup$ Just doing some more testing now. I've found that the sigma gradient is pretty far it's finite difference calculation. In doing some more simplifying, I think it's my modelling for sigma = exp(v) that may be causing problems. $\endgroup$ – Cam.Davidson.Pilon Jan 31 at 22:10

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.