scipy n'est pas de l'optimisation et de retours “Souhaité erreur n'est pas nécessairement dû à la perte de précision due”

J'ai le code suivant qui tente de minimiser un journal d'une fonction de vraisemblance.

#!/usr/bin/python
import math
import random
import numpy as np
from scipy.optimize import minimize

def loglikelihood(params, data):
    (mu, alpha, beta) = params
    tlist = np.array(data)
    r = np.zeros(len(tlist))
    for i in xrange(1,len(tlist)):
        r[i] = math.exp(-beta*(tlist[i]-tlist[i-1]))*(1+r[i-1])
    loglik  = -tlist[-1]*mu
    loglik = loglik+alpha/beta*sum(np.exp(-beta*(tlist[-1]-tlist))-1)
    loglik = loglik+np.sum(np.log(mu+alpha*r))
    return -loglik

atimes = [ 148.98894201,  149.70253172,  151.13717804,  160.35968355,
        160.98322609,  161.21331798,  163.60755544,  163.68994973,
        164.26131871,  228.79436067]
a= 0.01
alpha = 0.5
beta = 0.6
print loglikelihood((a, alpha, beta), atimes)

res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'BFGS',args = (atimes,))
print res

Il me donne 

28.3136498357
./test.py:17: RuntimeWarning: invalid value encountered in log
  loglik = loglik+np.sum(np.log(mu+alpha*r))
   status: 2
  success: False
     njev: 14
     nfev: 72
 hess_inv: array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1]])
      fun: 32.131359359964378
        x: array([ 0.01,  0.1 ,  0.1 ])
  message: 'Desired error not necessarily achieved due to precision loss.'
      jac: array([ -2.8051672 ,  13.06962156, -48.97879982])

Avis qu'il n'a pas réussi à optimiser les paramètres à tous et réduit la valeur 32 est plus grand que 28, qui est ce que vous obtenez avec un= 0,01, alpha = 0.5, bêta = 0.6 . Il est possible que ce problème pourrait être évité en choisissant mieux conjectures initiales, mais si oui, comment puis-je faire cela automatiquement?

  • Je pense que vous voulez maximiser LL, pas le minimiser. Si vous êtes à la réduction d'une somme de carrés, l'optimisation de LL.
  • Oui. Avis de la fonction renvoie -loglik qui gère cela.
  • Juste une remarque - une fois, j'ai eu un problème qui a partagé les mêmes symptômes que les vôtres, mais la cause était totalement différent. Il s'est avéré que j'avais un bug dans mon gradient de la fonction, donc quand je l'ai adoptée dans le la routine via le jac paramètre, la routine ne pouvait pas travailler. Les erreurs ont été cryptique et c'est seulement sur la ré-inspection de mon code que j'ai identifié le bug. Cela dit, la réponse ci-dessous qui utilise Nelder-Mead vraiment aidé parce qu'il pourrait optimiser sans le dégradé, et a donné la bonne réponse pour moi, de m'aider à me rendre compte que le problème était avec le bug dans mon gradient de la fonction.
InformationsquelleAutor felix | 2014-07-15
  1. 14

    J'ai copié votre exemple et essayé un peu. Ressemble si vous vous en tenez à BFGS solveur, après quelques itération, le mu+ alpha * r aura un certain nombres négatifs, et c'est ainsi que vous obtenez la RuntimeWarning.

    Le plus simple correctif, je pense, de passer à Nelder Mead solveur.

    res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'Nelder-Mead',args = (atimes,))

    Et il vous donnera ce résultat:

    28.3136498357
  status: 0
    nfev: 159
 success: True
     fun: 27.982451280648817
       x: array([ 0.01410906,  0.68346023,  0.90837568])
 message: 'Optimization terminated successfully.'
     nit: 92
    InformationsquelleAutor JimmyK
  2. 3

    Regarder pour les valeurs négatives de la fonction log (), à les résoudre et de dire à l'optimiseur qu'ils sont mauvais, par l'ajout d'une pénalité:

    #!/usr/bin/python
import math
import random
import numpy as np
from scipy.optimize import minimize

def loglikelihood(params, data):
    (mu, alpha, beta) = params
    tlist = np.array(data)
    r = np.zeros(len(tlist))
    for i in xrange(1,len(tlist)):
        r[i] = math.exp(-beta*(tlist[i]-tlist[i-1]))*(1+r[i-1])
    loglik = -tlist[-1]*mu
    loglik += alpha/beta*sum(np.exp(-beta*(tlist[-1]-tlist))-1)
    argument = mu + alpha * r
    limit = 1e-6
    if np.min(argument) < limit:
        # add a penalty for too small argument of log
        loglik += np.sum(np.minimum(0.0, argument - limit)) / limit
        # keep argument of log above the limit
        argument = np.maximum(argument, limit)
    loglik += np.sum(np.log(argument))
    return -loglik

atimes = [ 148.98894201,  149.70253172,  151.13717804,  160.35968355,
        160.98322609,  161.21331798,  163.60755544,  163.68994973,
        164.26131871,  228.79436067]
a= 0.01
alpha = 0.5
beta = 0.6
print loglikelihood((a, alpha, beta), atimes)

res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'BFGS',args = (atimes,))
print res
    InformationsquelleAutor optimizer
  3. 2

    Face à la même avertissement, je l'ai résolu par une réécriture de la log-vraisemblance de la fonction pour obtenir log(params) et log(data) comme arguments, au lieu de paramètres et de données.

    Donc, je ne peux pas utiliser np.log() dans la fonction de vraisemblance ou de la Jacobienne, si possible.

    InformationsquelleAutor Sahar
  4. 1

    Une autre solution (qui a fonctionné pour moi) est à l'échelle de votre fonction (et dégradés) à des valeurs proches de 0. Par exemple, mon problème est venu quand j'ai eu à évaluer la log-vraisemblance de 60k points. Cela voulait dire que mes log-vraisemblance a été d'un très grand nombre. Sur le plan conceptuel, la log-vraisemblance est une très très spikey fonction.

    Les gradients commencé large (à monter cette spikey de montagne), puis est devenu modérément petits, mais jamais inférieure à la valeur par défaut gtol paramètre dans la BGFS routine (qui est le seuil que tous les gradients doivent être en dessous de résiliation). Aussi, à cette époque, j'avais essentiellement arrivé au les valeurs correctes (j'ai été en utilisant les données générées donc je savais que les vraies valeurs).

    Ce qui se passait, c'était que mon gradients ont été environ. 60k * average individual gradient value, et même si le average individual gradient value était petite, disons moins de 1e-8, 60 k * 1e-8 > gtol. Donc, je n'ai jamais été satisfaisant, le seuil, même si j'étais arrivé à la solution.

    Sur le plan conceptuel, à cause de cette très spikey de montagne, l'algorithme a été de faire de petites étapes, mais l'intensification sur le vrai minimum et n'a jamais atteint average individual gradient << 1e-8 qui implique mon gradients ne suis jamais allé en vertu de gtol.

    Deux solutions:

    1) à l'Échelle de votre log-vraisemblance et dégradés par un facteur, comme 1/n où n est le nombre d'échantillons.

    2) à l'Échelle de votre gtol: par exemple "gtol": 1e-7 * n

    InformationsquelleAutor Cam.Davidson.Pilon