update readme

demo/demo_irt.py

@@ -3,11 +3,17 @@
 from __future__ import print_function, division, unicode_literals
 from psy import Irt, data
 
+# score = data['lsat.dat']
+# model = Irt(scores=score, link='logit', params_type='1PL')
+# res = model.fit()
+# print(res)
+#
+# model = Irt(scores=score, link='logit')
+# res = model.fit()
+# print(res)
+
 score = data['lsat.dat']
-model = Irt(scores=score, link='probit')
+model = Irt(scores=score, link='logit', params_type='3PL')
 res = model.fit()
 print(res)
 
-model = Irt(scores=score, link='logit')
-res = model.fit()
-print(res)

demo/demo_mirt.py

@@ -4,6 +4,6 @@
 from psy import Mirt, data
 
 score = data['lsat.dat']
-res = Mirt(scores=score, dim_size=2).em()
+res = Mirt(scores=score, dim_size=2).fit()
 print(res[2])
 print(res[0])

docs/conf.py

@@ -84,7 +84,10 @@
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
-html_theme = 'alabaster'
+# html_theme = 'alabaster'
+import sphinx_rtd_theme
+html_theme = "sphinx_rtd_theme"
+html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
 # Theme options are theme-specific and customize the look and feel of a
 # theme further.  For a list of options available for each theme, see the

psy/__init__.py

@@ -5,7 +5,8 @@
 __version__ = '0.0.1'
 
 from psy.cdm.irm import McmcHoDina, McmcDina, EmDina, MlDina
-from psy.irt.irm import Mirt, Irt
+from psy.irt.irm import Irt
+from psy.irt.mirm import Mirt
 from psy.irt.grm import Grm
 from psy.cat.tirt import SimAdaptiveTirt
 from psy.fa.rotations import GPForth

psy/irt/base.py

@@ -3,34 +3,44 @@
 import numpy as np
 
 
-class ZMixin(object):
+class _ZMixin(object):
 
-    def z(self, slop, threshold, theta):
+    def z(self, theta, threshold, slop=None):
         """
-
+        z function
         :param slop:
         :param threshold:
         :param theta:
         :return:
         """
+        if slop is None:
+            return theta + threshold
         return slop * theta + threshold
 
 
-class RaschZMixin(object):
+class _MirtZMixin(object):
 
-    def z(self, threshold, theta):
-        # z函数
-        return theta + threshold
+    def z(self, slop, threshold, theta):
+        z_val = np.dot(theta, slop) + threshold
+        z_val[z_val > 35] = 35
+        z_val[z_val < -35] = -35
+        return z_val
 
 
-class ProbitMixin(object):
+class _ProbitMixin(object):
+    """
+    probit模型
+    """
 
     def p(self, z):
         # 回答正确的概率函数
         return norm.cdf(z)
 
 
-class LogitMixin(object):
+class _LogitMixin(object):
+    """
+    logit模型
+    """
 
     def p(self, z):
         # 回答正确的概率函数
@@ -39,16 +49,75 @@ def p(self, z):
         return p_val
 
 
-class GuessLogitMixin(LogitMixin):
-
-    def p(self, guess, *args, **kwargs):
-        p_val = super(GuessLogitMixin, self).p(*args, **kwargs)
-        return guess + (1 - guess) * p_val
-
-
-class GuessProbitMixin(ProbitMixin):
-
-    def p(self, guess, *args, **kwargs):
-        p_val = super(GuessProbitMixin, self).p(*args, **kwargs)
-        return guess + (1 - guess) * p_val
-
+class _GuessMixin(object):
+    """
+    带有猜测参数的模型
+    """
+    def p_guess(self, guess, *args, **kwargs):
+        z_val = self.z(*args, **kwargs)
+        p_val = self.p(z_val)
+        return {'p_guess_val': guess + (1 - guess) * p_val, 'p_val': p_val}
+
+
+class LogitMixin(_LogitMixin, _ZMixin):
+    """
+    1参数和2参数logit模型
+    """
+
+
+class ProbitMixin(_ProbitMixin, _ZMixin):
+    """
+    1参数和2参数probit模型
+    """
+
+
+class Logit3PLMixin(_LogitMixin, _GuessMixin, _ZMixin):
+    """
+    3参数logit模型mixin
+    """
+
+
+class Probit3PLMixin(_ProbitMixin, _GuessMixin, _ZMixin):
+    """
+    3参数probit模型mixin
+    """
+
+
+class MirtLogit2PLMixin(_LogitMixin, _MirtZMixin):
+    """
+    2参数多维logit模型
+    """
+
+
+class BaseEmIrt(object):
+
+    def __init__(self, scores=None, max_iter=1000, tol=1e-5):
+        self.scores = scores
+        self._max_iter = max_iter
+        self._tol = tol
+        self.item_size = scores.shape[1]
+
+    def _e_step(self, p_val, weights):
+        # 计算theta的分布人数
+        scores = self.scores
+        lik_wt = self._lik(p_val) * weights
+        # 归一化
+        lik_wt_sum = np.sum(lik_wt, axis=0)
+        _temp = lik_wt / lik_wt_sum
+        # theta的人数分布
+        full_dis = np.sum(_temp, axis=1)
+        # theta下回答正确的人数分布
+        right_dis = np.dot(_temp, scores)
+        full_dis.shape = full_dis.shape[0], 1
+        # 对数似然值
+        loglik_val = np.sum(np.log(lik_wt_sum))
+        return full_dis, right_dis, loglik_val
+
+    def _lik(self, p_val):
+        # 似然函数
+        scores = self.scores
+        p_val[p_val <= 0] = 1e-10
+        p_val[p_val >= 1] = 1 - 1e-10
+        loglik_val = np.dot(np.log(p_val), scores.transpose()) + \
+                     np.dot(np.log(1 - p_val), (1 - scores).transpose())
+        return np.exp(loglik_val)

