1、IRT模型概述
IRT（item response theory 項目反映理論）模型。IRT模型用來描述被試者能力和項目特性之間的關系。在現實生活中，由於被試者的能力不能通過可觀測的數據進行描述，所以IRT模型用一個潛變量 $ \theta $ 來表示，並考慮與項目相關的一組參數來分析正確回答測試項目的概率。目前常見的IRT模型有2-PL模型和3-PL模型。其具體表達式如下：
2-PL模型的表達式如下：
$ p_{i,j}(\theta_i) = \frac {1} {1 + \exp\,[-D\,a_j\,(\theta_i - b_j)]} $
其種 $ \theta_i $ 是被試者能力的參數，$ a_j $ 和 $ b_j $ 分別代表的是題目的區分度參數和難度參數，D是為1.7的常數。
3-PL模型是在模型中引入了預測參數 $ c_j $ ，該參數的含義是描述被試者在沒有任何先驗知識的情況下，回答正確某項目的概率。常見的例子有學生在做選擇題時，即使對該問題沒有任何相關知識的獲取，也有一定的概率答對該題目。
3-PL模型的表達式如下：
$ p_{i,j}(\theta_i) = c_j + \frac {1 - c_j} {1 + \exp\,[-D\,a_j\,(\theta_i - b_j)]} $
$ c_j $ 表示的是預測參數，其余的參數含義和2-PL模型中的一致。
IRT模型滿足三條基本假設：
1）潛在單調性，IRT模型是連續嚴格單調的函數。
2）條件獨立性，IRT模型認為給定 $ \theta_i $，對於第 $ i $ 個被試者，$ Y_{i,j} $ 是獨立的；而對於不同的被試者，其各自的答案 $ Y_i $ 也是相互獨立的。
3）單維性假設，IRT模型認為某個測試的所有項目都是測量同一個潛在特質。
IRT模型的曲線又稱為項目相應曲線或者項目特征曲線（ICC），曲線圖如下：
區分度參數 $ a $ 的值是ICC曲線拐點處的斜率，即斜率的最大值。項目難度參數 $ b $ 一般表示在ICC曲線最陡的那一點所對應的 $ \theta $ 的值。改變 $ b $ 的值不會改線曲線的形狀，只會導致曲線左右移動。預測參數 $ c $ 代表了ICC曲線的下限，如上圖 $ c $ 的值為0.
在這里為了方便討論，我們用2-PL模型來作為示例。而實際情況中為了簡化模型，也可以將預測參數設置為一個常數，例如對於主觀題設置為0，對於選擇題設置為0.25.
2、EM算法求解IRT模型
EM算法的具體講解見這一篇。EM算法是一種迭代算法，用於含有隱變量的概率模型的極大似然估計，或者說是極大后驗概率估計。EM算法分為 E步和 M步
E步：計算聯合分布的條件概率期望。
M步：極大化對數似然函數的條件期望求解參數。
首先來分析我們的觀測數據和隱變量。
觀測數據：我們把被試者 $ i $ 對項目 $ j $ 的作答反映記為 $ y_{i, j} $ ，又記向量 $ y_i = (y_{i1}, y_{i2}, ....., y_{im}) $ ，稱為觀測數據，其中 $ i = 1, 2, ..., N, j = 1, 2, ..., m $。對於兩級計分模型（在這里我們只考慮兩級計分模型，即模型輸出的結果只有二分類），$ y_{ij} $ 的取值有0 和 1，分別表示被試者答錯和答對題目。除了兩級計分模型，還會有多級計分模型（即輸出的結果為多分類）。
隱變量（缺失數據）：我們把每個被試者的潛在的不可觀測的能力值稱為缺失數據，記為 $ \theta = (\theta_1, \theta_2, ......, \theta_N) $ ，其中 $ \theta_i $ 是被試者 $ i $ 的潛在能力值。
完全數據：完全數據對每一個被試者來說就是觀察數據加缺失數據，記作 $ [(y_1, \theta_1), (y_2, \theta_2), ..., (y_N, \theta_N)] $ 。
將EM算法應用到IRT模型中來，則E步和M步可以描述為：
E步：即在給定缺失數據的分布，觀察數據和參數初值時，求完全數據的對數似然函數的條件期望。
M步：即使用E步計算出的完全數據充分統計量的條件期望值，極大化完全數據的對數似然函數的條件期望求解參數的值。
不斷的循環迭代E步和M步，直到參數估計收斂。
在IRT模型當中，我們通常認為能力參數 $ \theta $ 是連續隨機變量，故可以取任意值。在EM算法估計參數的過程中，我們是能力參數 $ \theta $ 為離散分布。能力值只能取 $ q_1, q_2, ......, q_K $, K個值中的一個，且 $ P(\theta = q_k) = \pi_k $ 。
在給定了反應矩陣 $ Y $ 的情況下，假設項目參數 $ \Delta = [\delta_1, ......, \delta_J] $ ，能力分布參數 $ \pi = (\Pi_1, ......, \Pi_K) $ 。則E步中計算的條件期望推到如下：
對上面的公式做一些轉換表示為：
其中 $ \phi(\Delta) $ 和 $ \psi(\pi) $ 的表達式如下：
其中：
在上面式子中的 $ n_k^{(s)} $ 可以理解為在 $ N $ 個被試者中能力水平為 $ q_k $ 的被試數目的期望（即能力為 $ q_k $ 的被試者的個數），$ r_{jk}^{(s)} $ 可以理解為在 $ N $ 個被試者中具有能力水平 $ q_k $ 的被試者答對第 $ j $ 個項目的個數。 $ n_k^{(s)} $ 和 $ r_{jk}^{(s)} $ 都是人工數據。
應用EM算法的第 $ s $ 步迭代中
E步：利用第 $ s-1 $ 步得到的參數估計值 $ \Delta^{(s)} $ 和 $ \pi_k^{(s)} $ 計算 $ n_k^{(s)} $ 和 $ r_{jk}^{(s)} $ （首次迭代時初始化 $ \Delta^{(0)} $ 和 $ \pi_k^{(0)} $ 的值）。
M步：將E步計算出的 $ n_k^{(s)} $ 和 $ r_{jk}^{(s)} $ 代入到 $ \phi(\Delta) $ 和 $ \psi(\pi) $ 中，對兩項分別極大化可得參數 $ \Delta $ 和 $ \pi $ 的估計值 $ \Delta^{(s+1)} $ 和 $ \pi_k^{(s+1)} $ .具體的求解如下：
$ \pi_k^{(s+1)} $ 的值可以直接用上述表達式求出
而對於 $ \Delta^{(s+1)} $ 在這里要用牛頓-拉弗遜方法求解。
EM算法估計IRT模型的步驟如下：
1）E步：
首先確定 $ q_k $ 和 $ \pi_k $ ；用前一次迭代參數 $ \Delta^{(s)} $ 和 $ \pi_k^{(s)} $ 求出 $ n_k^{(s)} $ 和 $ r_{jk}^{(s)} $ 。
2）M步：
計算 $ \delta_j^{(s+1)} $ 和 $ \pi^{(s+1)} $ .
3) 重復E步和M步直到項目參數收斂為止。
3、MCMC算法求解IRT模型
MCMC算法的具體詳解見這一篇，MCMC算法全稱是馬爾科夫-蒙特卡洛算法。MCMC的基礎理論為馬爾可夫過程，在MCMC算法中，為了在一個指定的分布上采樣，根據馬爾可夫過程，首先從任一狀態出發，模擬馬爾可夫過程，不斷進行狀態轉移，最終收斂到平穩分布。用MCMC算法在這里估計參數，事實上就是建立一條參數的馬爾科夫鏈，根據參數的狀態轉移矩陣來輸出馬爾科夫鏈上的樣本，當馬爾科夫鏈到一定長度時就會開始收斂於某一值。
首先來看看MCMC算法在2-PL模型上的參數估計步驟：
1）取模型參數的先驗分布：$ \theta~N(0, 1)，\log(a)~N(0, 1)，b~N(0, 1) $ ，則初始化被試能力參數，項目參數的初始值 $ \theta_0 = 0，a_0 = 1，b_0 = 0 $ 。
2）根據項目參數初值 $ a_0，b_0 $ ，估計被試者的能力參數 $ \theta_1 $ 。
各被試的能力參數 $ \theta_* $ 獨立地從建議性分布 $ q_\theta $ 中選取，我們去被試能力參數的建議分布為 $ \theta_*~N(\theta_0, c_\theta^2) $ 。其中 一般$ c_\theta = 1.1 $ 。
計算從狀態 $ \theta_0 $ 轉移到狀態 $ \theta_1 $ 的接受概率 $ a(\theta_0, \theta_*) = min(1, R_\theta^0) $，由它來決定是否發生狀態的轉移。
其中
生成隨機數 $ r_1~U(0, 1) $ ，比較 $ r_1 $ 與接受概率 $ a(\theta_0, \theta_*) $ 的大小，進行狀態轉移判斷：
若 $ a(\theta_0, \theta_*) \ge r_1 $ ，則有 $ \theta_1 = \theta_* $ ；否則 $ \theta_1 = \theta_* $ 。
3）根據步驟（2）計算出來的被試能力參數 $\theta_1 $ ，估計項目參數 $ a_1, b_1 $ 。
各項目的區分度參數和難度參數 $ a_*, b_* $ 分別獨立地從建議分布 $ q_a, q_b $ 中選取，我們取區分度參數和難度參數的建議分布為 $ \log(a_*)~N(a_0, c_a^2), b*~N(b_0, c_b^2) $ ，其中 $ c_a = 0.3, c_b = 0.3 $ 。
計算從狀態 $ (a_0, b_0) $ 轉移至狀態 $ (a_*, b_*) $ 的接收概率 $ \alpha(a_0, b_0, a_*, b_*) = min(1, R_{a,b}^0) $ ，由它來決定是否發生狀態的轉移。
其中
生成隨機數 $ r_2~U(0, 1) $ ，比較 $ r_2 $ 與 $ \alpha(a_0, b_0, a_*, b_*) $ 的大小，進行狀態轉移判斷：
若 $ \alpha(a_0, b_0, a_*, b_*) \ge r_2 $ ，則 $ a_1 = a_*, b_1 = b_* $ ；否則 $ a_1 = a_0, b_1 = b_0 $ 。
4）重復步驟（2）和（3）$ n $ 次，刪除為首的 $ w $ 次，取剩余的 $ m = n - w $ 次迭代所得結果的均值即為參數的估計值。
5）步驟（4）會生成一條長度為 $ n $ 的 $ Markov $ 鏈，重復步驟（4）$ i $ 次（ $ i $ 次一般小於5），即可以得到 $ i $ 條 $ Markov $ 鏈，將這 $ i $ 條鏈得到的參數估計值的均值為最終的參數估計值。
MCMC 算法在求解3-PL模型的參數時和求2-PL模型的參數的方法一致。但是在同等的數據量下，3-PL模型的參數估計精度要低於2-PL模型的參數估計的精度（參數越多，模型越復雜，在同等量的參數下就更容易過擬合）。因此一般我們會將預測參數 $ c_j $ 標為常數值。
4、4-PL模型
在3-PL模型中引入預測參數 $ c_j $ 是為了描述被試者在沒有任何相關項目能力地情況下，答對項目的概率。然而現實中還會存在一種情況就是被試者在充分掌握了項目知識的情況下，可以認為理論上被試者關於此項目的能力是完全可以答對該項目的，但實際中卻答錯了該問題。為了解決這一問題，引入了參數 $ \gamma_j $ ，此時模型可以表示為：
$ p_{i,j}(\theta_i) = c_j + \frac {\gamma_j - c_j} {1 + \exp\,[-D\,a_j\,(\theta_i - b_j)]} $
對於4-PL模型的參數估計依然可以用估計2-PL參數的方法，此時項目參數不再只是 $ (a_j, b_j) $ ，而是 $ (a_j, b_j, c_j, \gamma_j) $ 。