python概率密度函数参数估计_Python与项目反应理论：基于EM和MCMC的参数估计算法...

项目反应理论的开端

早在上世纪初，智力测验的发明者比奈(也可能是西蒙)便发现了一条神奇的曲线，这条曲线的x轴是智力水平，y轴是试题正确率，而这是项目反应理论(以下简称IRT)的最初雏形。上世界五六十年代，ETS的统计学家Lord经过一系列的工作，正式开创了IRT理论。

为什么在经典测验理论(以下简称CTT)存在的情况下，还要继续IRT理论呢？因为CTT的统计模型毛病超级多。首先，CTT的核心概念——信度，无法计算，这是最致命的；其次，CTT统计模型过于简单，像猜测度、失误、时间等等参数都很难纳入模型。所以抛弃了信度概念的IRT应运而生。同时IRT可以解决CTT完全应付不了的问题，比如组卷，IRT可以很好的控制误差和等值，CTT想正儿八经的控制误差和等值难于登天，比如理想点反应机制，IRT有理想点模型处理，CTT无法处理，比如排序数据(CTT假设真分数和误差独立，但在排序数据中，真分数和误差明显不是独立的，而IRT的瑟斯顿模型可以解决这个问题)。

IRT的函数

对于最初的IRT来说，有三条基本假设，一是单维性，二是局部独立性，三是项目反应函数假设，但其实前两条假设可有可无(例如mirt理论打破了单维性假设，题组反应理论打破了局部独立性假设)，核心的是第三个假设。

对于单维度IRT来说，最常用的两个函数，一个是正态肩形曲线函数，一个是logistic函数(跟茴香豆的N种叫法一样，你也可以叫它sigmoid函数)。我们这里主要说logistic函数。logisitc的推导没什么好说的，无非是假设试题作答正确概率为

， 则答错的概率为

，取自然对数比，则

，而

是潜在特质

的线性函数

，于是我们可以得到函数

(这个推导其实是搞笑的，正儿八经的推导是基于正态肩曲线的)，其中

可以称为区分度(或斜率)，

可以称为阈值(或截距)，而

则是难度(或通俗度)。注意，有时候这个logistic函数上会有个特殊参数D，例如

，D通常等于1.702，这是为了让logistic函数更接近正态肩形曲线，可以证明

，但是加D不加D意义不大。

上面的项目反应函数中，如果区分度

恒等于1，那么这个函数也称为Rasch函数，是由丹麦数学家Rasch独立提出的一种统计模型，其中

呈现随机效应，而

呈现固定效应，学过线性模型的人能看出这是一个混合模型。

双参数二级计分模型的参数估计

我们先对较为简单的双参数二级计分模型进行参数估计，模型即为之前的

。

先定义一个基础类

from __future__ import print_function, division

import numpy as np

import warnings

class BaseIrt(object):

def __init__(self, scores=None):

self.scores = scores

@staticmethod

def p(z):

# 回答正确的概率函数

e = np.exp(z)

p = e / (1.0 + e)

return p

然后定义一个Irt2PL继承这个类，并加上了计算z函数的静态方法，对于z函数的值，我们加了一些限制，这是为了防止数值溢出。

class Irt2PL(BaseIrt):

''

@staticmethod

def z(slop, threshold, theta):

# z函数

_z = slop * theta + threshold

_z[_z > 35] = 35

_z[_z < -35] = -35

return _z

由于

和

均为未知，我们可以采用EM算法(当然，我们也可以用MCMC算法)，把

当缺失数据。E步，计算

下的样本量分布(人数)和答对试题的样本量分布(人数)，M步，极大似然法求解

和

的值 。由此，Irt2PL类的构造函数如下

class Irt2PL(BaseIrt):

# EM算法求解

def __init__(self, init_slop=None, init_threshold=None, max_iter=10000, tol=1e-5, gp_size=11,

m_step_method='newton', *args, **kwargs):

""":param init_slop: 斜率初值:param init_threshold: 阈值初值:param max_iter: EM算法最大迭代次数:param tol: 精度:param gp_size: Gauss–Hermite积分点数"""

super(Irt2PL, self).__init__(*args, **kwargs)

# 斜率初值

if init_slop is not None:

self._init_slop = init_slop

else:

self._init_slop = np.ones(self.scores.shape[1])

# 阈值初值

if init_threshold is not None:

self._init_threshold = init_threshold

else:

self._init_threshold = np.zeros(self.scores.shape[1])

self._max_iter = max_iter

self._tol = tol

self._m_step_method = '_{0}'.format(m_step_method)

self.x_nodes, self.x_weights = self.get_gh_point(gp_size)

E步的求解

从直觉的角度，依据贝叶斯法则，

，所以

下样本量分布(人数分布)为

，其中

是概率密度函数，

是试题作答模式

下的似然函数，而

下答对试题的样本量分布(人数分布)为

，

代表的是第

个作答模式下第

道题的答题情况。很明显，前面的公式为连续变量，很难求解，需要用数值积分的方法，考虑到

通常假设为正态分布的概率密度函数，所以我们可以用最简单的Gauss–Hermite积分求解(当然，最好的方法是自适应积分)。

Gauss–Hermite积分

Gauss–Hermite积分形式如下

我们假设

，于是我们的Gauss–Hermite积分形式为

，我们在Irt2PL类上添加一个静态方法，处理Gauss–Hermite积分

class Irt2PL(BaseIrt):

''

@staticmethod