psy/irt/grm.py

@@ -78,7 +78,7 @@ def _lik(self, p_val_dt):
             for j in range(rep_len):
                 p_pre = 1 if j == 0 else p_val_dt[i][:, j - 1]
                 p = 0 if j == rep_len - 1 else p_val_dt[i][:, j]
-                loglik_val += np.dot(np.log(p_pre - p + 1e-200)[:, np.newaxis], scores[i][:,  j][np.newaxis])
+                loglik_val += np.dot(np.log(p_pre - p + 1e-10)[:, np.newaxis], scores[i][:,  j][np.newaxis])
         return np.exp(loglik_val)
 
     def _e_step(self, p_val_dt, weights):
@@ -110,10 +110,10 @@ def _item_jac(p_val, pq_val, right_dis, len_threshold, rep_len, theta):
             temp1 = _theta * right_dis[:, i] * (1 - p_pre - p)
             dloglik_val[-1] += np.sum(temp1)
             if i < rep_len - 1:
-                temp2 = right_dis[:, i] * pq / (p - p_pre + 1e-200)
+                temp2 = right_dis[:, i] * pq / (p - p_pre + 1e-10)
                 dloglik_val[i] += np.sum(temp2)
             if i > 0:
-                temp3 = right_dis[:, i] * pq_pre / (p_pre - p + 1e-200)
+                temp3 = right_dis[:, i] * pq_pre / (p_pre - p + 1e-10)
                 dloglik_val[i - 1] += np.sum(temp3)
         return dloglik_val
 
@@ -126,17 +126,17 @@ def _item_hess(p_val, pq_val, full_dis, len_threshold, rep_len, theta):
             p_pre, dp_pre = (1, 0) if i == 0 else (p_val[:, i - 1], pq_val[:, i - 1])
             p, dp = (0, 0) if i == rep_len - 1 else (p_val[:, i], pq_val[:, i])
             if i < rep_len - 1:
-                temp1 = full_dis * _theta * dp * (dp_pre - dp) / (p_pre - p + 1e-200)
+                temp1 = full_dis * _theta * dp * (dp_pre - dp) / (p_pre - p + 1e-10)
                 ddloglik_val[len_threshold:, i] += np.sum(temp1)
-                temp2 = full_dis * dp ** 2 / (p_pre - p + 1e-200)
+                temp2 = full_dis * dp ** 2 / (p_pre - p + 1e-10)
                 ddloglik_val[i, i] += -np.sum(temp2)
             if i > 0:
-                temp3 = full_dis * _theta * dp_pre * (dp - dp_pre) / (p_pre - p + 1e-200)
+                temp3 = full_dis * _theta * dp_pre * (dp - dp_pre) / (p_pre - p + 1e-10)
                 ddloglik_val[len_threshold:, i - 1] += np.sum(temp3, axis=0)
-                temp4 = full_dis * dp_pre ** 2 / (p_pre - p + 1e-200)
+                temp4 = full_dis * dp_pre ** 2 / (p_pre - p + 1e-10)
                 ddloglik_val[i - 1, i - 1] += -np.sum(temp4)
             if 0 < i < rep_len - 1:
-                ddloglik_val[i, i - 1] = np.sum(full_dis * dp * dp_pre / (p_pre - p + 1e-200))
+                ddloglik_val[i, i - 1] = np.sum(full_dis * dp * dp_pre / (p_pre - p + 1e-10))
             temp5 = full_dis * _theta ** 2 * (dp_pre - dp) ** 2 / (p - p_pre)
             ddloglik_val[-1, -1] += np.sum(temp5, axis=0)
         ddloglik_val += ddloglik_val.transpose() - np.diag(ddloglik_val.diagonal())