bifactor graded response theory(IRT)

#EM算法，a_ini为初始输入矩阵，为简单因子分析的结果。代码算法需要Numpy包支持，基于(Gibbons等,2007,《Full-information item bifactor analysis of graded response data》,"Applied Psychological Measurement",31(1),4-19.)

# copyright @ ZhuoJun Gu，仅供研究之用，复制与传播必须遵守GNU V3.0协议 
# http://www.gnu.org/licenses/gpl-3.0.html

#EM算法 - step 0.1 参数初值计算-初始斜率
a_ini = np.minimum(a_ini,0.999)
a_ini = np.maximum(a_ini,0.001)
a_ini = a_ini * np.concatenate((np.ones(size[0])[:,np.newaxis],ujk),axis=1)
a_ini = a_ini / (1 - a_ini ** 2) ** 0.5

#EM算法 - step 0.2 参数初值计算-初始截距和类别参数（cj和dt）
c_item_ini = np.zeros((scores.shape[1],scores_trans[0].shape[1]))
for key in scores_trans:
    proportion = np.mean(scores_trans[key], axis=0)
    int_key = int(key)
    c_item_ini[int_key,:] = c_item_ini[int_key,:] + proportion
c_item_ini = c_item_ini - 10e-50
c_item_ini = abs(c_item_ini)
c_item_ini = inverse_logistic(c_item_ini) / D

#EM算法 - step 0.3 计算项目独有的参数cj
cj_item_ini = np.zeros((scores.shape[1],1))
for i in range(scores.shape[1]):
    cj_item_ini[i,0] =  np.sum(c_item_ini[i,:],axis=0) / scores_trans[0].shape[1]

#EM算法 - step 0.3 计算等级参数dt
dt_item_ini = np.zeros((scores_trans[0].shape[1],1))
for i in range(scores_trans[0].shape[1]):
    dt_item_ini[i,0] = (np.sum(c_item_ini[:,i]) - np.sum(cj_item_ini,axis=0)) / scores.shape[1]
dt_item_ini[0,0] = - np.sum(dt_item_ini[1:,0])


############EM E步和M步算法函数############

