[PYTHON] Mathematical Examination Part 5 (Examinee paramter estimation)

This is a continuation of Test Mathematical Part 4 (Implementation of Problem paramter Estimate).

Last time, it was "Implementation of problem paramter estimation of 3PL model". This time, it is "estimation of the ability parameter of the examinee".

The environment used is

• python 3.8
• numpy 1.19.2
• matplotlib 3.3.1

is.

theory

Problem setting

When the question parameters $ a, b, c $ are known, the 3PL model test taker parameter $ \ theta_j $ is estimated from the response vector $ u_j $ to the question.

Objective function

Here, we will solve by maximum likelihood estimation. Probability that test taker $ j $ answers correctly to question $ i $ in 3PL model of $ P_ {ij} $, likelihood is when $ Q_ {ij} = 1 --P_ {ij} $

\begin{align}
\Pr\{U_j = u_j|\theta_j, a, b, c\} &= \prod_{i} \Pr\{U_{ij} = u_{ij}|\theta, a, b, c\}\\
&= \prod_{i}P_{ij}^{u_{ij}}Q_{ij}^{1-u_{ij}}
\end{align}

So the log-likelihood is

\begin{align}
\mathcal{L}(u_j, \theta_j, a, b, c) &= \ln \Pr\{U_j = u_j|\theta_j, a, b, c\}\\
&= \sum_{i}u_{ij}\ln P_{ij} + (1-u_{ij})\ln Q_{ij}
\end{align}

It will be. Where the desired parameter $ \ hat {\ theta j} $ is

\begin{align}
\hat{\theta_j} = \arg \max_{\theta}\mathcal{L}(u_j, \theta, a, b, c)
\end{align}

It will be.

Derivative coefficient

Calculate the fine factor to solve the problem. When $ \ partial_ \ theta: = \ frac {\ partial} {\ partial \ theta} $

\begin{align}
\partial_\theta\mathcal{L}(u_j, \theta, a, b, c) &= \sum_i a_i\left(\frac{u_{ij}}{P_{ij}} - 1\right)P_{ij}^* \\
\partial_\theta^2\mathcal{L}(u_j, \theta, a, b, c) &= \sum_i a_i^2\left(\frac{c_iu_{ij}}{P_{ij}^2} - 1\right)P_{ij}^* Q_{ij}^*
\end{align}

It will be. Here, $ P ^ \ * \ _ {ij}, Q ^ \ * \ _ {ij} $ is the probability that the examinee $ j $ will answer the question $ i $ correctly in the 2PL model. From here, the solution should be $ \ partial_ \ theta \ mathcal {L} = 0, \ partial_ \ theta ^ 2 \ mathcal {L} <0 $.

solution

Here, the calculation is performed using Newton Raphson's method. The specific procedure is as follows.

1. Determine the value with a random number as $ \ theta $.
2. $ \ theta ^ {(new)} = \ theta ^ {(old)}-\ partial_ \ theta \ mathcal {L} (\ theta ^ {(old)}) / \ partial_ \ theta ^ 2 \ mathcal { The value is updated sequentially with L} (\ theta ^ {(old)}) $, and ends when $ \ partial_ \ theta \ mathcal {L} <\ delta $.
3. Determine if $ \ partial_ \ theta ^ 2 \ mathcal {L} <0 $ and if the value is divergent, and if the conditions are not met, start over from 1.
4. If the value is not found even after repeating $ N $, no solution is given.

Implementation

I implemented it using python as follows.

import numpy as np
from functools import partial

#Definition of constant delta:Convergence test, Theta_:Divergence judgment, N:Repeat judgment
delta = 0.0001
Theta_ = 20
N = 20
def search_theta(u_i, item_params, trial=0):
    #Returns None when the number of repetitions reaches the specified number
    if trial == N:
        return None
    #Determine theta with random numbers.
    theta = np.random.uniform(-2, 2)
    #NR method
    for _ in range(100):
        P3_i = np.array([partial(L3P, x=theta)(*ip) for ip in item_params])
        P2_i = np.array([partial(L2P, x=theta)(*ip) for ip in item_params])
        a_i = item_params[:, 0]
        c_i = item_params[:, 2]
        res = (a_i *( u_i / P3_i - 1) * P2_i).sum()
        d_res = (a_i * a_i * ( u_i  * c_i/ P3_i / P3_i  - 1) * P2_i * (1 - P2_i)).sum()
        theta -= res/d_res
        #Convergence test
        if abs(res) < delta:
            if d_res < -delta and -Theta_ < theta < Theta_:
                return theta
            else:
                break
    return search_theta(u_i, item_params, trial + 1)

