Die Welt ist voller verschiedener Phänomene. Um sie zu verstehen, muss ein Modell erstellt werden, das ihr Verhalten erklärt. Wenn das Modell die Daten zur Beobachtung eines bestimmten Phänomens vorhersagen kann, kann gesagt werden, dass das Modell ein bestimmtes Phänomen versteht.

Es gibt zwei Möglichkeiten, das Modell an die Daten anzupassen.

** Methode der kleinsten Quadrate **
** Maximum-Likelihood-Methode **

Welches sollte wann verwendet werden? Dieses Mal werde ich versuchen, ** Kurvenanpassung ** anhand eines einfachen Beispiels.

Aufgabe

Nehmen wir das Beispiel eines Studiums des klassischen Bewusstseins. Lassen Sie das Motiv vor dem Monitor sitzen und den visuellen Reiz auf dem Monitor präsentieren. Visuelle Reize umfassen Signale und Rauschen, und die Probanden werden gebeten zu melden, ob sie das Signal für jeden Versuch gesehen haben.

Intuitiv ist es umso einfacher zu erkennen, je höher das Signalverhältnis als das Rauschen ist. Wenn sich die Signalstärke beispielsweise auf der horizontalen Achse befindet und die Erkennungsleistung des Subjekts auf der vertikalen Achse liegt, ist dies wie folgt.

Dies wird als psychometrische Funktion bezeichnet, spielt aber jetzt keine Rolle, und Sie können die ** nichtlineare ** Beziehung deutlich erkennen. Erstellen wir ein Modell, das diese nichtlineare Beziehung darstellt, und passen es in die Daten ein. Wenn Sie dies tun können, können Sie die Erkennungsleistung zu diesem Zeitpunkt vorhersagen, auch wenn die Signalstärke im Experiment nicht verwendet wird. Wenn Sie das können, verstehen Sie dieses nichtlineare Phänomen.

Die diesmal verwendeten Daten sind nachstehend zusammengefasst.

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

Die für die Anpassung erforderlichen Informationen sind die Leistung, die der Stärke jedes Signals und der Anzahl der Versuche mit den Stimuli entspricht, als die Probanden aufgefordert wurden, alle 200 Versuche durchzuführen.

Dieses Mal verwenden wir Python. Importieren Sie also die erforderlichen Bibliotheken und fügen Sie die Daten in die Liste ein.

`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]

Modell

Es gibt viele nichtlineare Modelle, aber wir haben festgestellt, dass ** kumulativer Weibull ** für Erkennungsaufgaben wie diese gut funktioniert.

Wenn die Stärke eines Signals $ s $ als Stimulus dargestellt wird, wird die Wahrscheinlichkeit, dass das Subjekt richtig $ p $ antwortet, wie folgt ausgedrückt.


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

$ f $ ist $ f (x) = (1 --e ^ {-x}) $. Es gibt drei freie Parameter, die die Form des kumulativen Löschvorgangs bestimmen: $ \ lambda $, $ \ theta $, $ \ beta $. Da wir die Signalstärke $ s $ und die Subjektleistung $ p $ im Experiment kennen, werden wir diese verwenden, um ** freie Parameter zu finden, die zu den Daten passen **.

Bei der eigentlichen Programmierung werden alle freien Parameter in eine Listenvariable namens $ x $ eingefügt. $ x [0] = \ lambda $, $ x [1] = \ theta $, $ x [2] = \ beta $ und so weiter.

`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]))

Modellanpassung (Minimum-Square-Methode)

Versuchen Sie zuerst die ** Methode der kleinsten Quadrate **. Die Methode der kleinsten Quadrate ist eine Methode zum Auffinden freier Parameter wie ** "Die Differenz zwischen dem vorhergesagten Wert des Modells und den tatsächlichen Daten wird minimiert" **. Das Quadrat wird einfach verwendet, weil Sie nicht über das Zeichen nachdenken müssen.

Um die freien Parameter zu finden, erstellen Sie eine sogenannte ** Kostenfunktion **. Der Ablauf besteht darin, ihn in den Algorithmus zu werfen und den Algorithmus einen freien Parameter finden zu lassen, der seine Verlustfunktion minimiert.

