Quelles exigences spécifiques la fonction transmise à scipy.optimize.curve_fit doit-elle remplir pour fonctionner?

Question

Je travaille sur l'ajustement des modèles statistiques aux distributions en utilisant la fonction hist de matplotlib. Par exemple mon code correspond à une distribution exponentielle en utilisant le code suivant:

try: 

     def expDist(x, a, x0): 
      return a*(exp(-(x/x0))/x0) 

     self.n, self.bins, patches = plt.hist(self.getDataSet(), self.getDatasetSize()/10, normed=1, facecolor='blue', alpha = 0.55) 
     popt,pcov = curve_fit(expDist,self.bins[:-1], self.n, p0=[1,mean]) 
     print "Fitted gaussian curve to data with params a %f, x0 %f" % (popt[0], popt[1]) 
     self.a = popt[0] 
     self.x0 = popt[1] 

     self.fitted = True 
    except RuntimeError: 
     print "Unable to fit data to exponential curve" 

Ce qui fonctionne très bien, mais quand je le modifier pour faire la même chose pour une distribution uniforme entre a & b,

def uniDist(x, a, b): 
     if((x >= a)and(x <= b)): 
      return float(1.0/float(b-a)) 
     else: 
      return 0.000 

    try: 



     self.n, self.bins, patches = plt.hist(self.getDataSet(), self.getDatasetSize()/10, normed=1, facecolor='blue', alpha = 0.55) 
     popt,pcov = curve_fit(uniDist,self.bins[:-1], self.n, p0=[a, b]) 
     print "Fitted uniform distribution curve to data with params a %f, b %f" % (popt[0], popt[1]) 
     self.a = popt[0] 
     self.b = popt[1] 

     self.fitted = True 
    except RuntimeError: 
     print "Unable to fit data to uniform distribution pdf curve" 

les accidents de code, avec

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

La question semble être que quelque part dans curve_fit, e La fonction essaye d'appeler la fonction à équiper (expDist, et uniDist dans ce cas) avec un ensemble de valeurs itérable, mais je n'arrive pas à comprendre comment la fonction expDist est capable de prendre n'importe quoi itérative sans se bloquer?

Answer 1

Votre soupçon est en partie correct. curve_fit passe en effet une itérable à la fonction, mais pas n'importe laquelle itérable: numpy.ndarray. Ceux-ci arrivent à posséder des opérateurs arithmétiques vectorisés, donc

a*(exp(-(x/x0))/x0) 

les tableaux d'entrée fonctionnera tout simplement élément par élément sur sans erreur (et avec une sortie correcte). Il n'y a même pas beaucoup de magie impliquée: pour chaque évaluation de la fonction, les paramètres a et x0 seront des scalaires, seul x est un tableau.

Maintenant, le problème avec uniDist est qu'il ne contient pas seulement des opérateurs arithmétiques: il contient également des opérateurs de comparaison. Ces travaux très bien tant que seule une seule matrice est comparée à un scalaire:

>>> import numpy as np 
>>> a = np.arange(5) 
>>> a 
array([0, 1, 2, 3, 4]) 
>>> a>2 
array([False, False, False, True, True], dtype=bool) 

Ce qui précède démontre que l'utilisation des opérateurs de comparaison sur un tableau et un scalaire à nouveau des résultats élément par élément. L'erreur que vous voyez se pose lorsque vous essayez d'appliquer un opérateur logique à deux de ces tableaux booléens:

>>> a>2 
array([False, False, False, True, True], dtype=bool) 
>>> a<4 
array([ True, True, True, True, False], dtype=bool) 
>>> (a>2) and (a<4) 
Traceback (most recent call last): 
    File "<stdin>", line 1, in <module> 
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() 

Le message d'erreur est un peu déroutant. On peut remonter au fait que python essayait de trouver un résultat unique pour array1 and array2 (qui en python natif retournerait soit un tableau basé sur leur vide). Cependant, on soupçonne que ce n'est pas ce que vous voulez faire, et résiste à la tentation de deviner.

Puisque vous voulez que votre fonction fonctionne par éléments sur les deux tableaux booléens (qui proviennent de l'opération de comparaison), vous devrez utiliser l'opérateur &. Ceci est "binaire et" en python natif, mais pour les tableaux numpy cela vous donne le "logique et" des tableaux. Vous pouvez également utiliser numpy.logical_and (ou dans votre cas scipy.logical_and) d'être plus explicite:

>>> (a>2) & (a<4) 
array([False, False, False, True, False], dtype=bool) 
>>> np.logical_and(a>2,a<4) 
array([False, False, False, True, False], dtype=bool) 

Notez que pour le & cas, vous devez toujours parenthésée vos comparaisons, car encore une fois a>2&a<4 serait ambiguë (pour le programmeur) et faux (compte tenu de ce que vous veut faire).Puisque le "binary et" des booléens se comporteront exactement comme vous l'attendez, il est prudent de réécrire votre fonction pour utiliser & au lieu de and pour comparer les deux comparaisons. Cependant, il reste encore une étape à modifier: dans le cas des entrées ndarray, le if se comportera également différemment. Python ne peut s'empêcher de faire un seul choix dans un if, ce qui est également vrai si vous y mettez un tableau. Mais ce que vous voulez réellement faire est de contraindre les éléments de votre sortie par élément (à nouveau). Vous devez donc soit boucler votre tableau (ne le faites pas), soit refaire ce choix de façon vectorisée. Ce dernier est idiomatiques en utilisant numpy/scipy:

import scipy as sp 
def uniDist(x, a, b): 
    return sp.where((a<=x) & (x<=b), 1.0/(b-a), 0.0) 

Ce (à savoir numpy.where) retourne un tableau de la même taille que x. Pour les éléments dont la condition est True, la valeur de la sortie sera 1/(b-a). Pour le reste, la sortie est 0. Pour scalaire x, la valeur de retour est un scalaire numpy. Notez que j'ai supprimé la conversion float dans l'exemple ci-dessus, puisque avoir 1.0 dans le numérateur vous donnera certainement la vraie division, malgré votre utilisation de python 2. Bien que je suggère d'utiliser python 3, ou au moins from __future__ import division.

---

Note mineure: même pour un cas scalaire, je suggère en utilisant de l'opérateur de python enchaînant pour la comparaison, qui se prête à cet effet. Ce que je veux dire, c'est que vous pouvez simplement faire if a <= x <= b: ..., et contrairement à la plupart des langues, cela sera fonctionnellement équivalent à ce que vous avez écrit (mais plus joli).