def EM_calculation(a_ini, cj_item_ini, dt_item_ini, subject_response):
    # EM算法- step 1.0 E步 更新变量r和f（或n）的期望(Gibbons,2007):

    #数据输入合法性修改:nan,inf,-inf
    a_ini = np.nan_to_num(a_ini)
    cj_item_ini = np.nan_to_num(cj_item_ini)
    dt_item_ini = np.nan_to_num(dt_item_ini) 
    
    # EM算法- step 1.1 计算rjtk,nk-张量运算

    def rjtk_nk_para_esti(subject_response):

        # 数值准备，与下面计算jacobian矩阵定义temp1不同，待类编程改进
        temp1 = np.zeros((scores.shape[1], scores_trans[0].shape[1], a_ini.shape[1] - 1, gp_size, gp_size))

        # 给theta1赋值
        temp2 = a_ini[:, 0] * x_nodes
        temp2 = temp2.transpose()
        temp1 = temp1 + temp2[:, None, None, :, None]

        # 把thetak加上去
        for i in range(a_ini.shape[1] - 1):
            temp2 = a_ini[:, i + 1] * x_nodes
            temp2 = temp2.transpose()
            temp1[:, :, i, :, :] = temp1[:, :, i, :, :] + temp2[:, None, None, :]

            # 加上截距并乘以D，然后计算Phi
        temp2 = dt_item_ini[:, 0]
        temp1 = temp1 + temp2[None, :, None, None, None]
        temp2 = cj_item_ini[:, 0]
        temp1 = D * (temp1 + temp2[:, None, None, None, None])
        temp1 = np.maximum(temp1,-70) #防止溢出
        temp1 = np.minimum(temp1,70)
        temp1 = np.exp(temp1)
        temp1 = temp1 / (1 + temp1)

        #temp1 = np.nan_to_num(temp1) + 10e-100

        # 计算likelihood
        for i in range(c_item_ini.shape[1]):
            if i == c_item_ini.shape[1] - 1:
                temp1[:, i, :, :, :] = temp1[:, i, :, :, :]
            else:
                temp1[:, i, :, :, :] = temp1[:, i, :, :, :] - temp1[:, i + 1, :, :, :]

        temp1 = np.maximum(temp1,10e-50) #不能用10e-200，会溢出
        temp1 = np.minimum(temp1,1-10e-50)

        # 根据公式里Phi外面的元素，计算pkeppa，顺便计算lik和eik
        temp2 = temp1 * subject_response[:, :, None, None, None]  # 计算选项连乘01：因为是1或0的次方，实际是相乘再求和
        temp2 = np.sum(temp2, axis=1)  # 计算选项连乘02
        temp2 = temp2 * ujk[:, :, None, None]  # 计算题项连乘01：与选项同  
        temp2 = np.sum(temp2, axis=0)  # 计算题项连乘02   
        #temp2 = np.nan_to_num(temp2) + 10e-100
        lik = temp2
        temp2 = temp2 * x_weights[:, 0][None, None, :]
        temp2 = np.sum(temp2, axis=2)
        #temp2 = np.nan_to_num(temp2) + 10e-100
        eik = temp2
        temp2 = np.prod(temp2, axis=0)
        eik = temp2 / eik
        temp2 = temp2 * x_weights[:, 0]
        temp2 = np.sum(temp2, axis=0)
        #temp2 = np.nan_to_num(temp2) + 10e-100
        pi = temp2

        return lik, eik, pi

    # 遍历每个被测获得参数
    subject_pkeppa = np.zeros((scores.shape[0]))
    subject_eik = np.zeros((scores.shape[0], a_ini.shape[1] - 1, gp_size))
    subject_lik = np.zeros((scores.shape[0], a_ini.shape[1] - 1, gp_size, gp_size))
    for i in subject_response.keys():
        subject_lik[i, :, :, :], subject_eik[i, :, :], subject_pkeppa[i] = rjtk_nk_para_esti(subject_response[i])
        if i % 10 == 0:
            print('请等待，正在计算1个循环的E步指标，当前被试编号:', i + 1)

    # 防止Pi出现负值、零值、大于1
    subject_pkeppa = np.maximum(subject_pkeppa,10e-50)
    subject_pkeppa = np.minimum(subject_pkeppa,1-10e-50)

    # 计算rjtk和nk
    temp1 = subject_lik * subject_eik[:, :, :, None] / subject_pkeppa[:, None, None, None]
    #temp1 = np.nan_to_num(temp1) + 10e-100
    if scores.shape[1] == 2:
        nk = temp1 * subject_ri[:,None,None,None]
        nk = np.sum(temp1,axis=0)
    else:
        nk = np.sum(temp1, axis=0)
    rjtk = np.zeros((scores.shape[1], a_ini.shape[1] - 1, scores_trans[0].shape[1], gp_size, gp_size))
    temp2 = np.ones((scores.shape[1], a_ini.shape[1] - 1, scores_trans[0].shape[1], gp_size, gp_size))
    for i in range(scores.shape[0]):
        if scores.shape[1] == 2:
            temp2 = temp2 * subject_response[i][:, None, :, None, None] * subject_ri[i]
            temp2 = temp2 * temp1[i, None, :, None, :, :]
            rjtk = rjtk + temp2
        else:
            temp2 = temp2 * subject_response[i][:, None, :, None, None]
            temp2 = temp2 * temp1[i, None, :, None, :, :]
            rjtk = rjtk + temp2
    #rjtk = np.nan_to_num(rjtk) + 10e-100

    # EM算法- step 2.0 M步：r和f（或n）传入，做N-S迭代（后续可考虑采用加权最小二乘，绕过矩阵计算）(Gibbons,2007):

    # EM算法- step 2.1 初始化每题雅可比矩阵、信息矩阵（海塞矩阵的期望，解决海塞病态的问题）和似然方程组，初始化各个a,c初值用于N-S迭代
    NS_dimention = 5  # 对于单个方程求解，雅可比矩阵固定为1 X 5，信息矩阵固定为5 X 5，NS意指牛顿拉夫逊迭代
    NS_dt = 2  # 由于dt和dt+1计算方式与ac不同，必须单独计算，并用np.concatenate合并
    NS_ac = 3  # a1,ak,cj的矩阵大小参数

    # 2.1a 计算Jacobian-ac部分

    # Jacobian准备
    temp1 = np.zeros((scores.shape[1], NS_ac, a_ini.shape[1] - 1, scores_trans[0].shape[1], gp_size, gp_size))
    temp2 = a_ini[:, 0] * x_nodes
    temp2 = temp2.transpose()
    temp1 = temp1 + temp2[:, None, None, None, :, None]
    for i in range(a_ini.shape[1] - 1):
        temp2 = a_ini[:, i + 1] * x_nodes
        temp2 = temp2.transpose()
        temp1[:, :, i, :, :, :] = temp1[:, :, i, :, :, :] + temp2[:, None, None, None, :]
    temp2 = dt_item_ini[:, 0]
    temp1 = temp1 + temp2[None, None, None, :, None, None]
    temp2 = cj_item_ini[:, 0]
    temp1 = D * (temp1 + temp2[:, None, None, None, None, None])
    temp1 = np.maximum(temp1,-25) #防止溢出
    temp1 = np.minimum(temp1,25)
    temp1 = np.exp(temp1)
    temp1 = temp1 / (1 + temp1)
    #temp1 = np.nan_to_num(temp1) + 10e-100
    pjt = temp1[:, 0, :, :, :, :]

    # 求Pjt
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            pjt[:, :, i, :, :] = pjt[:, :, i, :, :]
        else:
            pjt[:, :, i, :, :] = pjt[:, :, i, :, :] - pjt[:, :, i + 1, :, :]
            # 对Pjt做概率修正，因为Pk-Pk+1结果有0值或负值不能被P/v除。具体原因未知。后续进一步做数值计算研究。此情况前面也有。
    pjt = np.maximum(pjt,10e-50)
    pjt = np.minimum(pjt,1-10e-50)

    # 求一阶偏导P/v
    temp2 = np.zeros((NS_ac, a_ini.shape[1] - 1, gp_size, gp_size))
    temp3 = x_nodes[:, 0]
    for i in range(a_ini.shape[1] - 1):
        # 构造每道题的Jacobian矩阵
        temp2[0, i, :, :] = (temp2[0, i, :, :] + temp3[:, None]) * D  # J矩阵a1偏导：theta 1 * D
        temp2[1, i, :, :] = (temp2[1, i, :, :] + temp3[None, :]) * D  # J矩阵ak偏导：theta k * D
        temp2[2, i, :, :] = temp2[2, i, :, :] + D  # J矩阵cj的偏导：D
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            temp1[:, :, :, i, :, :] = temp1[:, :, :, i, :, :] * temp2[None, :, :, :, :]
        else:
            temp1[:, :, :, i, :, :] = temp1[:, :, :, i, :, :] * temp2[None, :, :, :, :] - temp1[:, :, :, i + 1, :,
                                                                                          :] * temp2[None, :, :, :,
                                                                                               :]  # 这里加截距的张量阶数与E步不同，因为jacobian矩阵定义所致，后期类编程注意调整

    # 一阶偏导P/v除以Pjt,并计算jacobian
    temp1 = temp1 / pjt[:, None, :, :, :, :]
    temp1 = temp1 * rjtk[:, None, :, :, :, :]
    temp2 = x_weights[:, 0]
    temp1 = temp1 * temp2[None, None, None, None, None, :]
    temp1 = temp1 * temp2[None, None, None, None, :, None]
    temp1 = np.sum(temp1, axis=3)
    temp1 = np.sum(temp1, axis=4)
    temp1 = temp1 * ujk[:,None,:, None]
    temp1 = np.sum(temp1, axis=2)
    jacobian = np.sum(temp1, axis=2)
    temp1 = np.zeros((scores.shape[1], scores_trans[0].shape[1], NS_ac))
    temp1 = temp1 + jacobian[:,None,:]
    #temp1 = np.nan_to_num(temp1) + 10e-100
    jacobian = temp1

    # 2.1b 计算Jacobian-dt部分

    # Jacobian准备
    temp1 = np.zeros((scores.shape[1], scores_trans[0].shape[1], NS_dt, a_ini.shape[1] - 1, gp_size, gp_size))
    temp2 = a_ini[:, 0] * x_nodes
    temp2 = temp2.transpose()
    temp1 = temp1 + temp2[:, None, None, None, :, None]
    for i in range(a_ini.shape[1] - 1):
        temp2 = a_ini[:, i + 1] * x_nodes
        temp2 = temp2.transpose()
        temp1[:, :, :, i, :, :] = temp1[:, :, :, i, :, :] + temp2[:, None, None, None, :]
    temp2 = dt_item_ini[:, 0]
    temp1 = temp1 + temp2[None, :, None, None, None, None]
    temp2 = cj_item_ini[:, 0]
    temp1 = D * (temp1 + temp2[:, None, None, None, None, None])
    temp1 = np.maximum(temp1,-25) #防止溢出
    temp1 = np.minimum(temp1,25) 
    temp1 = np.exp(temp1)
    temp1 = temp1 / (1 + temp1)
    #temp1 = np.nan_to_num(temp1) + 10e-100
    pjt = temp1[:, :, 0, :, :, :]

    # 求Pjt
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            pjt[:, i, :, :, :] = pjt[:, i, :, :, :]
        else:
            pjt[:, i, :, :, :] = pjt[:, i, :, :, :] - pjt[:, i + 1, :, :, :]
            # 对Pjt做概率修正，因为Pk-Pk+1结果有0值或负值不能被P/v除。具体原因未知。后续进一步做数值计算研究。此情况前面也有。
    pjt = np.maximum(pjt,10e-50)
    pjt = np.minimum(pjt,1-10e-50)

    # 求一阶偏导P/v
    temp2 = np.zeros((NS_dt, a_ini.shape[1] - 1, gp_size, gp_size))
    temp3 = x_nodes[:, 0]
    for i in range(a_ini.shape[1] - 1):
        # 构造每道题的Jacobian矩阵
        temp2[0, i, :, :] = temp2[0, i, :, :] + D  # J矩阵dt偏导：D
        temp2[1, i, :, :] = temp2[1, i, :, :] + D  # J矩阵dt+1偏导：D
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            temp1[:, i, :, :, :, :] = temp1[:, i, :, :, :, :] * temp2[None, :, :, :, :]
        else:
            temp1[:, i, :, :, :, :] = temp1[:, i, :, :, :, :] * temp2[None, :, :, :, :] - temp1[:, i + 1, :, :, :,
                                                                                          :] * temp2[None, :, :, :,
                                                                                               :]  # 这里加截距的张量阶数与E步不同，因为jacobian矩阵定义所致，后期类编程注意调整

    # 一阶偏导P/v除以Pjt乘以rjkt,并计算jacobian
    temp2 = np.zeros((scores.shape[1], c_item_ini.shape[1], a_ini.shape[1] - 1, gp_size, gp_size))
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            temp2[:, i, :, :, :] = rjtk[:, :, i, :, :] / pjt[:, i, :, :, :]
        else:
            temp2[:, i, :, :, :] = rjtk[:, :, i, :, :] / pjt[:, i, :, :, :] - rjtk[:, :, i + 1, :, :] / pjt[:, i + 1, :,
                                                                                                        :, :]
    temp1 = temp1 * temp2[:, :, None, :, :, :]
    temp2 = x_weights[:, 0]
    temp1 = temp1 * temp2[None, None, None, None, None, :]
    temp1 = temp1 * temp2[None, None, None, None, :, None]
    temp1 = np.sum(temp1, axis=0)
    temp1 = np.sum(temp1, axis=4)
    temp2 = np.zeros((ujk.shape[0], temp1.shape[0], temp1.shape[1], temp1.shape[2], temp1.shape[3]))
    temp2 = temp2 + ujk[:, None, None, :, None]
    temp2 = temp2 * temp1[None, :, :, :, :]
    temp1 = temp2
    temp1 = np.sum(temp1, axis=3)
    temp1 = np.sum(temp1, axis=3)
    #temp1 = np.nan_to_num(temp1) + 10e-100

    # 合并
    jacobian = np.concatenate((jacobian, temp1), axis=2)

    # 2.1c 计算Information_EHessian-参数ac部分

    # information准备
    temp1 = np.zeros((scores.shape[1], NS_ac, NS_ac, scores_trans[0].shape[1], a_ini.shape[1] - 1, gp_size, gp_size))
    temp2 = a_ini[:, 0] * x_nodes
    temp2 = temp2.transpose()
    temp1 = temp1 + temp2[:, None, None, None, None, :, None]
    for i in range(a_ini.shape[1] - 1):
        temp2 = a_ini[:, i + 1] * x_nodes
        temp2 = temp2.transpose()
        temp1[:, :, :, :, i, :, :] = temp1[:, :, :, :, i, :, :] + temp2[:, None, None, None, None, :]
    temp2 = dt_item_ini[:, 0]
    temp1 = temp1 + temp2[None, None, None, :, None, None, None]
    temp2 = cj_item_ini[:, 0]
    temp1 = D * (temp1 + temp2[:, None, None, None, None, None, None])
    temp1 = np.maximum(temp1,-25) #防止溢出
    temp1 = np.minimum(temp1,25)
    temp1 = np.exp(temp1)
    temp1 = temp1 / (1 + temp1)
    #temp1 = np.nan_to_num(temp1) + 10e-100

    # 求一阶偏导P/v
    temp2 = np.zeros((scores.shape[1], NS_ac, NS_ac, a_ini.shape[1] - 1, gp_size, gp_size))
    temp3 = x_nodes[:, 0]
    for i in range(a_ini.shape[1] - 1):
        # 构造每道题目的信息矩阵
        temp2[:, 0, 0, i, :, :] = (temp2[:, 0, 0, i, :, :] + temp3[None, :, None]) ** 2 * D ** 2  # 信息矩阵：theta 1^2 * D^2
        temp2[:, 0, 1, i, :, :] = (temp2[:, 0, 1, i, :, :] + temp3[None, :, None]) * (
                    temp2[:, 0, 1, i, :, :] + temp3[None, None, :]) * D ** 2  # 信息矩阵：theta 1 * theta k * D^2
        temp2[:, 0, 2, i, :, :] = (temp2[:, 0, 2, i, :, :] + temp3[None, :, None]) * D ** 2  # 信息矩阵：theta 1 * D^2
        temp2[:, 1, 0, i, :, :] = temp2[:, 0, 1, i, :, :]
        temp2[:, 1, 1, i, :, :] = (temp2[:, 0, 0, i, :, :] + temp3[None, None, :]) ** 2 * D ** 2  # 信息矩阵：theta k^2 * D^2
        temp2[:, 1, 2, i, :, :] = (temp2[:, 0, 1, i, :, :] + temp3[None, None, :]) * D ** 2  # 信息矩阵：theta k * D^2
        temp2[:, 2, 0, i, :, :] = temp2[:, 0, 2, i, :, :]
        temp2[:, 2, 1, i, :, :] = temp2[:, 1, 2, i, :, :]
        temp2[:, 2, 2, i, :, :] = D ** 2  # 信息矩阵：D^2
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            temp1[:, :, :, i, :, :, :] = temp1[:, :, :, i, :, :, :] ** 2 * temp2[:, :, :, :, :, :]
        else:
            temp1[:, :, :, i, :, :, :] = (temp1[:, :, :, i, :, :, :] - temp1[:, :, :, i + 1, :, :, :]) ** 2 * temp2[:,
                                                                                                              :, :, :,
                                                                                                              :, :]

    # 一阶偏导P/v除以Pjt,并计算information
    temp1 = temp1 / pjt[:, None, None, :, :, :, :]
    temp1 = np.sum(temp1, axis=3)
    temp1 = temp1 * nk[None, None, None, :, :, :]
    temp1 = np.sum(temp1, axis=5)
    temp1 = np.sum(temp1, axis=4)
    #temp1 = np.nan_to_num(temp1) + 10e-100
    information = np.zeros((scores.shape[1], scores_trans[0].shape[1], NS_ac, NS_ac, a_ini.shape[1] - 1))
    for i in range(c_item_ini.shape[1]):
        information[:, i, :, :, :] = information[:, i, :, :, :] + temp1[:, :, :, :]

    # 2.1d 计算Information_EHessian-参数dt部分

    # information准备
    # ac_dt偏导部分
    temp1 = np.zeros((scores.shape[1], NS_dt, NS_ac, scores_trans[0].shape[1], a_ini.shape[1] - 1, gp_size, gp_size))
    temp2 = a_ini[:, 0] * x_nodes
    temp2 = temp2.transpose()
    temp1 = temp1 + temp2[:, None, None, None, None, :, None]
    for i in range(a_ini.shape[1] - 1):
        temp2 = a_ini[:, i + 1] * x_nodes
        temp2 = temp2.transpose()
        temp1[:, :, :, :, i, :, :] = temp1[:, :, :, :, i, :, :] + temp2[:, None, None, None, None, :]
    temp2 = dt_item_ini[:, 0]
    temp1 = temp1 + temp2[None, None, None, :, None, None, None]
    temp2 = cj_item_ini[:, 0]
    temp1 = D * (temp1 + temp2[:, None, None, None, None, None, None])
    temp1 = np.maximum(temp1,-25) #防止溢出
    temp1 = np.minimum(temp1,25)
    temp1 = np.exp(temp1)
    temp1 = temp1 / (1 + temp1)
    #temp1 = np.nan_to_num(temp1) + 10e-100

    # 求一阶偏导P/d
    p_d = temp1[:, 0, 0, :, :, :, :]
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            p_d[:, i, :, :, :] = p_d[:, i, :, :, :] * D
        else:
            p_d[:, i, :, :, :] = (p_d[:, i, :, :, :] - p_d[:, i + 1, :, :, :]) * D
    #p_d = np.nan_to_num(p_d) + 10e-100

    # 求一阶偏导P/v
    temp2 = np.zeros((scores.shape[1], NS_dt, NS_ac, a_ini.shape[1] - 1, gp_size, gp_size))
    temp3 = x_nodes[:, 0]
    for i in range(a_ini.shape[1] - 1):
        temp2[:, 0, 0, i, :, :] = (temp2[:, 0, 0, i, :, :] + temp3[None, :, None]) * D
        temp2[:, 0, 1, i, :, :] = (temp2[:, 0, 1, i, :, :] + temp3[None, None, :]) * D
        temp2[:, 0, 2, i, :, :] = temp2[:, 0, 2, i, :, :] + D ** 2
        temp2[:, 1, 0, i, :, :] = temp2[:, 0, 0, i, :, :]
        temp2[:, 1, 1, i, :, :] = temp2[:, 0, 1, i, :, :]
        temp2[:, 1, 2, i, :, :] = temp2[:, 0, 2, i, :, :]
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            temp1[:, :, :, i, :, :, :] = temp1[:, :, :, i, :, :, :] * temp2[:, :, :, :, :, :]
        else:
            temp1[:, :, :, i, :, :, :] = temp1[:, :, :, i, :, :, :] * temp2[:, :, :, :, :, :] - temp1[:, :, :, i + 1, :,
                                                                                                :, :] * temp2[:, :, :,
                                                                                                        :, :, :]

    # 一阶偏导P/v除以Pjt,并计算information
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            temp1[:, :, :, i, :, :, :] = temp1[:, :, :, i, :, :, :] / pjt[:, None, None, i, :, :, :]
        else:
            temp1[:, :, :, i, :, :, :] = temp1[:, :, :, i, :, :, :] / pjt[:, None, None, i, :, :, :] - temp1[:, :, :,
                                                                                                       i + 1, :, :,
                                                                                                       :] / pjt[:, None,
                                                                                                            None, i + 1,
                                                                                                            :, :, :]
    temp1 = temp1 * p_d[:, None, None, :, :, :, :]
    temp1 = np.sum(temp1, axis=0)
    temp1 = temp1 * nk[None, None, None, :, :, :]
    temp1 = np.sum(temp1, axis=5)
    temp1 = - np.sum(temp1, axis=4)
    #temp1 = np.nan_to_num(temp1) + 10e-100
    information_ac_dt = np.zeros((scores.shape[1], scores_trans[0].shape[1], NS_dt, NS_ac, a_ini.shape[1] - 1))
    for i in range(c_item_ini.shape[1]):
        information_ac_dt[:, i, :, :, :] = information_ac_dt[:, i, :, :, :] + temp1[None, :, :, i, :]

    # dt偏导部分
    temp1 = np.zeros((scores.shape[1], NS_dt, NS_dt, scores_trans[0].shape[1], a_ini.shape[1] - 1, gp_size, gp_size))
    for i in range(NS_dt):
        temp1[:, i, i, :, :, :, :] = temp1[:, i, i, :, :, :, :] + p_d[:, :, :, :, :] ** 2
        for j in range(c_item_ini.shape[1]):
            if j == c_item_ini.shape[1] - 1:
                temp1[:, i, i, j, :, :, :] = temp1[:, i, i, j, :, :, :] / pjt[:, j, :, :, :]
            else:
                temp1[:, i, i, j, :, :, :] = temp1[:, i, i, j, :, :, :] * (
                            1 / pjt[:, j, :, :, :] + 1 / pjt[:, j + 1, :, :, :])
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            temp1[:, 0, 1, i, :, :, :] = temp1[:, 0, 1, i, :, :, :] + p_d[:, i, :, :, :]
        else:
            temp1[:, 0, 1, i, :, :, :] = temp1[:, 0, 1, i, :, :, :] + p_d[:, i, :, :, :] * p_d[:, i + 1, :, :, :]
    temp1[:, 1, 0, :, :, :, :] = temp1[:, 0, 1, :, :, :, :]
    temp1 = np.sum(temp1, axis=0)
    temp1 = temp1 * nk[None, None, None, :, :, :]
    temp1 = np.sum(temp1, axis=5)
    temp1 = - np.sum(temp1, axis=4)
    #temp1 = np.nan_to_num(temp1) + 10e-100
    information_dt_dt = np.zeros((scores.shape[1], scores_trans[0].shape[1], NS_dt, NS_dt, a_ini.shape[1] - 1))
    for i in range(c_item_ini.shape[1]):
        information_dt_dt[:, i, :, :, :] = information_dt_dt[:, i, :, :, :] + temp1[None, :, :, i, :]

    # 合并
    temp1 = information
    temp1 = np.concatenate((temp1, np.transpose(information_ac_dt, axes=[0, 1, 3, 2, 4])), axis=3)
    temp2 = np.concatenate((information_ac_dt, information_dt_dt), axis=3)
    information = np.concatenate((temp1, temp2), axis=2)

    # EM算法- step 2.2 牛顿-拉夫逊迭代（N-S迭代，因计算效率，仅迭代1次）

    parameter = np.zeros((scores.shape[1], scores_trans[0].shape[1], NS_dimention, a_ini.shape[1] - 1))
    for i in range(a_ini.shape[1] - 1):
        parameter[:, :, :, i] = np.linalg.solve(information[:, :, :, :, i], jacobian)
    parameter = parameter * ujk[:, None, None, :]  # 对有变化的无关参数置0
    #parameter = np.nan_to_num(parameter) + 10e-100

    # 取a1
    temp1 = np.sum(parameter, axis=3)
    temp1 = np.mean(temp1, axis=1)
    temp1 = temp1[:, 0]
    #temp1 = np.nan_to_num(temp1) + 10e-100
    a_ini[:, 0] = a_ini[:, 0] - temp1

    # 取ak
    for i in range(a_ini.shape[1] - 1):
        temp1 = parameter[:, :, :, i]
        temp1 = np.mean(temp1, axis=1)
        temp1 = temp1[:, 1]
        #temp1 = np.nan_to_num(temp1) + 10e-100
        a_ini[:, i + 1] = a_ini[:, i + 1] - temp1

    #防止区分度过大，但与事实不符
    #a_ini = np.minimum(a_ini,100) 
    #a_ini = np.maximum(a_ini,-100)
    #a_ini = a_ini * np.concatenate((np.ones(size[0])[:,np.newaxis],ujk),axis=1)

    # 取cj
    temp1 = np.sum(parameter, axis=3)
    temp1 = np.mean(temp1, axis=1)
    temp1 = temp1[:, 2]
    #temp1 = np.nan_to_num(temp1) + 10e-100
    cj_item_ini[:,0] = cj_item_ini[:,0] - temp1

    # 取dt
    temp1 = np.sum(parameter, axis=3)
    temp1 = np.mean(temp1, axis=0)
    temp1 = temp1[:, 3]
    #temp1 = np.nan_to_num(temp1) + 10e-100
    dt_item_ini[:,0] = dt_item_ini[:,0] - temp1
    dt_item_ini[0,0] = - np.sum(dt_item_ini[1:,0])

    #数据输出合法性修改:nan,inf,-inf
    a_ini = np.nan_to_num(a_ini)
    cj_item_ini = np.nan_to_num(cj_item_ini)
    dt_item_ini = np.nan_to_num(dt_item_ini) 

    return a_ini, cj_item_ini, dt_item_ini
    

