項目反応理論のモデル群#

正規累積モデル#

\(i\)番目の被験者の\(j\)番目の項目の値\(y_{ij}\)が二値\(\{0, 1\}\)であるとする(例えば正解・不正解だったり、アンケートの「あてはまる」「あてはまらない」という2件法など)。

\(y_{ij}\)の背後には潜在的な能力の連続量\(\theta_i \in \mathbb{R}\)が存在し、\(\theta_i\)が閾値\(b_j\)を超えていたら1、超えていなければ0が観測されるとする。つまり\(y_{ij}\)が以下のように決まるとする。

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } \theta_i < b_j\\ 1 & \text{ if } \theta_i \geq b_j\\ \end{cases} \end{split}\]

1パラメータ正規累積モデル#

しかし、実際には被験者\(i\)の体調や運(たまたま正解できた)などにより、常にこのようにきれいに正解・不正解が決まるわけではないと考えられる。こうした誤差を表すパラメータ\(\varepsilon_{ij} \sim N(0, \sigma^2_{\varepsilon})\)も追加して

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } (\theta_i - \varepsilon_{ij}) < b_j\\ 1 & \text{ if } (\theta_i - \varepsilon_{ij}) \geq b_j\\ \end{cases} \end{split}\]

とする。誤差が確率変数のため、\(y_{ij}\)のとる値も確率変数として考えることができるようになる。\(b_j\)を移項すると

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } (\theta_i - \varepsilon_{ij} - b_j) < 0\\ 1 & \text{ if } (\theta_i - \varepsilon_{ij} - b_j) \geq 0\\ \end{cases} \end{split}\]

となる。\(\theta_i - \varepsilon_{ij} - b_j \sim N(\theta_i - b_j, \sigma^2_{\varepsilon})\)である。 \(\varepsilon_{ij}\)を移項すれば

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } (\theta_i - b_j) < \varepsilon_{ij}\\ 1 & \text{ if } (\theta_i - b_j) \geq \varepsilon_{ij}\\ \end{cases} \end{split}\]

でもあるので「\(y_{ij}=1\)となるのは誤差\(\varepsilon_{ij}\)\(\theta_i - b_j\)以下のとき」とわかる。

仮に\(\varepsilon_{ij}\)が標準正規分布(\(\sigma^2_{\varepsilon} = 1\)の正規分布)に従うならば、特性値\(\theta_i\)の人が項目\(j\)に当てはまると回答する確率は

\[\begin{split} \begin{aligned} P(y_{ij} = 1) &= P(\varepsilon_{ij} \leq \theta_i - b_j)\\ &= \int_{-\infty}^{\left(\theta_i-b_j\right)} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^2}{2}\right) d z \end{aligned} \end{split}\]

となる(最後のは、\(\varepsilon_{ij}\)が従う標準正規分布のうち \(-\infty\) から \(k\theta_i-b_j\) までの範囲の面積が\(P(\varepsilon_{ij} \leq \theta_i - b_j)\)ということ)。

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
epsilon = np.linspace(-4, 4, 1000)
y = norm.pdf(epsilon)
fig, ax = plt.subplots(figsize=[4,2])
ax.plot(epsilon, y, 'k-')

c = 1

ax.fill_between(epsilon, y, where=(epsilon < c), color='blue', alpha=0.3, label=r'$P(y_{ij} = 1)$')
ax.fill_between(epsilon, y, where=(epsilon >= c), color='red', alpha=0.3, label=r'$P(y_{ij} = 0)$')
ax.axvline(x=c, color='grey', linestyle='--', label=r'$\theta_i - b_j$')
ax.set(
    title=r'$\varepsilon_{ij} \sim N(0, 1)$',
    xlabel=r"$\varepsilon_{ij}$",
    ylabel="density"
)
ax.legend()
ax.grid(True)
fig.show()
../../_images/a949ca679157a2003c80109acfdf5f0e790b3ede5573e318890f4a5e003c0900.png

2パラメータ正規累積モデル#

\(\sigma^2_{\varepsilon}\)が項目ごとに異なる場合を考える。\(\sigma^2_{\varepsilon}=1/a_j\)とすると、誤差の確率分布は

\[ \varepsilon_{ij} \sim N\left(0, \frac{1}{a_j}\right) \]

