Le monde regorge de phénomènes variés. Les comprendre, c'est créer un modèle qui explique leur comportement. Si le modèle peut prédire les données d'observation d'un certain phénomène, on peut dire que le modèle comprend un certain phénomène.

Il existe deux manières principales d'ajuster le modèle aux données.

** méthode des moindres carrés **
** méthode du maximum de vraisemblance **

En fait, lequel doit être utilisé quand? Cette fois, je vais essayer ** l'ajustement de courbe ** à travers un exemple simple.

Tâche

Prenons l'exemple d'une étude de la conscience classique. Faites asseoir le sujet devant le moniteur et présentez le stimulus visuel sur le moniteur. Les stimuli visuels comprennent les signaux et le bruit, et les sujets sont invités à indiquer s'ils ont vu le signal pour chaque essai.

Intuitivement, plus le rapport de signal est élevé par rapport au bruit, plus il est facile à détecter. En fait, si la force du signal est sur l'axe horizontal et que les performances de détection du sujet sont sur l'axe vertical, par exemple, ce sera comme suit.

C'est ce qu'on appelle la fonction psychométrique, mais cela n'a pas d'importance maintenant, et vous pouvez clairement voir la relation ** non linéaire **. Créons un modèle qui représente cette relation non linéaire et insérons-le dans les données. Si vous pouvez faire cela, vous pouvez prédire les performances de détection à ce moment même si la force du signal n'est pas utilisée dans l'expérience. Si vous pouvez faire cela, vous comprenez ce phénomène non linéaire.

Les données utilisées cette fois sont résumées ci-dessous.

signal	performance	number of trials
0	0.5	50
0.1	0.73	45
0.2	0.84	40
0.4	0.92	35
0.8	0.98	30

Les informations requises pour l'ajustement sont les performances correspondant à la force de chaque signal et le nombre d'essais des stimuli lorsque les sujets ont été invités à effectuer les 200 essais.

Cette fois, nous utiliserons python, alors importez les bibliothèques requises et mettez les données dans la liste.

`curvefitting.py`



# import libraries
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import optimize

# experimental data ========
# signal strength
signal = [0, 0.1, 0.2, 0.4, 0.8]

# accuracy
accuracy = [0.5, 0.73, 0.84, 0.92, 0.98]

# number of trials
ntr = [50, 45, 40, 35, 30]

modèle

Il existe de nombreux modèles non linéaires, mais nous avons constaté que ** Weibull cumulatif ** fonctionne bien pour des tâches de détection telles que celle-ci.

Lorsque la force d'un signal $ s $ est présentée comme un stimulus, la probabilité que le sujet réponde correctement $ p $ est exprimée comme suit.


p = 0.5 + (0.5 - \lambda)f((\frac{s}{\theta})^\beta)

$ f $ est $ f (x) = (1 --e ^ {-x}) $. Il existe trois paramètres gratuits qui déterminent la forme de l'effacement cumulatif, $ \ lambda $, $ \ theta $, $ \ beta $. Puisque nous connaissons la force du signal $ s $ et les performances du sujet $ p $ dans l'expérience, nous les utiliserons pour trouver des ** paramètres gratuits qui correspondent aux données **.

Dans la programmation réelle, tous les paramètres libres sont placés dans une variable de liste appelée $ x $. $ x [0] = \ lambda $, $ x [1] = \ theta $, $ x [2] = \ beta $ et ainsi de suite.

`curvefitting.py`



# model for performance ========
def probabilityCorrect(x, signal):
    # x: model parameters
    # s: signal strength
    
    # probability of correct 
    return 0.5 + (0.5 - x[0])*(1 - np.exp(-(signal/x[1])**x[2]))

Ajustement du modèle (méthode du carré minimum)

Tout d'abord, essayez la ** méthode des moindres carrés **. La méthode des moindres carrés est une méthode de recherche de paramètres libres tels que ** "La différence entre la valeur prédite du modèle et les données réelles est minimisée" **. Le carré est utilisé simplement parce que vous n'avez pas à penser au signe.

Pour trouver les paramètres gratuits, créez quelque chose appelé ** fonction de coût **. Le flux consiste à le lancer dans l'algorithme et à ce que l'algorithme trouve un paramètre libre qui minimise sa fonction de perte.

Avec la méthode des moindres carrés, la fonction de perte est très intuitive. La valeur prévue du modèle et la différence carrée des données réelles doivent être calculées en fonction de la force de chaque signal, et additionnées pour donner une valeur numérique.