Execution and results

Numerical calculation was performed using the above function. As a problem parameter for execution previous

Was used.

params = np.array(
    [[ 3.49633348,  1.12766137,  0.35744497],
     [ 2.06354365,  1.03621881,  0.20507606],
     [ 1.64406087, -0.39145998,  0.48094315],
     [ 1.47999466,  0.9592384 ,  0.18384673],
     [ 1.44474336,  1.12406269,  0.31475672],
     [ 3.91285332, -1.09218709,  0.18379076],
     [ 1.44498535,  1.50705016,  0.20601461],
     [ 2.37497907,  1.61937999,  0.10503096],
     [ 3.10840278,  1.69962392,  0.22051818],
     [ 1.79969976,  0.06053145,  0.29944448]]
)

Now, regarding the estimation of the examinee paramter, in fact, there are only 10 questions this time, and the possible examination results as input are from [0, 0, ..., 0] to [1, 1, ... There are only 1024 ways of, 1] . Actually, there are only two correct and incorrect responses, so if the number of questions is $ I $, the possible input will be $ 2 ^ I $. Therefore, it is impossible to see the real parameter $ \ theta $ with complete precision. For example, let's say a test taker's parameter is around 0.1224, and the result is[0, 1, 1, 0, 1, 1, 1, 0, 0, 1]. here,

search_theta(np.array([0, 1, 1, 0, 1, 1, 1, 0, 0, 1]), params)

If you calculate, the result will be about 0.8167. If you plot on the graph, it will look like this:

import matplotlib.pyplot as plt
y = []
y_ = []
x = np.linspace(-3, 3, 201)
u_ = np.array([0, 1, 1, 0, 1, 1, 1, 0, 0, 1])
for t in x:
    theta = t
    P3_i = np.array([partial(L3P, x=theta)(*ip) for ip in params])
    P2_i = np.array([partial(L2P, x=theta)(*ip) for ip in params])
    a_i = params[:, 0]
    #Log likelihood
    res = (u_ * np.log(P3_i) + (1 - u_) * np.log(1 - P3_i)).sum()
    #Derivative coefficient
    res_ = (a_i *( u_ / P3_i - 1) * P2_i).sum()
    y.append(res)
    y_.append(res_)


# plot
#Derivative coefficient
plt.plot(x, y_, label="dL")
plt.plot(x, x * 0, label="0")
plt.xlabel("theta")
plt.ylabel("val")
plt.legend()
plt.show()
#Log likelihood
plt.plot(x, y, label="L")
plt.xlabel("theta")
plt.ylabel("val")
plt.legend()
plt.show()

Actually, it was about 0.1224, so it is a fairly high number. This is because I was lucky enough to answer Q2 and Q5 correctly. This is probably due to the low number of questions, so the higher the number of questions, the easier it will be to get a more accurate estimate.

bonus

Binary used for input, for example

n_digits = 10
num = 106
u_ = np.array([int(i) for i in f'{num:0{n_digits}b}'])

Then you can get ʻarray ([0, 0, 0, 1, 1, 0, 1, 0, 1, 0])and so on.0 <= num <= 1023 Now you can see the optimal value for all inputs. However, some, such as[1, 1, ..., 1]`, do not give decent values.

Conclusion

This is the "Mathematical Examination" series that I wrote over 5 articles, but this time it is complete. Thank you for reading this article and the previous articles. The item reaction theory is deep, for example, what about data that does not answer all (toothless)? What if it's not binary in the first place, it's multi-step, or the data comes out in real numbers? Various applications and developments are possible. Also, if I have the opportunity, I would like to study about that and write an article. If you have any questions or concerns, please feel free to comment or request corrections.

References

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