##########EM算法 - whole step： 初始化各种参数开始迭代，或与上一次似然值比较，并评估迭代次数，并更新参数以继续迭代

def pi(subject_response):
    # 数值准备，与下面计算jacobian矩阵定义temp1不同，待类编程改进
    temp1 = np.zeros((scores.shape[1], scores_trans[0].shape[1], a_ini.shape[1] - 1, gp_size, gp_size))

    # 给theta1赋值
    temp2 = a_ini[:, 0] * x_nodes
    temp2 = temp2.transpose()
    temp1 = temp1 + temp2[:, None, None, :, None]

    # 把thetak加上去
    for i in range(a_ini.shape[1] - 1):
        temp2 = a_ini[:, i + 1] * x_nodes
        temp2 = temp2.transpose()
        temp1[:, :, i, :, :] = temp1[:, :, i, :, :] + temp2[:, None, None, :]

        # 加上截距并乘以D，然后计算Phi
    temp2 = dt_item_ini[:, 0]
    temp1 = temp1 + temp2[None, :, None, None, None]
    temp2 = cj_item_ini[:, 0]
    temp1 = D * (temp1 + temp2[:, None, None, None, None])
    temp1 = np.maximum(temp1,-25) #防止溢出
    temp1 = np.minimum(temp1,25)
    temp1 = np.exp(temp1)
    temp1 = temp1 / (1 + temp1)
    #temp1 = np.nan_to_num(temp1) + 10e-100

    # 计算likelihood
    for i in range(c_item_ini.shape[1]):
        if i == c_item_ini.shape[1] - 1:
            temp1[:, i, :, :, :] = temp1[:, i, :, :, :]
        else:
            temp1[:, i, :, :, :] = temp1[:, i, :, :, :] - temp1[:, i + 1, :, :, :]
    temp1 = np.maximum(temp1,10e-50)
    temp1 = np.minimum(temp1,1-10e-50)

    # 根据公式里Phi外面的元素计算pkeppa
    temp2 = temp1 * subject_response[:, :, None, None, None]  # 计算选项连乘01：因为是1或0的次方，实际是相乘再求和
    temp2 = np.sum(temp2, axis=1)  # 计算选项连乘02
    temp2 = temp2 * ujk[:, :, None, None]  # 计算题项连乘01：与选项同  
    temp2 = np.sum(temp2, axis=0)  # 计算题项连乘02   
    temp2 = temp2 * x_weights[:, 0][None, None, :]
    temp2 = np.sum(temp2, axis=2)
    temp2 = np.prod(temp2, axis=0)
    temp2 = temp2 * x_weights[:, 0]
    temp2 = np.sum(temp2, axis=0)
    pi = np.log(temp2)
    pi = np.maximum(pi,10e-50)
    pi = np.minimum(pi,1-10e-50)
    return pi