Voici un exemple de mise en œuvre.

`curvefitting.py`



def cost_LSE(x, signal, performance, ntr):       
    # compute cost
    c = 0
    for n in np.arange(len(ntr)):
        # model output 
        p = probabilityCorrect(x, signal[n])
        
        # cost as the summed squared error
        c += ntr[n]*(performance[n] - p)**2
        
    return c

Maintenant que nous avons une fonction de perte, lancez-la dans l'algorithme pour trouver un paramètre libre qui prend la valeur minimale. En python, scipy inclut un package appelé Optimize, alors utilisez-le.

`curvefitting.py`



# initial guess and boundary
x0 = [0, 0.2, 1]
bound = [(0, None),(0.0001, None),(0,None)]

# least squared error
params_LSE = optimize.minimize(cost_LSE,x0,args=(signal,performance,ntr,),method='l-bfgs-b',\
                        jac=None, bounds=bound, tol=None, callback=None,\
                        options={'disp': None, 'maxls': 20, 'iprint': -1,\
                                 'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 15000,\
                                 'ftol': 2.220446049250313e-09, 'maxcor': 10,\
                                 'maxfun': 15000})

Les arguments peuvent prêter à confusion, mais fondamentalement Site officiel Il s'agit du copier-coller par défaut répertorié dans. Veuillez noter que ** vous devez décider de la valeur initiale du paramètre gratuit $ x0 $ **.

L'algorithme utilisé cette fois est L-BFGS-B, qui est un algorithme polyvalent convivial pour la mémoire, mais tout algorithme répète les calculs numériques (numériquement) pour trouver la valeur minimale, de sorte que la valeur initiale est proche de la réponse. C'est idéal. Si la valeur initiale n'est pas pertinente, le calcul de l'algorithme peut ** ne pas converger ** et l'ajustement souhaité ne sera pas possible. Je n'ai pas d'autre choix que d'essayer diverses choses, mais cette fois, empiriquement, la valeur initiale est "Eh bien, cela se produira".

Un autre, ** vous pouvez prérégler la limite des paramètres gratuits $ \ lambda $, $ \ theta $, $ \ beta $ **. Avec des paramètres de plage appropriés, l'algorithme peut éviter d'errer dans un grand nombre d'océans. Cette fois, tous les paramètres libres sont positifs, et nous ne voulons pas que $ x [1] = \ theta $ soit $ 0 $ car il s'agit du dénominateur du modèle, nous définissons donc la borne comme ci-dessus. C'est réglé.

Maintenant, l'algorithme L-BFGS-B trouvera la valeur du paramètre libre qui minimise la fonction de perte.

Ajustement du modèle (méthode d'estimation la plus probable)

À titre de comparaison, essayez la ** méthode du maximum de vraisemblance **. La méthode d'estimation la plus probable est, en un mot, une méthode de recherche de paramètres libres tels que ** «Maximiser la probabilité (probabilité) que le modèle prédit des données réelles lorsqu'un certain paramètre libre est donné» **. est.

J'ai fait une expérience et j'ai obtenu des données. Les données proviennent de la distribution de probabilité d'un modèle, qui est déterminée par les paramètres libres du modèle. La justification est donc de trouver un paramètre libre qui maximise la probabilité (plausibilité, vraisemblance, vraisemblance ** de la valeur réelle des données sortant du modèle.

Plus précisément, ** la vraisemblance est définie par la probabilité conjointe de $ P (data_i | x) $ donnée au modèle par le paramètre libre $ x $ pour chaque événement d'observation $ data_i $. **. Comme nous l'avons appris au lycée, si chaque événement d'observation se produit indépendamment, la probabilité simultanée est exprimée comme le produit de la probabilité que chaque événement se produise.


likelihood = P(data_0 | x)P(data_1 | x)P(data_2 | x)...P(data_n | x)

Je veux créer une fonction qui maximise cela, mais la seule façon de trouver des paramètres gratuits est de créer une fonction de coût dont la valeur minimale peut être trouvée par un algorithme. Cependant, puisque je veux la valeur maximale cette fois, je transmets la probabilité avec un signe négatif à l'algorithme.

Une autre raison est que la probabilité est le produit de nombreuses probabilités, donc si le nombre d'observations est grand, la valeur se rapprochera progressivement de 0, ce qui rendra difficile pour l'algorithme de trouver la valeur maximale. Par conséquent, nous utilisons l'astuce mathématique de prendre un journal. Comme je l'ai appris au lycée, en prenant un journal, le produit des éléments peut être converti en somme des éléments.

En conséquence, la fonction de perte introduite dans l'algorithme par la méthode d'estimation la plus probable est la suivante, le log étant pris et minimisé.


- log(likelihood) = - (log(P(data_0 | x)) + log(P(data_1 | x)) + log(P(data_2 | x)) ... + log(P(data_n | x))

Voici un exemple de mise en œuvre.

`curvefitting.py`



def cost_MLE(x, signal, performance, ntr):       
    # compute cost
    c = 0
    for n in np.arange(len(ntr)):
        # model output 
        p = probabilityCorrect(x, signal[n])
        
        # cost as the negative likelihood
        if 0 < p <= 1:
            c += -ntr[n]*(performance[n]*np.log(p) + (1 - performance[n])*np.log(1 - p))
        
    return c

La valeur prévue du modèle doit être positive pour la commodité de la journalisation. De plus, puisqu'il s'agit d'une probabilité, elle ne peut pas être supérieure à 1. J'ai mis ces deux conditions dans $ if $. Notez que dans cet exemple, la probabilité que le sujet réponde correctement est modélisée comme $ p $, donc la probabilité que le sujet ne réponde pas correctement est $ 1-p $.

Maintenant que nous avons une fonction de perte, nous la jetons dans l'algorithme.

`curvefitting.py`



# initial guess and boundary
x0 = [0, 0.2, 1]
bound = [(0, None),(0.0001, None),(0,None)]

# maximum likelihood
params_MLE = optimize.minimize(cost_MLE,x0,args=(signal,performance,ntr,),method='l-bfgs-b',\
                        jac=None, bounds=bound, tol=None, callback=None,\
                        options={'disp': None, 'maxls': 20, 'iprint': -1,\
                                 'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 15000,\
                                 'ftol': 2.220446049250313e-09, 'maxcor': 10,\
                                 'maxfun': 15000})

Comparaison des résultats d'ajustement

Qu'en est-il du résultat approprié? Faisons un graphique.

`curvefitting.py`



# visualize
signals = np.linspace(0,0.8,100)
acc_mle = np.zeros(100)
acc_lse = np.zeros(100)

for i in np.arange(100):
    acc_mle[i] = probabilityCorrect(params_MLE.x,signals[i])
    acc_lse[i] = probabilityCorrect(params_LSE.x,signals[i])

fig = plt.figure()
fig.suptitle('psychometric function')

ax = fig.add_subplot(111)

ax.set_xlabel('signal strength')
ax.set_ylabel('performance')    
plt.plot(signal, performance, 'ko')
plt.plot(signals, acc_mle, '-b')
plt.plot(signals, acc_lse, '-r')

Tracez les raccords obtenus par la méthode des moindres carrés en rouge et les raccords obtenus par la méthode d'estimation la plus probable en bleu.

…… C'est presque la même chose w Les deux sont bien ajustés.

En fait, la qualité de l'ajustement peut être quantifiée avec l'indice de ** variation expliquée **. Ce qui suit est officiel.


varianceExplained = 1 - var(data - prediction)/var(data)

Nous quantifions dans quelle mesure la variation de la valeur prédite expliquait la variation des données réelles. Calculons cela.

`curvefitting.py`



# compute variance explained
fitted_acc_mle = np.zeros(len(signal))
fitted_acc_lse = np.zeros(len(signal))
for s in np.arange(len(signal)):
    fitted_acc_mle[s] = probabilityCorrect(params_MLE.x,signal[s])
    fitted_acc_lse[s] = probabilityCorrect(params_LSE.x,signal[s])
    
varexp_mle = 1 - (np.var(fitted_acc_mle - performance)/np.var(performance))
varexp_lse = 1 - (np.var(fitted_acc_lse - performance)/np.var(performance))

print('variance explained (MLE): ' + str(varexp_mle))
print('variance explained (LSE): ' + str(varexp_lse))

À la suite du calcul,

variance explained (MLE): 0.999339315626 variance explained (LSE): 0.999394193462

est devenu. Il semble que la méthode des moindres carrés soit un peu meilleure, mais la différence est négligeable. Vous voudrez peut-être essayer les deux pour le moment et choisir celui qui vous convient le mieux.

Résumé récapitulatif

--Utilisez la méthode des moindres carrés ou la méthode d'estimation la plus probable pour ajuster le modèle.

Le flux de création d'un modèle, de création d'une fonction de perte et de le lancer dans un algorithme qui trouve la valeur minimale. ――La qualité de l'ajustement peut être calculée par variation expliquée.

À la fin

La modélisation est une idée essentielle dans divers domaines tels que l'apprentissage automatique. Le flux est fixe, donc même si cela semble difficile au début, si vous l'essayez vous-même, vous pouvez vous y habituer de manière inattendue. J'aimerais le modéliser, mais j'ai un peu peur et j'espère que cela aidera ces gens. Le code source est répertorié ci-dessous.

`curvefitting.py`



# -*- coding: utf-8 -*-

# import libraries
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import optimize

# model for performance ========
def probabilityCorrect(x, signal):
    # x: model parameters
    # s: signal strength
    
    # probability of correct 
    return 0.5 + (0.5 - x[0])*(1 - np.exp(-(signal/x[1])**x[2]))

def cost_MLE(x, signal, performance, ntr):       
    # compute cost
    c = 0
    for n in np.arange(len(ntr)):
        # model output 
        p = probabilityCorrect(x, signal[n])
        
        # cost as the negative likelihood
        if 0 < p <= 1:
            c += -ntr[n]*(performance[n]*np.log(p) + (1 - performance[n])*np.log(1 - p))
        
    return c 

def cost_LSE(x, signal, performance, ntr):       
    # compute cost
    c = 0
    for n in np.arange(len(ntr)):
        # model output 
        p = probabilityCorrect(x, signal[n])
        
        # cost as the summed squared error
        c += ntr[n]*(performance[n] - p)**2
        
    return c 

# signal strength
signal = [0, 0.1, 0.2, 0.4, 0.8]

# performance
performance = [0.5, 0.73, 0.84, 0.92, 0.98]

# number of trials
ntr = [50, 45, 40, 35, 30]

# initial guess and boundary
x0 = [0, 0.2, 1]
bound = [(0, None),(0.0001, None),(0,None)]

# maximum likelihood
params_MLE = optimize.minimize(cost_MLE,x0,args=(signal,performance,ntr,),method='l-bfgs-b',\
                        jac=None, bounds=bound, tol=None, callback=None,\
                        options={'disp': None, 'maxls': 20, 'iprint': -1,\
                                 'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 15000,\
                                 'ftol': 2.220446049250313e-09, 'maxcor': 10,\
                                 'maxfun': 15000})

# least squared error
params_LSE = optimize.minimize(cost_LSE,x0,args=(signal,performance,ntr,),method='l-bfgs-b',\
                        jac=None, bounds=bound, tol=None, callback=None,\
                        options={'disp': None, 'maxls': 20, 'iprint': -1,\
                                 'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 15000,\
                                 'ftol': 2.220446049250313e-09, 'maxcor': 10,\
                                 'maxfun': 15000})

# compute variance explained
fitted_acc_mle = np.zeros(len(signal))
fitted_acc_lse = np.zeros(len(signal))
for s in np.arange(len(signal)):
    fitted_acc_mle[s] = probabilityCorrect(params_MLE.x,signal[s])
    fitted_acc_lse[s] = probabilityCorrect(params_LSE.x,signal[s])
    
varexp_mle = 1 - (np.var(fitted_acc_mle - performance)/np.var(performance))
varexp_lse = 1 - (np.var(fitted_acc_lse - performance)/np.var(performance))

print('variance explained (MLE): ' + str(varexp_mle))
print('variance explained (LSE): ' + str(varexp_lse))    
    
# visualize
signals = np.linspace(0,0.8,100)
acc_mle = np.zeros(100)
acc_lse = np.zeros(100)

for i in np.arange(100):
    acc_mle[i] = probabilityCorrect(params_MLE.x,signals[i])
    acc_lse[i] = probabilityCorrect(params_LSE.x,signals[i])

fig = plt.figure()
fig.suptitle('psychometric function')

ax = fig.add_subplot(111)

ax.set_xlabel('signal strength')
ax.set_ylabel('performance')    
plt.plot(signal, performance, 'ko')
plt.plot(signals, acc_mle, '-b')
plt.plot(signals, acc_lse, '-r')

[PYTHON] Méthode du carré minimum et méthode d'estimation la plus probable (comparaison par ajustement du modèle)

Tâche

curvefitting.py

modèle

curvefitting.py

Ajustement du modèle (méthode du carré minimum)

curvefitting.py

curvefitting.py

Ajustement du modèle (méthode d'estimation la plus probable)

curvefitting.py

curvefitting.py

Comparaison des résultats d'ajustement

curvefitting.py

curvefitting.py

Résumé récapitulatif

À la fin

curvefitting.py

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`