Mit der Methode der kleinsten Quadrate ist die Verlustfunktion sehr intuitiv. Der vorhergesagte Wert des Modells und die quadratische Differenz der tatsächlichen Daten sollten basierend auf der Stärke jedes Signals berechnet und zu einem numerischen Wert addiert werden.

Das Folgende ist ein Implementierungsbeispiel.

`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

Nachdem wir nun eine Verlustfunktion haben, werfen Sie sie in den Algorithmus, um einen freien Parameter zu finden, der den Mindestwert annimmt. In Python enthält scipy ein Paket namens optimize. Verwenden Sie dieses Paket.

`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})

Die Argumente mögen verwirrend sein, aber im Grunde genommen Offizielle Website Dies ist das Standard-Kopieren und Einfügen in. Bitte beachten Sie, dass ** Sie den Anfangswert des freien Parameters $ x0 $ ** festlegen müssen.

Der diesmal verwendete Algorithmus ist L-BFGS-B, ein speicherfreundlicher Allzweckalgorithmus. Jeder Algorithmus wiederholt jedoch numerische Berechnungen (numerisch), um den Mindestwert zu ermitteln, sodass der Anfangswert nahe an der Antwort liegt. Es ist ideal. Wenn der Anfangswert irrelevant ist, kann die Algorithmusberechnung ** nicht konvergieren ** und die gewünschte Anpassung kann nicht erreicht werden. Ich habe keine andere Wahl, als verschiedene Dinge auszuprobieren, aber dieses Mal lautet der Anfangswert empirisch "Nun, das wird passieren".

Eine andere Sache ist, dass Sie die Grenzen der ** freien Parameter $ \ lambda $, $ \ theta $, $ \ beta $ ** voreinstellen können. Mit den richtigen Entfernungseinstellungen kann der Algorithmus vermeiden, durch eine große Anzahl von Ozeanen zu wandern. Dieses Mal sind alle freien Parameter positiv und wir möchten nicht, dass $ x [1] = \ theta $ $ 0 $ ist, da es sich um den Nenner des Modells handelt. Deshalb setzen wir die Grenze wie oben. Es ist eingestellt.

Jetzt findet der Algorithmus L-BFGS-B den Wert des freien Parameters, der die Verlustfunktion minimiert.

Modellanpassung (wahrscheinlichste Schätzmethode)

Versuchen Sie zum Vergleich die ** Maximum-Likelihood-Methode **. Kurz gesagt, die wahrscheinlichste Schätzmethode ist eine Methode zum Auffinden freier Parameter wie ** "Maximieren der Wahrscheinlichkeit (Wahrscheinlichkeit), dass das Modell tatsächliche Daten vorhersagt, wenn ein bestimmter freier Parameter angegeben wird" **. ist.

