[PYTHON] Examen Mathématiques Partie 4 (Implémentation de l'estimation des paramètres du problème)

Ceci est une suite de Test Mathematics Part 3 (optimisation du modèle 3PL).

La dernière fois, c'était "A propos des mathématiques d'optimisation du modèle 3PL". Cette fois, il s'agit de "A propos de la mise en œuvre du problème d'estimation des paramètres du modèle 3PL".

L'environnement utilisé est

• python 3.8
• numpy 1.19.2
• scipy 1.5.2

est.

Génération de données

Générer divers paramètres et matrices de réaction d'item pour l'expérience comme effectué dans Test Mathématiques Partie 1 (Définition des questions et génération de données).

import numpy as np
from functools import partial

#Gardez un petit ε pour qu'il ne diverge pas dans le calcul numérique
epsilon =  0.0001
#3 Définition du modèle logistique des paramètres
def L3P(a, b, c, x):
    return c + (1 - epsilon - c) / (1 + np.exp(-  a * (x - b)))

#2 Définition du modèle logistique des paramètres. Afin d'unifier le traitement, c est également pris comme argument.
def L2P(a, b, c, x):
    return (1 - epsilon) / (1 + np.exp(-  a * (x - b)))

#Définition du paramètre du modèle
#a est un nombre réel positif,b est un nombre réel,c doit être supérieur à 0 et inférieur à 1

a_min = 0.7
a_max = 4

b_min = -2
b_max = 2

c_min = 0
c_max = .4

#Combien de questions, combien de personnes, 10 questions 4000 personnes ci-dessous
num_items = 10
num_users = 4_000

#Générer un paramètre de problème
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

#Génération de paramètres candidats
user_params = np.random.normal(size=num_users)

#Créer une matrice de réaction d'élément, élément 1(Bonne réponse)Ou 0(Mauvaise réponse)
#Dans la ligne i et la colonne j, comment le candidat j a-t-il réagi à la question i?
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]
    )
)

Comme précédemment, nous utiliserons $ i $ comme indice pour la question et $ j $ comme indice pour le candidat. Nous porterons un grand nombre de candidats à 1 000 ou plus. En effet, il est empiriquement connu que l'estimation est stable dans une certaine mesure lorsqu'il y a 100 personnes ou plus pour l'estimation 1PL, 300 ou plus pour l'estimation 2PL et 1000 ou plus pour l'estimation 3PL. Si vous êtes intéressé, veuillez changer ici et expérimenter.

Ici, comme paramètre de problème

Eu Le but est d'estimer cela.

Prétraitement des données

Pour stabiliser la précision de l'estimation avant de faire l'estimation

• Problèmes avec un taux de réponse correct trop élevé
• Problèmes avec un pourcentage trop faible de bonnes réponses
• Problèmes de faible corrélation avec les scores bruts

Est supprimé de la cible d'estimation. En particulier

#Éliminez les problèmes difficiles à estimer.
#Le taux de réponse correct est trop élevé(> 0.9), Trop bas(0.1 <), Il y a trop peu de corrélation avec le score brut(< 0.3)Supprimer le problème
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

#Redéfinir la matrice de réaction de l'élément
irm_ij = ir_matrix_ij[filter_total]
num_items, num_users = irm_ij.shape

Préparation des fonctions et des constantes pour l'estimation

Dernière fois Comme vous pouvez le voir, chaque problème $ 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}

Est-ce que l'étape E pour calculer

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

Calculer

\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}

Le pas M est de mettre à 0. Ici, $ X_k $ et $ g_k = g (X_k) $ sont préparés comme des points représentatifs de $ \ theta $, car la marginalisation et l'intégration ci-dessus sont résolues d'une manière produit segmenté.

Par conséquent, préparez ce qui suit.

#Définissez la plage possible de paramètres du candidat.
X_k = np.linspace(-4, 4, 41)
#Définissez la distribution des paramètres du candidat. Ici, la distribution normale de scipy est utilisée.
from scipy.stats import norm
g_k = norm.pdf(X_k)
#Fonction étape E
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
#Fonction de score pour l'étape M
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)

Effectuer une estimation

Maintenant, effectuons l'estimation. Étant donné que l'estimation dépend du paramètre initial et que la stabilité ne devrait pas être aussi bonne, préparez au hasard certains paramètres initiaux et adoptez la valeur médiane parmi les plus stables comme résultat de l'estimation. Pour l'optimisation de l'étape M, nous utiliserons minimiser de scipy.

from scipy.optimize import minimize
#Définissez les contraintes pour minimiser.
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],
    ]
}
#Préparez un paramètre pour la génération initiale des paramètres.
a_min, a_max = 0.1, 8.0
b_min, b_max = -4.0, 4.0
c_min, c_max = epsilon, 0.6

#Oarameter pour l'exécution de l'estimation
#Condition de fin de répétition de l'algorithme EM
delta = 0.001
#Nombre maximum de répétitions de l'algorithme EM
max_num_of_itr = 1000

#Calculer plusieurs fois la stabilité numérique et adopter la valeur médiane parmi les stables
p_data = []
for n_try in range(10):
    #Définissez la valeur initiale de l'estimation.
    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 :Calcul du paramètre exp
        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 = []
        #Résoudre les problèmes d'optimisation pour chaque problème
        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)
        #Le calcul se termine lorsque la différence moyenne par rapport à l'heure précédente tombe en dessous d'une certaine valeur
        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, :]
    #Élimine les extrêmes dans les résultats de calcul
    filter_1 = t_[:, 1] < b_max - epsilon 
    filter_2 = t_[:, 1] > b_min + epsilon
    #La médiane restante est utilisée comme résultat du calcul.
    result_.append(np.median(t_[filter_1 & filter_2], axis=0))
    
result = np.array(result_)

Résultat expérimental

Cette fois, j'ai obtenu ce qui suit à la suite d'un calcul. ** Résultat du calcul (objectif) **

Il y a des valeurs étranges telles que b (difficulté) en Q3, mais elles peuvent être estimées avec une bonne précision.

la prochaine fois

La prochaine fois, nous évaluerons le paramètre de capacité $ \ theta_j $ du candidat en utilisant le résultat du paramètre en question. Examen Mathématiques Partie 5 (Estimation du paramètre de l'examinateur)

Les références

• Livre Item Reaction Theory / Theory-Test Mathematics (Statistical Library) (amazon.co.jp)
• Livre Introduction to Item Reaction Theory (2nd Edition) (Statistics Library) (amazon.co.jp)