となる。両辺を\(a_j\)倍すると、\(a_j \varepsilon_{ij} \sim N(0, 1)\)と表すことができ、引き続き標準正規分布を使うことができる。そのためモデルは\(a_i\)が追加され

\[\begin{split} \begin{aligned} P(y_{ij} = 1) &= P(a_j \varepsilon_{ij} \leq \theta_i - b_j)\\ &= \int_{-\infty}^{ a_j (\theta_i-b_j)} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^2}{2}\right) d z \end{aligned} \end{split}\]

となる。

パラメータの意味

\(b_j\)が大きくなると\(\theta_i - b_j\)の値は小さくなり、\(P(y_{ij} = 1)\)の面積が小さくなる。\(y_{ij} = 1\)が正解を表しているとするなら、正答率が低くなる方向に作用する。そのため\(b_j\)項目困難度(item difficulty) と呼ばれる。

項目特性曲線(ICC)との関係でいうと、ICCが0.5になる点が\(b_j = \theta_i\)となる。つまり、困難度が高いほど正答率50%の点にくる\(\theta_i\)も高くなる

また\(a_j\)は値が大きくなると\(\varepsilon_{ij}\)の分散を下げて分布がより尖っていく。また横軸に\(\theta_i - b_j\)、縦軸に\(P(y_{ij} = 1)\)のグラフを書くとき、この曲線の傾きを急にして、\(\theta_i\)が低い人と高い人の間で\(P(y_{ij} = 1)\)の変化を大きくする。例えば、\(\theta_i\)のある点を境にして正解/不正解がはっきり分かれるデータならICCの傾きは急になり、\(a_j\)は高くなる。 そのため\(a_j\)項目識別力(item discrimination) と呼ばれる。

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
x = np.linspace(-4, 4, 1000)

fig, axes = plt.subplots(figsize=[8, 2], ncols=2)


a = 1
for beta in [-1, 0, 1]:
    axes[0].plot(x, norm.cdf(x, loc=beta, scale=1/a), label=r"$b$ = " + f"{beta}")
axes[0].set(xlabel=r"$\theta_i$", ylabel=r"$P(y_{ij} = 1)$", xticklabels=[], yticklabels=[])
axes[0].legend()
axes[0].grid(True)

beta = 0
for a in [0.5, 1, 2]:
    axes[1].plot(x, norm.cdf(x, loc=beta, scale=1/a), label=r"$a$ = " + f"{a}")
axes[1].set(xlabel=r"$\theta_i$", ylabel=r"$P(y_{ij} = 1)$", xticklabels=[], yticklabels=[])
axes[1].legend()
axes[1].grid(True)

fig.show()
../../_images/2a0f992e5a789a6d6cf6c1eb60807988965f6af257b479e51cdd13fa02820755.png

なお横軸に\(\theta_i\)、縦軸に\(P(y_{ij}=1)\)をとったグラフは 項目特性曲線 (item characteristic curve: ICC) と呼ばれる。

因子分析と2P正規累積モデルは等価

因子分析と2パラメータ正規累積モデルは数学的に等価であると知られている。標準化していない(切片が0でない)1因子モデルは、切片\(\tau_j\)、因子負荷量\(a_j\)、因子スコア\(f_i\)

\[ y_{ij} = \tau_j + a_j f_i - \varepsilon_{ij} \quad (\varepsilon_{ij}\sim N(0, \sigma_{\varepsilon})) \]

となる(IRTの説明に合わせて誤差の符号をマイナスにしている)

カテゴリカル因子分析では「離散的な観測変数\(y_{ij}\)はその背後にある潜在的な連続変数によって決まる」という考え方をするため、IRTの冒頭の説明と同じで

\[\begin{split} y_{ij} = \begin{cases} 0 & \text{ if } (\tau_j + a_j f_i) < \varepsilon_{ij}\\ 1 & \text{ if } (\tau_j + a_j f_i) \geq \varepsilon_{ij}\\ \end{cases} \end{split}\]

つまり標準正規分布に従う誤差\(\varepsilon_{ij}\)\(\tau_j + a_j f_i\)より小さいときに\(y_{ij}=1\)となると考えるため、その確率\(P(y_{ij}=1)\)

