[PYTHON] Prüfungsmathematik Teil 4 (Implementierung der Problemparameterschätzung)

Dies ist eine Fortsetzung von Testmathematik Teil 3 (3PL-Modelloptimierung).

Letztes Mal war es "Über die Mathematik der Optimierung des 3PL-Modells". Diesmal geht es um "die Implementierung der Problemparameterschätzung des 3PL-Modells".

Die verwendete Umgebung ist

• python 3.8
• numpy 1.19.2
• scipy 1.5.2

ist.

Datengenerierung

Generieren Sie verschiedene Parameter und Item-Reaktionsmatrizen für das Experiment, wie in Testmathematik Teil 1 (Frageneinstellung und Datengenerierung) durchgeführt.

import numpy as np
from functools import partial

#Halten Sie ein kleines ε, damit es bei der numerischen Berechnung nicht abweicht
epsilon =  0.0001
#3 Definition des Parameterlogistikmodells
def L3P(a, b, c, x):
    return c + (1 - epsilon - c) / (1 + np.exp(-  a * (x - b)))

#2 Definition des Parameterlogistikmodells. Um die Verarbeitung zu vereinheitlichen, wird auch c als Argument verwendet.
def L2P(a, b, c, x):
    return (1 - epsilon) / (1 + np.exp(-  a * (x - b)))

#Definition des Modellparameters
#a ist eine positive reelle Zahl,b ist eine reelle Zahl,c sollte größer als 0 und kleiner als 1 sein

a_min = 0.7
a_max = 4

b_min = -2
b_max = 2

c_min = 0
c_max = .4

#Wie viele Fragen, wie viele Personen, 10 Fragen 4000 Personen unten
num_items = 10
num_users = 4_000

#Problemparameter generieren
item_params = np.array(
    [np.random.uniform(a_min, a_max, num_items),
     np.random.uniform(b_min, b_max, num_items),
     np.random.uniform(c_min, c_max, num_items)]
).T

#Generierung von Kandidatenparametern
user_params = np.random.normal(size=num_users)

#Erstellen Sie eine Artikelreaktionsmatrix, Element 1(Richtige Antwort)Oder 0(Falsche Antwort)
#Wie hat der Prüfling j in Zeile i und Spalte j auf Frage i reagiert?
ir_matrix_ij = np.vectorize(int)(
    np.array(
        [partial(L3P, *ip)(user_params) > np.random.uniform(0, 1, num_users) for ip in item_params]
    )
)

Nach wie vor verwenden wir $ i $ als Index für die Frage und $ j $ als Index für den Kandidaten. Wir werden eine große Anzahl von Prüflingen auf 1000 oder mehr bringen. Dies liegt daran, dass empirisch bekannt ist, dass die Schätzung bis zu einem gewissen Grad stabil ist, wenn 100 oder mehr Personen für die 1PL-Schätzung, 300 oder mehr für die 2PL-Schätzung und 1000 oder mehr für die 3PL-Schätzung anwesend sind. Wenn Sie interessiert sind, ändern Sie bitte hier und experimentieren Sie.

Hier als Problemparameter

Habe Der Zweck ist es, dies abzuschätzen.

Datenvorverarbeitung

Stabilisierung der Genauigkeit der Schätzung vor der Schätzung

• Probleme mit einer zu hohen korrekten Antwortrate
• Probleme mit zu geringem Prozentsatz korrekter Antworten
• Probleme mit geringer Korrelation mit Rohwerten

Wird aus dem Schätzziel entfernt. Speziell

#Beseitigen Sie schwer abzuschätzende Probleme.
#Die richtige Antwortrate ist zu hoch(> 0.9), Zu niedrig(0.1 <)Es gibt zu wenig Korrelation mit der Rohbewertung(< 0.3)Entfernen Sie das Problem
filter_a = ir_matrix_ij.sum(axis=1) / num_users < 0.9
filter_b = ir_matrix_ij.sum(axis=1) / num_users > 0.1
filter_c = np.corrcoef(np.vstack((ir_matrix_ij, ir_matrix_ij.sum(axis=0))))[num_items][:-1] >= 0.3
filter_total = filter_a & filter_b & filter_c

#Definieren Sie die Item-Reaktionsmatrix neu
irm_ij = ir_matrix_ij[filter_total]
num_items, num_users = irm_ij.shape

Vorbereitung von Funktionen und Konstanten zur Schätzung

Letztes Mal Wie Sie sehen können, ist jedes Problem $ i $

\begin{align}
r_{i\theta} &:= 
\sum_{j}u_{ij} \Pr\{\theta| U_j = u_j, a^{old}, b^{old}, c^{old}\} \\
f_\theta &= 
\sum_{j} \Pr\{\theta| U_j = u_j, a^{old}, b^{old}, c^{old}\} 
\end{align}

Ist der E-Schritt zu berechnen

\begin{align}
R_{i\theta} :=\frac{r_{i\theta}}{P_{i\theta}} - f_\theta
\end{align}

Berechnung

\begin{align}
\partial^a\mathcal{\tilde{L}}(a, b, c) 
&= \int R_{i\theta}(\theta - b_i)P_{i\theta}^*d\theta\\
\partial^b\mathcal{\tilde{L}}(a, b, c) 
&= -a_i \int R_{i\theta}P_{i\theta}^*d\theta\\
\partial^c\mathcal{\tilde{L}}(a, b, c) 
&= \frac{1}{(1 - c_i)} \int R_{i\theta} d\theta
\end{align}

