Comment puis-je transmettre des arguments à d'autres fonctions (généralement et via scipy) ?

Question

Je essaie de minimiser une fonction qui produit du chi-square via scipy et trouver les mu, sigma, normc qui fournissent le meilleur ajustement pour une superposition gaussienne.

from math import exp
from math import pi
from scipy.integrate import quad
from scipy.optimize import minimize
from scipy.stats import chisquare
import numpy as np

# devinez les valeurs initiales pour le chi-square minimisé
mu, sigma = np.mean(mydata), np.std(mydata) # mes points de données sont dans mydata
normc = 1/(sigma * (2*pi)**(1/2)) 

gauss = lambda x: normc * exp( (-1) * (x - mu)**2 / ( 2 * (sigma **2) ) ) # Distribution gaussienne

# supposons que j'ai des limites de classe pré-définies en tant que liste appelée binbound

def expvalperbin(binbound,mu,sigma,normc):
    # calcule la valeur d'attente par classe
    ans = []
    for index in range(len(binbound)):
        if index != len(binbound)-1:
            ans.append( quad( gauss, binbound[index], binbound[index+1])[0] )
    return ans

expvalguess = expvalperbin(binbound,mu,sigma,normc)
obsval = countperbin(binbound,mydata)
arglist = [mu,sigma,normc]

def chisquareopt(obslist,explist):
    return chisquare(obslist,explist)[0]

chisquareguess = chisquareopt((obsval,expvalguess), expvalguess, args=arglist)

result = minimize( chisquareopt(obsval,expvalguess), chisquareguess   )
print(result)

En exécutant ce code, cela me donne cette erreur :

TypeError: chisquareopt() got an unexpected keyword argument 'args'

J'ai quelques questions :

1) Comment puis-je écrire une fonction pour permettre aux arguments d'être passés à ma fonction chisquareopt ?

2) Comment puis-je dire si scipy optimisera les paramètres [mu, sigma, normc] qui donnent le chi-square minimum ? Comment pourrais-je trouver ces paramètres à partir de l'optimisation ?

3) Il est difficile de savoir si je progresse ici ou non. Suis-je sur la bonne voie ?

EDIT : Si c'est pertinent, j'ai une fonction qui prend en entrée [mu, sigma, normc] et renvoie une liste de sous-listes, chaque sous-liste contenant une combinaison possible de [mu, sigma, normc] (où la liste externe couvre toutes les combinaisons possibles de paramètres dans des plages spécifiées).

Answer 1

Généralement, ces fonctions scipy transmettent le tuple args de valeurs à votre code inchangé. Je devrais vérifier le code, mais avec

minimize(myfunc, x0, args=(y,z))

def myfunc(x, y, z): 

minimize prend la valeur actuelle de la variable x (scalaire ou tableau, selon ce que x0 ressemble), et le paramètre args, et construit

args = tuple(x) + args
myfunc(*args)

En d'autres termes, il joint le tuple args avec la variable d'itération et le transmet à votre fonction. Ainsi, toute définition intermédiaire de fonction doit fonctionner avec ce motif.

Pour illustrer, définissez une fonction qui prend un tuple args générique.

In [665]: from scipy.optimize import minimize
In [666]: def myfunc(*args):
     ...:     print(args)
     ...:     return np.abs(args[0])**2
     ...: 
In [667]: myfunc(1,2,3)
(1, 2, 3)
Out[667]: 1
In [668]: myfunc(2,2,3)
(2, 2, 3)
Out[668]: 4
In [669]: minimize(myfunc, 10, args=(2,3))
(array([ 10.]), 2, 3)
(array([ 10.00000001]), 2, 3)
(array([ 10.]), 2, 3)
(array([ 8.99]), 2, 3)
....
(array([-0.00000003]), 2, 3)
Out[669]: 
      fun: 1.7161984122524196e-15
 hess_inv: array([[ 0.50000001]])
      jac: array([-0.00000007])
  message: 'Optimization terminated successfully.'
     nfev: 15
      nit: 4
     njev: 5
   status: 0
  success: True
        x: array([-0.00000004])

---

(discussion supprimée sur la confusion concernant les paramètres à minimiser. Voir l'autre réponse ou mon historique de modifications)

Answer 2

J'ai simplifié quelque peu votre problème pour vous donner une idée de votre question 2).

Tout d'abord, j'ai défini de manière statique votre histogramme obslist et le nombre de points de données N en tant que variables globales (ce qui simplifie un peu les signatures de fonction). Ensuite, j'ai défini de manière statique les limites des bacs dans expvalperbin, en supposant 9 bacs avec une largeur fixe de 5 et en commençant par le premier bac à 30 (donc l'histogramme s'étend de 30 à 75).

Ensuite, j'utilise optimize.fmin (Nelder-Mead) au lieu de optimize.minimize. La raison d'utiliser fmin au lieu de minimize est que le passage de paramètres supplémentaires via args=(x,y) ne semble pas fonctionner dans le sens où les paramètres supplémentaires restent aux valeurs fixes depuis la toute première invocation. Ce n'est pas ce que vous voulez : vous voulez optimiser mu et sigma simultanément.

Grâce à ces simplifications, nous avons le script suivant (certainement très peu pythonique) :

from math import exp
from math import pi
from scipy.integrate import quad
from scipy.optimize import fmin
from scipy.stats import chisquare

obslist = [12, 51, 144, 268, 264, 166, 75, 18, 2] # histogramme, 1000 observations
N = 1000 # nombre de points de données

def gauss(x, mu, sigma):
    return 1/(sigma * (2*pi)**(1/2)) * exp( (-1) * (x - mu)**2 / ( 2 * (sigma **2) ) )

def expvalperbin(mu, sigma):
    e = []
    # limites de bac codées en dur
    for i in range(30, 75, 5):
        e.append(quad(gauss, i, i + 5, args=(mu, sigma))[0] * N)
    return e

def chisquareopt(args):
    # args[0] = mu
    # args[1] = sigma
    return chisquare(obslist, expvalperbin(args[0], args[1]))[0]

# hypothèses initiales
initial_mu = 35.5
initial_sigma = 14

result = fmin(chisquareopt, [initial_mu, initial_sigma])

print(result)

Optimisation terminée avec succès.

Valeur actuelle de la fonction : 2.010966

Itérations : 49

Évaluations de fonction : 95

[ 50.57590239 7.01857529]

En passant, l'histogramme obslist est un échantillon aléatoire de 1000 points d'une distribution normale N(50.5, 7.0). N'oubliez pas que ce sont mes toutes premières lignes de code Python, alors s'il vous plaît ne me jugez pas sur le style. Je voulais juste vous donner une idée de la structure générale du problème.