\[\begin{split} \begin{aligned} P(y_{ij} = 1) &= P(\varepsilon_{ij} \leq \tau_j + a_j f_i)\\ &= \int_{-\infty}^{ (\tau_j + a_j f_i)} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^2}{2}\right) d z \end{aligned} \end{split}\]

となる。2つのモデルのパラメータを\((f_i, a_j, \tau_j) = (\theta_i, a_j, a_j b_j)\)と対応させると

\[\begin{split} \begin{aligned} \tau_j + a_j f_i &= a_j \left( f_i + \frac{\tau_j}{a_j} \right) \\ &= a_j (\theta_i - b_i) \end{aligned} \end{split}\]

となる。

ロジスティックモデル#

正規累積モデルはプロビット回帰と同様のことをするので、コンピュータで積分計算をするのがやや難しいという問題がある。そこでロジスティック分布に置き換えたものが使われる。

ロジスティック分布の確率密度関数と累積分布関数は

\[ f(x)=\frac{\exp (-x)}{[1+\exp (-x)]^2}, \quad F(x)=\frac{1}{1+\exp (-x)} \]

となる。とくに\(x\)を約1.7倍したロジスティック分布は累積分布関数が正規分布と非常に近くなることが知られている。

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, logistic
x = np.linspace(-10, 10, 1000)

fig, axes = plt.subplots(figsize=[12, 3], dpi=90, ncols=2)
ax = axes[0]
ax.plot(x, norm.pdf(x), 'k-', label='Normal Distribution')
ax.plot(x, logistic.pdf(x),  'orange', alpha=0.9, label='Logistic Distribution')

ax.set_title('Normal Distribution vs Logistic Distribution')
ax.set_xlabel('x')
ax.set_ylabel('Probability Density')
ax.legend()
ax.grid(True)

ax = axes[1]
ax.plot(x, norm.cdf(x), 'k-', label=r'CDF of $Normal(x)$')
ax.plot(x, logistic.cdf(x), 'orange', alpha=0.9, label=r'CDF of $Logistic(x)$')
ax.plot(x, logistic.cdf(x * 1.704), 'red', alpha=0.9, label=r'CDF of $Logistic(1.7 x)$')
ax.set_title('Normal Distribution vs Logistic Distribution (CDF)')
ax.set_xlabel('x')
ax.set_ylabel('Probability')
ax.legend()
ax.grid(True)
plt.show()
../../_images/be1480d3d184ef42d63b818995b7efdef53b7fb5a1b31b8b494df28edcfe3682.png

2PLモデル#

正規分布の代わりにロジスティック分布を使った 2パラメータロジスティック(2PL)モデル は以下のように表される。

2PLモデル

\[ P(y_{ij} = 1) = \frac{1}{1+ \exp(-D a_j \big( \theta_i - b_j) \big)} \]
  • \(a_j\):項目識別力

  • \(b_j\):項目困難度

なお\(D\)はロジスティック・シグモイド関数を正規累積モデルの関数に近づけるための定数(通常は\(D=1.7\)\(D=1\)にする)なので、正規累積モデルと比較する必要がなければ不要(\(D=1\)でいい)。

3PLモデル#

例えば4択問題では、正解がわからなくて適当に選んだとしても1/4は当たることになる。こうした影響を「当て推量」パラメータ\(c_j\)として取り入れたモデル。

3PLモデル

\[ P(y_{ij} = 1) = c_j + \frac{1 - c_j}{1 - \exp(-a_j \big( \theta_i - b_j) \big)} \]
  • \(a_j\):項目識別力

  • \(b_j\):項目困難度

  • \(c_j\):当て推量

\(c_j\)は項目特性曲線の下限となる。\(\theta_i\)がどんなに低い人でも必ず\(c_j\)以上の\(P(y_{ij} = 1)\)になるということ。

4PLモデル#

項目特性曲線の上限を表すパラメータ\(d_j\)を追加したもの。\(\theta_i\)がどんなに高い人でも100%の正答率にはできない高難度な状況(運ゲー)を想定したモデル。

4PLモデル

\[ P(y_{ij} = 1) = c_j + \frac{d_j - c_j}{1 - \exp(-a_j \big( \theta_i - b_j) \big)} \]
  • \(a_j\):項目識別力

  • \(b_j\):項目困難度

  • \(c_j\):当て推量。項目特性曲線の下限

  • \(d_j\):項目特性曲線の上限

