Dans Scipy comment et pourquoi curve_fit calculer la covariance des estimations des paramètres

J'ai été en utilisant scipy.optimize.leastsq pour s'adapter à certaines données. J'aimerais obtenir des intervalles de confiance de ces estimations, donc, je regarde dans le cov_x de sortie, mais la documentation ne précise pas ce que c'est et comment obtenir la matrice de covariance pour mes paramètres de cette.

Tout d'abord, il dit que c'est un Jacobien, mais dans le notes il dit aussi que "cov_x est un Jacobien approximation de la Hessienne" de sorte qu'il n'est pas réellement un Jacobien mais un de Hesse en utilisant une approximation de la Jacobienne. Laquelle de ces affirmations est correcte?

Deuxièmement, cette phrase m'est déroutant:

Cette matrice doit être multiplié par la variance résiduelle pour obtenir la covariance des estimations des paramètres – voir curve_fit.

J'ai en effet aller regarder le code source pour curve_fit où ils n':

s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq

qui correspond à la multiplication des cov_x par s_sq mais je ne trouve pas cette équation dans toute référence. Quelqu'un peut m'expliquer pourquoi cette équation est correcte?
Mon intuition me dit qu'il devrait être dans l'autre sens depuis cov_x est censé être un dérivé (Jacobien ou de Hesse) alors que je pensais:
cov_x * covariance(parameters) = sum of errors(residuals) où sigma(parameters) est la chose que je veux.

Comment puis-je connecter la chose curve_fit est en train de faire avec ce que je vois à, par exemple. wikipedia:
http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations

InformationsquelleAutor HansHarhoff | 2013-02-13
  1. 26

    OK, je crois que j'ai trouvé la réponse. D'abord la solution:
    cov_x*s_sq est tout simplement la covariance des paramètres qui est ce que vous voulez. En prenant la racine carrée de la diagonale éléments vont vous donner l'écart-type (mais attention à la covariances!).

    Variance résiduelle = réduction de chi carré = s_sq = somme[(f(x)-y)^2]/(N-n), où N est le nombre de points de données et n est le nombre de paramètres d'appareillage. Réduction de chi carré.

    La raison de ma confusion, c'est que cov_x donnée par leastsq n'est pas réellement ce qu'on appelle cov(x) dans d'autres endroits, c'est plutôt la réduction des cov(x) ou de fractions de cov(x). La raison pour laquelle il n'apparaît pas dans les autres références, c'est que c'est un simple changement d'échelle qui est utile dans les calculs numériques, mais n'est pas pertinent pour un manuel.

    À propos de Hesse contre Jacobien, la documentation est mal formulé. C'est la Hesse, qui est calculé dans les deux cas, comme il est évident depuis le Jacobien est nul au minimum. Ce qu'ils veulent dire, c'est qu'ils sont en utilisant une approximation de la Jacobienne pour trouver la toile de jute.

    Une note supplémentaire. Il semble que la curve_fit résultat ne tient pas réellement compte de la taille absolue de l'erreur, mais seulement de tenir compte de la taille relative de la sigmas fourni. Cela signifie que la pcov retourné ne change pas même si le errorbars changer par un facteur d'un million. Ce n'est évidemment pas la droite, mais semble être le standard de la pratique de l'ie. Matlab fait la même chose lors de l'utilisation de leur ajustement de la Courbe boîte à outils. La procédure est décrite ici: https://en.wikipedia.org/wiki/Linear_least_squares_(mathématiques)#Parameter_errors_and_correlation

    Il semble assez simple à faire une fois l'optimum a été trouvé, au moins pour le Linéaire des moindres carrés.

    • La formule pour le chi-carré devrait être multiplié par 1.0 / erreurs^2, je pense?
    • “Il semble que la curve_fit résultat ne tient pas réellement compte de la taille absolue de l'erreur, mais seulement de tenir compte de la taille relative de la sigmas fourni.” Il y a un drapeau pour cette: absolute_sigma. Si elle est désactivée (par défaut), puis curve_fit, l'estimation de la var(y) basée sur vos données; sinon, il faudra utiliser votre condition sigma valeurs.
    InformationsquelleAutor HansHarhoff
  2. 7

    J'ai trouvé cette solution au cours de ma recherche pour une question similaire, et j'ai seulement une petite amélioration sur HansHarhoff de réponse. La sortie complète de leastsq fournit une valeur de retour infodict, qui contient infodict['fvec'] = f(x) -y. Ainsi, pour calculer la réduction de chi carré = (dans les notations ci-dessus)

    s_sq = (infodict['fvec']**2).sum()/(N-n)

    BTW. Grâce HansHarhoff pour faire la plupart de levage lourds de résoudre ce problème.

    • Très sympa! Exactement ce que je cherchais. Specificially: result=scipy.optimize.leastsq(...,, full_output=True);s_sq = (result]['fvec']**2).sum()/(len(result[2]['fvec'])-len(result[0]))
    • quelle version avez-vous à l'appui de la full_output?
    • Je ne sais pas ce que la version minimum serait d'appuyer le full_output option, mais j'utilise scipy 0.13.3 et 0.14.1.
    • il vous manque un [2 là ... s_sq = (result[2]['fvec']**2).sum()/(len(result[2]['fvec'])-len(result[0]))
    InformationsquelleAutor Jim Parker