#获取第i位被测的回答模式
subject_response = {}
for i in range(scores.shape[0]):
    subject_response[i] = np.zeros((c_item_ini.shape[0],c_item_ini.shape[1]))
    for j in scores_trans.keys():
        subject_response[i][j,:] = scores_trans[j][i,:]

iteration = 100
change = 10e-5
log_likelihood_previous = np.zeros((scores.shape[0],1))
log_likelihood_next = np.zeros((scores.shape[0],1))

for i in range(iteration):
    for j in range(scores.shape[0]):
        if scores.shape[1] == 2:
            log_likelihood_previous = pi(subject_response[j]) * subject_ri[j]
        else:
            log_likelihood_previous = pi(subject_response[j])
    a_ini,cj_item_ini,dt_item_ini = EM_calculation(a_ini,cj_item_ini,dt_item_ini,subject_response)
    for j in range(scores.shape[0]):
        if scores.shape[1] == 2:
            log_likelihood_next = pi(subject_response[j]) * subject_ri[j]
        else:
            log_likelihood_next = pi(subject_response[j])
    change_temp = abs(np.sum(log_likelihood_previous,axis=0) - np.sum(log_likelihood_next,axis=0))
    print('请等待，这是EM算法第',i+1,'次循环，最多',iteration,'次循环')
    print('当前残差:',change_temp,'，与EM算法停止残差标准：',change)
    if change_temp < change:
        break

#取因素负荷
alpha_factorloading = a_ini + 1e-6
alpha_factorloading = ((alpha_factorloading ** -4 + 4) ** 0.5 - alpha_factorloading ** -2) / 2
alpha_factorloading = alpha_factorloading * np.concatenate((np.ones(size[0])[:,np.newaxis],ujk),axis=1)
temp1 = np.var(scores)
alpha_factorloading = alpha_factorloading / temp1


############标准化后做斜交旋转############
#alpha_factorloading = np.maximum(alpha_factorloading,0.001)
#alpha_factorloading = np.minimum(alpha_factorloading,0.999)
#alpha_factorloading = alpha_factorloading * np.concatenate((np.ones(size[0])[:,np.newaxis],ujk),axis=1)
#alpha_factorloading = np.round(alpha_factorloading,decimals=3)

np.savetxt('d:/alpha.csv',alpha_factorloading,delimiter=',',newline='\n')

###end估计参数###