4PLMになるとかなりモデルが複雑になりパラメータの推定も不安定になるので、1~3PLMほど一般的ではない。対応していないライブラリも多い。

5PLモデル#

「非対称性」のパラメータ\(e_j\)を追加したもの。4PLまでは項目特性曲線の動き方が0.5を中心に対称になっている。5PLでは「最初は\(\theta_i\)があがるほど急激に\(P(y_{ij}=1)\)が上がるが、徐々に上がりにくくなる」などの状況を表すことができる。

5PLモデル

\[ P(y_{ij} = 1) = c_j + \frac{d_j - c_j}{\left[ 1 - \exp(-a_j \big( \theta_i - b_j) \big) \right]^{e_j}} \]
  • \(a_j\):項目識別力

  • \(b_j\):項目困難度

  • \(c_j\):当て推量。項目特性曲線の下限

  • \(d_j\):項目特性曲線の上限

  • \(e_j\):非対称性

多値型モデル#

項目の反応\(y_{ij}\)が多値になった場合のモデルも存在する。

段階反応モデル#

段階反応モデル(graded response model: GRM) は複数の二値IRTモデルを組み合わせて多値反応を表現する。

回答者\(i\)の項目\(j\)に対する回答\(y_{ij} = k \quad (k=1,2,\dots,K)\)について、「\(k\)以上のカテゴリを選ぶ確率」を考えると、これはまだ「\(k\)未満 or \(k\)以上」の二値なので2PLなどで表せる。例えば以下のようになる。

\[ P(y_{ij} \geq k) = \frac{1}{1+ \exp \big(-a_j ( \theta_i - b_{jk}) \big)} \]

なお、困難度は項目\(j\)のカテゴリ\(k\)ごとに用意されるため\(b_{jk}\)に変更している。

このモデルを組み合わせると、「ちょうど\(k\)番目のカテゴリを選ぶ確率」は

\[ P(y_{ij} = k) = P(y_{ij} \geq k) - P(y_{ij} \geq k + 1) \]

と表すことができる。ただし端のカテゴリは\(P(y_{ij} \geq 1) = 1, P(y_{ij} \geq K + 1) = 0\)とする。また確率100%の困難度は低くて当然なので\(b_{j1} = -\infty\)とする。

段階反応モデル

\[ P(y_{ij} = k) = \frac{1}{1+ \exp(-a_j \big( \theta_i - b_{jk}) \big)} - \frac{1}{1+ \exp(-a_j \big( \theta_i - b_{jk+1}) \big)} \]

名義反応モデル#

名義反応モデル(nominal response model: NRM) も段階反応モデルと同様に多値の回答にIRTを適用したモデル。softmax関数のような形で多値化する。

\[ P(y_{ij} = k) = \frac{ \exp (a_{jk} \theta_i + \gamma_{jk}) }{ \sum_{\kappa=1}^K \exp (- a_{j\kappa} \theta_i + \gamma_{j\kappa}) } \]

連続反応モデル#

連続反応モデル(Continuous Response Model, CRM)は連続的な反応(例:時間、強度、割合)を扱うことができる。

CRMは段階反応モデル(GRM)をさらに多段階に拡張していって連続値を扱う。

受験者\(i\)が項目\(j\)において\(x_j\)点以上をとる確率は、次のようにロジスティック関数で表すことができる

\[ p_i^*(x_j)=\frac{1}{ 1 + \exp \big(- a_j (\theta_i - b_{x_j}) \big)} \]

これを使うと、受験者\(i\)が項目\(j\)において\(x_j\)点をとる確率 \(p_i(x_j)\) を次のように定義できる

\[ p_i(x_j) = \lim _{\Delta x_j \to 0} \frac{p_i^*(x_j)-p_i^*(x_j+\Delta x_j)}{\Delta x_j} \]

\(b_{x_j}\)は項目\(j\)において\(x_j\)以上の点を取る困難度で、ロジスティック関数の逆関数を用いて以下のように定義される

\[ b_{x_j} = \beta_j + \frac{1}{\alpha_j} \log \frac{x_j}{K_j-x_j} \]