Der M-Schritt besteht darin, auf 0 zu setzen. Hier werden $ X_k $ und $ g_k = g (X_k) $ als repräsentative Punkte von $ \ theta $ vorbereitet, da die Marginalisierung und die obige Integration in einer segmentierten Produktweise gelöst werden.

Bereiten Sie daher Folgendes vor.

#Definieren Sie den möglichen Bereich der Testteilnehmerparameter.
X_k = np.linspace(-4, 4, 41)
#Definieren Sie die Verteilung der Testteilnehmerparameter. Hier wird die Normalverteilung von scipy verwendet.
from scipy.stats import norm
g_k = norm.pdf(X_k)
#E Schrittfunktion
def get_exp_params(irm_ij, g_k, P3_ik):
    Lg_jk = np.exp(irm_ij.T.dot(np.log(P3_ik)) + (1 - irm_ij).T.dot(np.log(1 - P3_ik)))* g_k
    n_Lg_jk = Lg_jk / Lg_jk.sum(axis=1)[:, np.newaxis]
    f_k = n_Lg_jk.sum(axis=0)
    r_ik = irm_ij.dot(n_Lg_jk)
    return f_k, r_ik
#Bewertungsfunktion für M-Schritt
def score_(param, f_k, r_k, X_k):
    a, b, c = param
    P3_k = partial(L3P, a, b, c)(X_k)
    P2_k = partial(L2P, a, b, c)(X_k)
    R_k = r_k / P3_k - f_k
    v = [
        ((X_k - b) * R_k * P2_k).sum(),
        - a * (R_k * P2_k).sum(),
        R_k.sum() / (1 - c)
    ]
    return np.linalg.norm(v)

Schätzung durchführen

Lassen Sie uns nun die Schätzung durchführen. Da die Schätzung vom Anfangsparameter abhängt und erwartet wird, dass die Stabilität nicht so gut ist, bereiten Sie einige Anfangsparameter zufällig vor und übernehmen Sie den Medianwert unter den stabilen als Schätzergebnis. Zur Optimierung des M-Schritts verwenden wir die Minimierung von Scipy.

from scipy.optimize import minimize
#Definieren Sie Einschränkungen für die Minimierung.
cons_ = {
    'type': 'ineq',
    'fun': lambda x:[
        x[0] - a_min,
        a_max - x[0],
        x[1] - b_min,
        b_max - x[1],
        x[2] - c_min,
        c_max - x[2],
    ]
}
#Bereiten Sie einen Parameter für die anfängliche Parametergenerierung vor.
a_min, a_max = 0.1, 8.0
b_min, b_max = -4.0, 4.0
c_min, c_max = epsilon, 0.6

#Oarameter für die Ausführung der Schätzung
#Wiederholungsbeendigungsbedingung des EM-Algorithmus
delta = 0.001
#Maximale Anzahl von Wiederholungen des EM-Algorithmus
max_num_of_itr = 1000

#Berechnen Sie mehrmals die numerische Stabilität und übernehmen Sie den Medianwert unter den stabilen
p_data = []
for n_try in range(10):
    #Definieren Sie den Anfangswert der Schätzung.
    item_params_ = np.array(
        [np.random.uniform(a_min, a_max, num_items),
         np.random.uniform(b_min, b_max, num_items),
         np.random.uniform(c_min, c_max, num_items)]
    ).T

    prev_item_params_ = item_params_
    for itr in range(max_num_of_itr):
        # E step :Berechnung von exp param
        P3_ik = np.array([partial(L3P, *ip)(X_k) for ip in item_params_])
        f_k, r_ik = get_exp_params(irm_ij, g_k, P3_ik)
        ip_n = []
        #Lösen Sie Optimierungsprobleme für jedes Problem
        for item_id in range(num_items):
            target = partial(score_, f_k=f_k, r_k=r_ik[item_id], X_k=X_k)
            result = minimize(target, x0=item_params_[item_id], constraints=cons_, method="slsqp")         
            ip_n.append(list(result.x))

        item_params_ = np.array(ip_n)
        #Die Berechnung endet, wenn die durchschnittliche Differenz zum vorherigen Zeitpunkt einen bestimmten Wert unterschreitet
        mean_diff = abs(prev_item_params_ - item_params_).sum() / item_params_.size
        if mean_diff < delta:
            break
        prev_item_params_ = item_params_

    p_data.append(item_params_)
    
p_data_ = np.array(p_data)
result_ = []
for idx in range(p_data_.shape[1]):
    t_ = np.array(p_data)[:, idx, :]
    #Beseitigen Sie Extreme in den Berechnungsergebnissen
    filter_1 = t_[:, 1] < b_max - epsilon 
    filter_2 = t_[:, 1] > b_min + epsilon
    #Der verbleibende Median wird als Berechnungsergebnis verwendet.
    result_.append(np.median(t_[filter_1 & filter_2], axis=0))
    
result = np.array(result_)

Versuchsergebnis

Dieses Mal habe ich als Ergebnis der Berechnung Folgendes erhalten. ** Berechnungsergebnis (Zweck) **

Es gibt einige seltsame Werte wie b (Schwierigkeit) in Q3, aber sie können mit ungefähr guter Genauigkeit geschätzt werden.

nächstes Mal

Das nächste Mal werden wir den Fähigkeitsparameter $ \ theta_j $ des Prüflings anhand des Ergebnisses des betreffenden Parameters schätzen. Prüfungsmathematik Teil 5 (Schätzung des Prüflingsparameters)

Verweise

• Buch Item Reaction Theory / Theory-Test Mathematics (Statistische Bibliothek) (amazon.co.jp)
• Buch Einführung in die Theorie der Gegenstandsreaktion (2. Ausgabe) (Statistikbibliothek) (amazon.co.jp)