Ich habe ein Experiment durchgeführt und einige Daten erhalten. Die Daten stammen aus der Wahrscheinlichkeitsverteilung eines Modells, die durch die freien Parameter des Modells bestimmt wird. Das Grundprinzip besteht also darin, einen freien Parameter zu finden, der die Wahrscheinlichkeit (Plausibilität, Wahrscheinlichkeit, Wahrscheinlichkeit ** des tatsächlichen Datenwerts, der aus dem Modell kommt, maximiert.

Insbesondere wird die ** Wahrscheinlichkeit durch die gemeinsame Wahrscheinlichkeit von $ P (data_i | x) $ definiert, die dem Modell durch den freien Parameter $ x $ für jedes Beobachtungsereignis $ data_i $ gegeben wird. **. Wie wir in der High School gelernt haben, wird die gleichzeitige Wahrscheinlichkeit als Produkt der Wahrscheinlichkeit ausgedrückt, dass jedes Ereignis eintreten wird, wenn jedes Beobachtungsereignis unabhängig auftritt.


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

Ich möchte eine Funktion erstellen, die dies maximiert, aber die einzige Möglichkeit, freie Parameter zu finden, besteht darin, eine Kostenfunktion zu erstellen, deren Mindestwert von einem Algorithmus ermittelt werden kann. Da ich diesmal jedoch den Maximalwert haben möchte, übergebe ich die Wahrscheinlichkeit mit einem negativen Vorzeichen an den Algorithmus.

Ein weiterer Grund ist, dass die Wahrscheinlichkeit das Produkt vieler Wahrscheinlichkeiten ist. Wenn also die Anzahl der Beobachtungen groß ist, nähert sich der Wert allmählich 0, was es für den Algorithmus schwierig macht, den Maximalwert zu finden. Daher verwenden wir den mathematischen Trick der Protokollierung. Wie ich in der High School gelernt habe, kann das Produkt von Elementen durch Erstellen eines Protokolls in die Summe der Elemente umgewandelt werden.

Infolgedessen ist die Verlustfunktion, die durch das wahrscheinlichste Schätzverfahren in den Algorithmus geworfen wird, wie folgt, wobei das Protokoll genommen und minimiert wird.


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

Das Folgende ist ein Implementierungsbeispiel.

`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

Der vorhergesagte Wert des Modells muss zur Erleichterung der Protokollierung positiv sein. Da es sich um eine Wahrscheinlichkeit handelt, kann sie auch nicht größer als 1 sein. Ich setze diese beiden Bedingungen in $ if $. Beachten Sie, dass in diesem Beispiel die Wahrscheinlichkeit, dass das Subjekt richtig antwortet, als $ p $ modelliert wird, sodass die Wahrscheinlichkeit, dass das Subjekt nicht richtig antwortet, $ 1-p $ beträgt.

Nachdem wir eine Verlustfunktion haben, werfen wir sie in den Algorithmus.

`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})

Vergleich der Anpassungsergebnisse

Was ist mit dem Anpassungsergebnis? Lassen Sie uns ein Diagramm erstellen.

`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')

Tragen Sie die mit der Methode der kleinsten Quadrate erhaltenen Fittings in Rot und die mit der wahrscheinlichsten Schätzmethode erhaltenen Fittings in Blau auf.

…… Es ist fast das gleiche w Beide sind ordentlich montiert.

In der Tat kann mit dem Index ** Variation erklärt ** quantifiziert werden, wie gut die Anpassung durchgeführt wird. Das Folgende ist offiziell.


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

Wir quantifizieren, inwieweit die Variation des vorhergesagten Werts die Variation der tatsächlichen Daten erklärt hat. Berechnen wir das.

`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))

Als Ergebnis der Berechnung,

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

ist geworden. Es scheint, dass die Methode der kleinsten Quadrate etwas besser ist, aber der Unterschied ist vernachlässigbar. Vielleicht möchten Sie vorerst beides ausprobieren und dasjenige auswählen, das besser passt.

Zusammenfassung Zusammenfassung

--Verwenden Sie die Methode der kleinsten Quadrate oder die Methode der Schätzung mit der größten Wahrscheinlichkeit, um das Modell anzupassen.

Der Ablauf des Erstellens eines Modells, des Erstellens einer Verlustfunktion und des Werfens in einen Algorithmus, der den Mindestwert findet. ――Die Qualität der Anpassung kann anhand der erklärten Variation berechnet werden.

Am Ende

Modellierung ist eine wesentliche Idee in verschiedenen Bereichen wie dem maschinellen Lernen. Der Fluss ist fest. Selbst wenn es auf den ersten Blick schwierig aussieht, können Sie sich unerwartet daran gewöhnen, wenn Sie es selbst versuchen. Ich würde es gerne modellieren, aber ich habe ein wenig Angst und hoffe, dass es diesen Menschen hilft. Der Quellcode ist unten aufgeführt.

`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] Minimum-Square-Methode und wahrscheinlichste Schätzmethode (Vergleich durch Modellanpassung)

Aufgabe

curvefitting.py

Modell

curvefitting.py

Modellanpassung (Minimum-Square-Methode)

curvefitting.py

curvefitting.py

Modellanpassung (wahrscheinlichste Schätzmethode)

curvefitting.py

curvefitting.py

Vergleich der Anpassungsergebnisse

curvefitting.py

curvefitting.py

Zusammenfassung Zusammenfassung

Am Ende

curvefitting.py

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`

`curvefitting.py`