ここで\(\alpha_j, \beta_j\)はロジスティック関数の逆関数における識別力と困難度を表現するパラメータ、\(K_j\)\(x_j\)のとる最大値(\(x_j \in [0, K_j]\)

多次元モデル#

多次元因子分析モデルでは項目反応\(y_{ij}\)

\[ y_{i j} = \tau_j + \sum_{t=1}^T f_{i t} b_{t j} - \varepsilon_{ij} \]

のように表す。ここから考えると、多次元IRTモデル(の2PLM)の項目反応関数は

\[ P(y_{ij}=1) = \frac{1}{ 1 + \exp \left[ -\sum_{t=1}^T \alpha_{j t} \theta_{i t} - \tau_j \right]} \]

のように考えられる。\(\tau_j\)は切片のパラメータで、これを困難度に変換するには、次元ごとの識別力を一つにまとめた 多次元識別力 である

\[ \operatorname{MDISC}_j=\sqrt{\sum_{t=1}^T \alpha_{t j}^2} \]

を使って

\[ \beta_j=\frac{-\tau_j}{\mathrm{MDISC}_j} \]

と変換する(豊田 2013 『項目反応理論[中級編]』) 。

補償型と非補償型#

ある次元の特性値が低くても、別の次元の特性値がとても高い場合は反応確率が高くなるようなモデルは 補償型 (compensatory) モデル と呼ばれる。例えば先述の

\[ P(y_{ij}=1) = \frac{1}{ 1 + \exp \left[ -\sum_{t=1}^T \alpha_{j t} \theta_{i t} - \tau_j \right]} \]

は補償型モデルである。

一方で複数の次元の特性値がそろわないと高くならない多次元IRTモデルは 非補償 (noncompensatory) モデル あるいは 部分補償(partially compensatory)モデル と呼ばれる(例:「総合学力」に対して「数学力スコア」と「英語力スコア」はそれぞれ独立していると考える)。最もシンプルな非補償モデルの項目反応関数は

\[ P(y_{ij}=1) = \prod_{t=1}^T \frac{1}{1+\exp \left[-\alpha_{j t}\left(\theta_{i t}-\beta_{j t}\right)\right]} \]

である。それぞれの次元の特性値\(\theta_{i t}\)で説明される確率の積にして同時確率のようになっている。

推定の難しさなどの観点から補償型のほうがよく使われるらしい

適合度の評価#

局所独立性の確認#

  • \(\chi^2\)統計量

個人適合度#

\(\theta_i\)が高い人なのに困難度が低い項目で間違えるのはおかしい」といった考え方から、特性\(\theta_i\)と困難度\(b_j\)の関係性を見る。

  • \(z_h\)\(l_z\))統計量 :ある回答者の反応パターン\(\boldsymbol{y}_i = (y_{i1}, y_{i2}, \dots, y_{iJ})\) と「考えられる全反応パターン」での尤度を比較する

IRTの前提条件・仮定#

局所独立性(local independence)の仮定#

\(\theta_i\)で条件づけたとき、項目\(j\)と項目\(l\)\(j \neq l\)) への回答は完全に独立であるという仮定。最尤推定の計算の簡易化のために必要になる。

仮定を満たすか確認する方法:

  • 例えばデータから\(\theta_i\)が特定の値の人を集めて(\(\theta_i\)で条件づけたデータを用意して)、項目\(j,l\)の回答のクロス表を作ったとき、関連性がない場合は局所独立性が満たされていることがわかる。

一次元性#

もともとのIRTは単一のテストに対する単一の能力を測定することが主たる目的のため、推定する因子数は1である(因子数を拡張したモデルも存在するが、基本は1)。そのため、すべての項目が同じ能力・特性のみを反映していることが必要となる。

仮定を満たすか確認する方法:

  • 例えば、1因子の確認的因子分析を行い、適合度指標(RMSEAやCFA、AGFIなど)が許容範囲内かどうかで確認する。

  • DIMTEST(Stout et al., 1996)や DETECT(Zhang & Stout, 1999)といった手法もある

PythonでのIRTの実装#

一般的なモデル → pyirtなどパッケージがある

複雑な、独自のモデル → PyStanやPyMCでベイズ推定するのが良さそう

2PL#

データの生成#

Hide code cell source

# ダミーデータの生成
# 出所: https://qiita.com/takuyakubo/items/43d56725952e67032b49
import numpy as np
from functools import partial

# 3 parameter logistic model の定義
def L3P(a, b, c, x):
    return c + (1 - c) / (1 + np.exp(-  a * (x - b)))

# model parameterの定義
# aは正の実数, bは実数, cは0より大きく1未満であれば良い
a_min = 0.3
a_max = 1

b_min = -2
b_max = 2

c_min = 0
c_max = .4

# 問題数と回答者数
num_items = 10
num_users = 1000

# 問題parameterの生成
item_params = np.array(
    [np.random.uniform(a_min, a_max, num_items),
     np.random.uniform(b_min, b_max, num_items),
     np.random.uniform(c_min, c_max, num_items)]
).T

# 受験者parameterの生成
user_params = np.random.normal(size=num_users)

# 項目反応行列の作成、 要素は1(正答)か0(誤答)
# i行j列は問iに受験者jがどう反応したか
ir_matrix_ij = np.vectorize(int)(
    np.array([partial(L3P, *ip)(user_params) > np.random.uniform(0, 1, num_users) for ip in item_params])
)

import pandas as pd

df = pd.DataFrame(ir_matrix_ij.T,
                  index=[f"user_{i+1}" for i in range(num_users)],
                  columns=[f"question_{j+1}" for j in range(num_items)])

df.head()
question_1 question_2 question_3 question_4 question_5 question_6 question_7 question_8 question_9 question_10
user_1 1 0 1 0 1 0 1 1 1 1
user_2 1 1 1 1 1 0 1 0 0 1
user_3 1 1 1 0 1 1 1 0 0 1
user_4 1 1 1 0 0 1 1 0 0 1
user_5 1 0 1 1 0 1 0 0 1 0
# 縦持ちへ変換
df_long = pd.melt(
    df.reset_index(),
    id_vars="index",
    var_name="item",
    value_name="response",
).rename(columns={"index": "user"})
df_long.head()
user item response
0 user_1 question_1 1
1 user_2 question_1 1
2 user_3 question_1 1
3 user_4 question_1 1
4 user_5 question_1 1

モデルの定義#

# indexと値の取得
user_idx, users = pd.factorize(df_long["user"])
item_idx, items = pd.factorize(df_long["item"])
responses = df_long["response"].to_numpy()

import pymc as pm
model = pm.Model(coords={"user": df.index, "item": df.columns})
with model:
    # 観測値の配列
    response_obs = pm.ConstantData("responses", responses)

    # 2PLM
    theta = pm.Normal("theta", mu=0.0, sigma=1.0, dims="user")
    a = pm.HalfNormal("a", sigma=1.0, dims="item")
    beta = pm.Normal("beta", mu=0.0, sigma=1.0, dims="item")
    logit_p = pm.Deterministic("logit_p", a[item_idx] * (theta[user_idx] - beta[item_idx]))

    # ベルヌーイ分布(pではなくlogit_pの引数に渡すことでシグモイド関数の計算をpymc側にまかせている)
    obs = pm.Bernoulli("obs", logit_p=logit_p, observed=response_obs)

g = pm.model_to_graphviz(model)
g
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[6], line 10
      7 model = pm.Model(coords={"user": df.index, "item": df.columns})
      8 with model:
      9     # 観測値の配列
---> 10     response_obs = pm.ConstantData("responses", responses)
     12     # 2PLM
     13     theta = pm.Normal("theta", mu=0.0, sigma=1.0, dims="user")

AttributeError: module 'pymc' has no attribute 'ConstantData'

推定#

%%time
with model:
    idata = pm.sample(random_seed=0)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [theta, a, beta]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 56 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
CPU times: user 12.9 s, sys: 718 ms, total: 13.6 s
Wall time: 1min 4s
idata
arviz.InferenceData
    • <xarray.Dataset> Size: 353MB
      Dimensions:        (chain: 4, draw: 1000, user: 1000, item: 10,
                          logit_p_dim_0: 10000)
      Coordinates:
        * chain          (chain) int64 32B 0 1 2 3
        * draw           (draw) int64 8kB 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
        * user           (user) <U9 36kB 'user_1' 'user_2' ... 'user_999' 'user_1000'
        * item           (item) <U11 440B 'question_1' 'question_2' ... 'question_10'
        * logit_p_dim_0  (logit_p_dim_0) int64 80kB 0 1 2 3 4 ... 9996 9997 9998 9999
      Data variables:
          theta          (chain, draw, user) float64 32MB 0.009313 0.2246 ... -0.2381
          beta           (chain, draw, item) float64 320kB -1.785 -0.8052 ... -1.193
          a              (chain, draw, item) float64 320kB 0.6816 0.1785 ... 0.2549
          logit_p        (chain, draw, logit_p_dim_0) float64 320MB 1.223 ... 0.2435
      Attributes:
          created_at:                 2025-03-03T08:40:11.496294+00:00
          arviz_version:              0.20.0
          inference_library:          pymc
          inference_library_version:  5.17.0
          sampling_time:              56.205814838409424
          tuning_steps:               1000

    • <xarray.Dataset> Size: 496kB
      Dimensions:                (chain: 4, draw: 1000)
      Coordinates:
        * chain                  (chain) int64 32B 0 1 2 3
        * draw                   (draw) int64 8kB 0 1 2 3 4 5 ... 995 996 997 998 999
      Data variables: (12/17)
          step_size_bar          (chain, draw) float64 32kB 0.2212 0.2212 ... 0.2179
          diverging              (chain, draw) bool 4kB False False ... False False
          energy                 (chain, draw) float64 32kB 8.004e+03 ... 8.084e+03
          perf_counter_start     (chain, draw) float64 32kB 1.098e+06 ... 1.098e+06
          index_in_trajectory    (chain, draw) int64 32kB 8 5 -4 12 -10 ... -5 13 9 -7
          max_energy_error       (chain, draw) float64 32kB 2.905 0.5263 ... 0.432
          ...                     ...
          acceptance_rate        (chain, draw) float64 32kB 0.3001 0.8956 ... 0.8464
          perf_counter_diff      (chain, draw) float64 32kB 0.01731 ... 0.01731
          n_steps                (chain, draw) float64 32kB 15.0 15.0 ... 15.0 15.0
          energy_error           (chain, draw) float64 32kB 0.2613 -0.1657 ... 0.1819
          process_time_diff      (chain, draw) float64 32kB 0.01732 ... 0.01731
          step_size              (chain, draw) float64 32kB 0.1867 0.1867 ... 0.2182
      Attributes:
          created_at:                 2025-03-03T08:40:11.515943+00:00
          arviz_version:              0.20.0
          inference_library:          pymc
          inference_library_version:  5.17.0
          sampling_time:              56.205814838409424
          tuning_steps:               1000

    • <xarray.Dataset> Size: 160kB
      Dimensions:    (obs_dim_0: 10000)
      Coordinates:
        * obs_dim_0  (obs_dim_0) int64 80kB 0 1 2 3 4 5 ... 9995 9996 9997 9998 9999
      Data variables:
          obs        (obs_dim_0) int64 80kB 1 1 0 0 1 1 1 1 0 1 ... 1 1 0 1 0 1 0 1 0
      Attributes:
          created_at:                 2025-03-03T08:40:11.525280+00:00
          arviz_version:              0.20.0
          inference_library:          pymc
          inference_library_version:  5.17.0

    • <xarray.Dataset> Size: 120kB
      Dimensions:          (responses_dim_0: 10000)
      Coordinates:
        * responses_dim_0  (responses_dim_0) int64 80kB 0 1 2 3 ... 9997 9998 9999
      Data variables:
          responses        (responses_dim_0) int32 40kB 1 1 0 0 1 1 1 ... 1 0 1 0 1 0
      Attributes:
          created_at:                 2025-03-03T08:40:11.527285+00:00
          arviz_version:              0.20.0
          inference_library:          pymc
          inference_library_version:  5.17.0

一部の項目の\(a_j, b_j\)

Hide code cell source

import arviz as az

query = {"item": ["question_1", "question_2"]}

az.plot_trace(idata, coords=query, var_names=["a", "beta"], figsize=[4, 4])
plt.tight_layout()
plt.show()
../../_images/847d3dc8c7accadbbb682ec47ac405b7c071f600e5278f78abf25459e4a09afc.png

一部の回答者の\(\theta_i\)

Hide code cell source

import arviz as az

query = {"user": ["user_1", "user_2"]}

az.plot_trace(idata, coords=query, var_names=["theta"], figsize=[4, 2])
plt.tight_layout()
plt.show()
../../_images/5c21829b51ba0494adc1c777272e996fc9c6c1370c41132a44abb67a4b52a643.png

参考#