段階反応モデル（Graded Response Model、GRM）は、5件法のように順序のある選択肢への回答から、回答者の特性θと項目の識別力・境界困難度を推定する項目反応理論（IRT）のモデルです。この記事では、PyMC 6.3.2 でGRMをベイズ推定するコードを、実際に走らせた出力つきで示します。
掲載したコードと数値は、Python 3.12・PyMC 6.3.2・ArviZ 1.3.0 で疑似データ（300人×8項目×5件法）を推定して得たものです。境界の順序制約を外したとき、識別力の符号を縛らなかったときに何が起きるかも、エラー文と診断値で載せています。
まとめ：段階反応モデルをPyMCで推定する要点
- GRMは「カテゴリk以上を選ぶ確率」を P(Y≥k|θ)=1/(1+exp(−a(θ−bk))) で表し、隣り合う累積確率の差をカテゴリ確率にします（Samejima, 1969）。
- PyMCでは
pm.OrderedLogisticに eta=a·θ、cutpoints=a·b を渡せばGRMになります。観測値は0始まり（5件法なら0〜4）に直します。 - 境界bには
transform=pm.distributions.transforms.ordered、識別力aには LogNormal など正の事前分布を置きます。どちらかを外すと初期化の時点でSamplingErrorになりました。 - 判定は R-hat<1.01、ESS>400、発散0件が目安です（Vehtari et al., 2021）。実測ではchains=4・draws=1000・target_accept=0.9 で R-hat最大1.004、ESS最小1,100、発散0件でした。
- 大標本で標準的な一次元GRMを当てるだけなら、Rの mirt（
itemtype='graded'）の周辺最尤推定を先に検討します。PyMCを選ぶ理由は、少標本・疎なカテゴリ・階層化など事前分布が効く場面です。
段階反応モデル（GRM）の定義と累積確率の考え方
GRMは Fumiko Samejima が1969年に Psychometrika のモノグラフ（Psychometric Monograph Supplement No.17、34巻S1号）で示したモデルです。二値の2パラメータロジスティックモデル（2PL）を、順序のある3つ以上のカテゴリへ広げた形になっています。
K件法の項目jでは、カテゴリ境界ごとに「k以上を選ぶ確率」を2PLの曲線で書きます。
P(Y_ij ≥ k | θ_i) = 1 / (1 + exp(−a_j (θ_i − b_jk))) k = 1, …, K−1
P(Y_ij = k | θ_i) = P(Y_ij ≥ k | θ_i) − P(Y_ij ≥ k+1 | θ_i)
（P(Y ≥ 0) = 1、P(Y ≥ K) = 0 とおく）パラメータは項目ごとに識別力aが1つ、境界困難度bがK−1個です。5件法なら1項目あたり5個、8項目で40個になり、これに回答者数ぶんのθが加わります。K=2にすると境界が1本になり、2PLと一致します。
カテゴリ確率を差で作るため、bj1<bj2<…と昇順でないと確率が負になります。PyMCで順序制約を必ず入れる理由はここにあります。
GRMとPCM・RSM・名義反応モデルの違い
順序カテゴリを扱うIRTモデルは境界の定義が異なり、同じデータでもパラメータの意味が変わります。
| モデル | 原典 | 境界の定義 | 識別力 | 順序の前提 |
|---|---|---|---|---|
| 段階反応モデル（GRM） | Samejima（1969） | 累積（k以上 vs 未満） | 項目ごと | あり |
| 部分採点モデル（PCM） | Masters（1982） | 隣接カテゴリ間 | 全項目で共通（Rasch系） | あり |
| 評定尺度モデル（RSM） | Andrich（1978） | 隣接カテゴリ間 | 全項目で共通（Rasch系） | あり・閾値構造も共通 |
| 名義反応モデル（NRM） | Bock（1972） | 選択肢ごと | 選択肢ごと | なし |
項目によって弁別の鋭さが違う心理尺度ならGRMが第一候補です。PCMはRSMを「項目ごとにカテゴリ数や閾値構造が違ってよい」形へ拡張したモデルで、Rasch系の性質（素点が十分統計量になる等）を尺度設計の前提にしたいときに選びます。選択肢に順序がない多肢選択の誤答分析なら名義反応モデルです。
ベイズ推定でGRMを扱う理由と最尤推定（mirt）との使い分け
GRMの推定法は、Rの mirt などが採る周辺最尤推定と、事後分布をMCMCで描くベイズ推定に分かれます。最尤推定の考え方は尤度とは？読み方（ゆうど）・確率との違いから尤度関数・最尤推定・尤度比までわかりやすく解説で扱っています。
ベイズ推定が効くのは、データが境界を決めるのに足りない場面です。今回の疑似データでも、項目0は300人中「カテゴリ1」を選んだのが15人、「カテゴリ3」が14人しかいませんでした。このような疎なカテゴリでは、最尤推定の境界値が大きく振れたり、ある境界より上を誰も選ばない項目で推定値が発散したりします。事前分布 b〜Normal(0, 2) は極端な境界値を抑える方向に働きます。ただし、この尺度で妥当な設定かどうかは、事前予測と事前分布を変えた感度分析で確認します。
結果が事後分布で返るので、θの区間をそのまま個人の判定や短縮版の項目選定に使える点も利点です。項目パラメータに階層事前分布を置いて項目間で情報を共有するといった拡張も、モデル記述を数行足すだけで済みます。
逆に、回答者が数千人規模で、一次元のGRMを当てて項目パラメータを読みたいだけなら、PyMCを選ぶ理由は薄くなります。今回の300人×8項目でも、4チェーンのサンプリングに51〜99秒（2回の実行）かかりました。周辺最尤推定はθを積分して消すので、推定対象は項目パラメータ（この例では40個）だけです。NUTSは人数ぶんのθも同時にサンプリングするため、回答者が増えるほど推定するパラメータ数が増えます。処理時間は同じデータと環境で比べて決めてください。
PyMC 6でGRMを組む実装手順
検証環境とPyMC 6系で変わった点
PyMCは2026年5月13日に6.0.0が公開され、5系は2026年5月1日の5.28.5が最終版です。この記事では2026年9月8日公開の6.3.2と、2026年8月11日公開の ArviZ 1.3.0 を使いました（どちらもPython 3.12以上が必要）。5系の解説記事のコードを写すときに引っかかる差は次の3点です。
pm.sample()の戻り値は xarray のDataTreeです。idata.posterior、idata.posterior_predictiveのように子ノードとして参照します。- NUTSの設定は
nuts={"target_accept": 0.9}で渡す形がdocstringで案内されており、nuts_sampler_kwargsは非推奨（Deprecated）と書かれています。6.3.2ではtarget_accept=0.9を直接渡してもエラーにはなりませんでした。 - ArviZ 1.3.0 の
az.summary()は、既定の区間が89%の等裾区間（列名eti89_lb・eti89_ub）です。0.x系のhdi_3%・hdi_97%列を前提にした後処理は書き換えが要ります。
Intel Mac（x86_64）に uv で入れたところ、依存の llvmlite 0.49.0 にはこの環境向けのホイールが無く、ソースビルドになって失敗しました。--only-binary llvmlite,numba でソースビルドを禁止すると、ホイールのある版へ解決し直され、numba 0.62.1 が入りました。
uv venv -p 3.12 .venv
uv pip install -p .venv/bin/python pymc==6.3.2 arviz==1.3.0 matplotlib --only-binary llvmlite,numbamatplotlib はPyMCの必須依存に入っていませんが、6.3.2の進捗バー表示が import するため、入れずに pm.sample() を呼ぶと ModuleNotFoundError: No module named 'matplotlib' が発生します。入れない場合は progressbar=False を付けます。
回答データをロング形式と0始まりへ整形する工程
PyMCの OrderedLogistic は観測値を0〜K−1の整数で受け取ります。1〜5で記録した5件法なら、1を引いてから渡します。
回答行列（人×項目）は、1行=1回答のロング形式へ展開します。回答者の番号 person、項目の番号 item、回答 y の3本の配列にします。欠測メカニズムを無視できる条件なら、無回答セルの行を除外し、観測された回答だけで尤度を構成できます。無回答が未観測の回答内容に依存する場合は、欠測過程のモデル化も検討します。逆転項目は、この段階で K−1−y に変換して向きを揃えます。
OrderedLogisticへ識別力と境界を渡すパラメータ化
PyMCの OrderedLogistic(eta, cutpoints) は、P(Y≥k)=logistic(eta−cutpoints[k−1]) という累積ロジットです。GRMの a(θ−bk) を展開すると aθ−abk なので、eta=a·θ、cutpoints=a·b を渡せば式が一致します。
このとき a>0 なら、bが昇順である限り a·b も昇順に保たれます。境界をbのまま事前分布に置くので、「境界困難度」をθと同じ尺度で解釈でき、事前分布も Normal(0, 2) のように決めやすくなります。cutpoints を直接パラメータにする傾き切片型も書けますが、その場合は後述する符号の不定性に注意が要ります。
次は疑似データの生成から要約までのコードです。後掲の事後予測チェックと項目情報量のコードは、このファイルのmain()内の末尾に、各行を4スペース字下げして追加してください。両コードはmain()内で定義した変数を使います。
import numpy as np
import pymc as pm
import arviz as az
def main():
# 1. 疑似データ（300人 x 8項目 x 5件法）
rng = np.random.default_rng(42)
N, J, K = 300, 8, 5
theta_true = rng.normal(0, 1, N)
a_true = rng.lognormal(0.2, 0.3, J)
b_true = np.sort(rng.normal(0, 1, (J, K - 1)), axis=1)
p_ge = 1 / (1 + np.exp(-a_true[None, :, None] * (theta_true[:, None, None] - b_true[None, :, :])))
u = rng.uniform(size=(N, J, 1))
Y = (u < p_ge).sum(axis=2) # 0..K-1 の整数
# 2. ロング形式へ（欠測があれば行ごと落とすだけで済む）
person, item = np.indices(Y.shape)
person, item, y = person.ravel(), item.ravel(), Y.ravel()
coords = {"person": np.arange(N), "item": np.arange(J), "cut": np.arange(K - 1)}
with pm.Model(coords=coords) as grm:
theta = pm.Normal("theta", 0.0, 1.0, dims="person")
a = pm.LogNormal("a", 0.0, 0.5, dims="item")
b = pm.Normal(
"b", 0.0, 2.0, dims=("item", "cut"),
transform=pm.distributions.transforms.ordered,
initval=np.tile(np.linspace(-1.5, 1.5, K - 1), (J, 1)),
)
eta = a[item] * theta[person]
cutpoints = a[item][:, None] * b[item]
pm.OrderedLogistic("y", eta=eta, cutpoints=cutpoints, observed=y)
idata = pm.sample(1000, tune=1000, chains=4, random_seed=1, nuts={"target_accept": 0.9})
summ = az.summary(idata, var_names=["a", "b"])
print(summ.head(12).to_string())
print("max r_hat", float(summ["r_hat"].max()))
print("min ess_bulk", float(summ["ess_bulk"].min()), "min ess_tail", float(summ["ess_tail"].min()))
print("divergences", idata["sample_stats"]["diverging"].sum().item())
if __name__ == "__main__":
main()θの事前分布を Normal(0, 1) に固定しているのは、尺度の原点と単位を決めるためです。GRMはθとbを同じ量だけ平行移動しても、θとbを同じ正の定数倍にしてaをその定数で割っても尤度が変わりません。尺度を識別する制約がないと、位置と尺度の不定性が推定を難しくします。本記事ではθの母集団分布を標準正規分布に固定しますが、項目パラメータを固定する識別方法もあります。
ordered変換の初期値 initval を厳密な昇順で与えるのは、隣接差が0以下になる初期値を避けるためです。等間隔にしているのは、ランダムな初期化で境界が極端な位置から始まるのを避けるためです。
NUTSの実行設定とmainガード
NUTSの target_accept は既定0.8です。上げるほどステップ幅が小さくなり、発散は減る一方で1回の反復は遅くなります。PyMCの公式ドキュメントは、ウォームアップ後に発散が1件でも残る場合は結果が信頼できない可能性があるとして、再パラメータ化するか0.9〜0.95へ上げるよう案内しています。GRMはθの数が多く事後分布の形が場所によって変わるので、最初から0.9で回しました。chains=4 は、Vehtari et al.（2021）が既定として推奨する本数です。
コードを if __name__ == "__main__": で包んでいるのには理由があります。chains=4 は別プロセスで並列に走るため、macOSのように fork 以外の方式で子プロセスを起動する環境では、ガードがないと子プロセスがスクリプト全体を再実行します。ガードなしのスクリプトで chains=2・cores=2 のサンプリングを実行すると、freeze_support() を案内する RuntimeError が出て止まりました。Windowsも fork を使わない起動方式なので、.py として実行するなら必須と考えてください。
収束診断と事後予測チェックの判定基準
MCMCの結果は、収束したかを数値で確かめてから読みます。MCMCそのものの仕組みと収束判定の考え方はMCMC（マルコフ連鎖モンテカルロ法）とは？仕組み・MH法とHMCの違い・収束判定をPythonで解説にまとめています。
R-hat・ESS・発散の3指標と実測値
Vehtari et al.（2021）は、順位正規化した R-hat が1.01未満のときだけ標本を使うこと、順位正規化した ESS が400を超えることを推奨しています。GelmanとRubin（1992）の1.1より厳しい基準で、1.1未満でも収束していない例が見つかったことが理由として挙げられています。ESSは分布の中心を見る ess_bulk と、区間の端を見る ess_tail の両方を確認します。
上のコードを実行した結果の先頭12行が次です。
mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat mcse_mean mcse_sd
a[0] 1.59 0.243 1.2 2 1462 2073 1.00 0.0064 0.0043
a[1] 0.904 0.154 0.67 1.2 1605 1750 1.00 0.0038 0.0022
a[2] 1.36 0.199 1.1 1.7 1103 1729 1.00 0.006 0.0034
a[3] 1.25 0.202 0.94 1.6 1566 1867 1.00 0.0051 0.0032
a[4] 0.821 0.125 0.63 1 1247 1295 1.00 0.0035 0.0021
a[5] 0.879 0.14 0.67 1.1 1217 1686 1.00 0.004 0.0024
a[6] 0.967 0.169 0.7 1.2 1442 1642 1.00 0.0044 0.0028
a[7] 1.47 0.212 1.1 1.8 1138 1635 1.00 0.0063 0.0039
b[0, 0] 0.577 0.118 0.4 0.78 1577 2173 1.00 0.003 0.0019
b[0, 1] 0.801 0.132 0.61 1 1621 2078 1.00 0.0033 0.0022
b[0, 2] 1.72 0.205 1.4 2.1 1574 1871 1.00 0.0052 0.0037
b[0, 3] 2.2 0.255 1.8 2.6 1638 2226 1.00 0.0064 0.0048
max r_hat 1.0038202526158047
min ess_bulk 1099.8272407270088 min ess_tail 1295.620991042581
divergences 0a・bの全40パラメータで R-hat は最大1.004、ESSは最小でも1,100、発散は0件で、3指標とも基準を満たしました。θ（300人分）も同じ設定の実行で R-hat最大1.006、ess_bulk最小3,087でした。個人のθを解釈するときは az.summary(idata, var_names=["theta"]) で同じ3指標を確かめます。r_hat 列は既定の自動表示で小数点以下2桁に丸めて表示されるため、1.01と表示された場合は実際の値が1.01未満かどうかを表示だけでは判断できません。上のように列の最大値を数値で出して確かめます。
同じモデルを draws=500・tune=500・target_accept 既定（0.8）で回すと、「The rhat statistic is larger than 1.01 for some parameters」という警告が出ました。a・bは R-hat最大1.006で基準内でしたが、θは300人中3人で1.01を超え、最大1.012でした。θは1人あたりの回答が8件しかなく事後分布が広いため、反復を削ると先にθ側が基準を割ります。
疑似データで確かめたパラメータの再現度
真値がわかっている疑似データなので、推定値との一致も確かめました。識別力aの事後平均と真値の相関は0.89、θの事後平均と真値の相関は0.87です。
| 項目 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| a の真値 | 2.05 | 0.77 | 1.58 | 1.11 | 1.20 | 0.89 | 1.10 | 1.80 |
| a の事後平均 | 1.59 | 0.90 | 1.36 | 1.25 | 0.82 | 0.88 | 0.97 | 1.47 |
項目0（真値2.05）と項目4（真値1.20）は、真値が89%区間（1.2〜2.0、0.63〜1.0）の外にありました。300人×8項目の規模では、識別力の点推定が真値から0.4程度ずれることがあります。項目の取捨を識別力の小数第1位で決めるような使い方は、この規模のデータでは避けるべきです。
事後予測チェックでカテゴリ度数を照合する手順
収束しても、モデルがデータに合っているとは限りません。事後予測チェックでは、推定したパラメータから回答を生成し直し、観測した度数がその分布の中に収まるかを見ます。
with grm:
pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=5)
yrep = idata.posterior_predictive["y"].values # (chain, draw, 回答数)
mask = item == 0
obs = np.bincount(y[mask], minlength=K)
rep = np.stack([np.bincount(r[mask], minlength=K) for r in yrep.reshape(-1, yrep.shape[-1])])
lo, hi = np.percentile(rep, [3, 97], axis=0)本番と同じ設定（draws=1000・target_accept=0.9）で回した項目0の結果は次のとおりで、5カテゴリすべてで観測度数が事後予測の94%区間に収まりました。
obs [200 15 47 14 24]
lo [178. 7. 33. 7. 14.]
hi [214. 26. 66. 26. 35.]実データで観測度数が区間から外れる項目があれば、カテゴリの使われ方がGRMの累積構造と合っていない（中間カテゴリだけ選ばれやすい等）か、多次元性が混じっている可能性を疑います。
推定が壊れる失敗パターンと実測した症状
GRMをPyMCで組むとき、制約の入れ忘れはエラーか収束不全のどちらかで表面化します。同じ疑似データで、制約を1つずつ外して実行した結果です。
境界の順序制約を外したときの初期化エラー
bから transform=pm.distributions.transforms.ordered を外して実行すると、サンプリングが始まる前に止まりました。
pymc.exceptions.SamplingError: Initial evaluation of model at starting point failed!
Logp initial evaluation results:
{'theta': np.float64(-321.61), 'a': np.float64(-11.74), 'b': np.float64(-53.38), 'y': np.float64(-inf)}初期値の境界が昇順になっていない項目があると、カテゴリ確率が負になり、観測の対数尤度 y が −inf になります。このエラー文で y だけが −inf なら、まず境界の順序を疑います。model.debug() を呼ぶと、どの値で計算が壊れたかを詳しく出せます。
識別力に負の値を許したときの符号の不定性
aの事前分布を LogNormal から Normal(0, 1) に替えると、本記事のパラメータ化（cutpoints=a·b）では、上と同じ SamplingError で止まりました。aが負の項目では a·b が降順になり、OrderedLogistic の前提が崩れるためです。
より厄介なのは、cutpoints を直接パラメータにした傾き切片型（P(Y≥k)=logistic(aθ−ck)）で a〜Normal(0, 1) としたときです。aとθの符号を同時に反転しても尤度が変わらないため、エラーは出ずにサンプリングが完走し、チェーンごとに別の山へ収まりました。
mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat
a[0] 0 1.65 -1.9 1.9 6 39 1.73
a[1] 0 1.021 -1.2 1.2 6 40 1.73
...
chain means a[0]: [ 1.64 -1.63 -1.62 1.63]4本のチェーンのうち2本は a[0]≒+1.64、2本は≒−1.63 に張り付き、識別力aのR-hatの表示値は1.73〜1.74、ess_bulkの表示値は全8項目で6でした。事後平均は0になり、「識別力0の項目」という誤った結論が出ます。R-hatが1.7前後に揃い、チェーン平均が正負で対称になっていたら、符号の不定性です。aを LogNormal や HalfNormal で正に縛れば解消します。
疎なカテゴリと境界の接近
項目0の境界は b[0,0]=0.58 と b[0,1]=0.80 で、差は0.22しかありません。カテゴリ1の選択確率は2本の累積曲線の差であり、θがこの区間に入る人数とは一致しません。近接した境界によって選択確率が低くなり、この疑似データでの観測度数は300人中15人でした。境界がここまで近い項目は、その選択肢が回答者に区別されていない可能性があります。尺度を改訂するなら、隣のカテゴリとの統合候補になります。
あるカテゴリを誰も選ばなかった項目では、関連する境界の推定が不安定になり、事前分布への依存が強まることがあります。ただし、隣接カテゴリなどの観測も尤度を通じて境界を制約します。推定前に度数表を作り、度数0の原因を確認します。隣接カテゴリとの統合は選択肢の意味と尺度の比較可能性を確認して判断し、項目ごとにKが変わる場合は尤度の実装も変更します。
推定結果の読み方：識別力・境界困難度・項目情報量
識別力aと境界困難度bの解釈
識別力aは、θが変わったときに回答カテゴリがどれだけ鋭く切り替わるかを表します。今回の推定では項目0のaが1.59、項目4が0.82で、項目4はθの違いを回答に反映しにくい項目です。ただし前節のとおり区間が広いので、各項目の89%区間で不確実性を確認し、識別力の大小を判断するときは、同じ事後標本でa同士の差を計算して、その区間や差が正になる確率を確認します。
境界困難度bは、その境界を越える確率が0.5になるθの位置です。項目0は4本の境界が0.58〜2.20とθの高い側に寄っており、事後平均のa・bを代入すると、平均的な回答者（θ=0）が最低カテゴリを選ぶ確率は約71%です。観測度数で300人中200人がカテゴリ0だったのはこのためです。
項目情報量で測定精度の高いθ帯を特定する計算
GRMの項目情報量は I(θ)=Σk Pk′(θ)2/Pk(θ) で、局所独立を仮定すると全項目の和がテスト情報量となり、その逆数の平方根は、項目パラメータを既知とする場合のθの最尤推定の標準誤差の漸近的な目安です。本記事で推定したθの事後標準偏差とは異なります。事後平均のa・bを入れて計算するコードです。
a = idata.posterior["a"].mean(("chain", "draw")).values
b = idata.posterior["b"].mean(("chain", "draw")).values
grid = np.linspace(-3, 3, 13)
def info(aj, bj):
ps = 1 / (1 + np.exp(-aj * (grid[:, None] - bj[None]))) # P(Y>=k)
Ps = np.hstack([np.ones((len(grid), 1)), ps, np.zeros((len(grid), 1))])
dPs = aj * Ps * (1 - Ps)
P = Ps[:, :-1] - Ps[:, 1:]
dP = dPs[:, :-1] - dPs[:, 1:]
return (dP ** 2 / P).sum(1)
test_info = sum(info(a[j], b[j]) for j in range(len(a)))
se = 1 / np.sqrt(test_info)| θ | −3 | −2 | −1 | 0 | 0.5 | 1 | 2 | 3 |
|---|---|---|---|---|---|---|---|---|
| テスト情報量 | 0.73 | 1.53 | 2.51 | 3.18 | 3.36 | 3.33 | 2.62 | 1.39 |
| 標準誤差 | 1.17 | 0.81 | 0.63 | 0.56 | 0.55 | 0.55 | 0.62 | 0.85 |
この8項目の尺度は θ=0.5〜1 付近で最も精度が高く、θ=−3では標準誤差が1.17と、θの標準偏差（1）を超えます。特性の低い層を見分けたい尺度なら、境界が低い側にある項目を足す必要があります。項目別に見ると、項目0の情報量は θ=1.5 で最大、項目2と項目6は θ=−0.5 で最大でした。短縮版を作るときは、測りたいθ帯に情報量の山を持つ項目を残します。
よくある質問
GRMとは何の略で、どんなデータに使いますか？
Graded Response Model（段階反応モデル）の略で、Samejima（1969）が提案した項目反応理論のモデルです。「全く当てはまらない」〜「非常に当てはまる」のように選択肢に順序がある質問紙や、論述を段階で採点したデータに使います。選択肢が2つならGRMは2PLと一致します。
ベイズ推定は何を推定しているのですか？
パラメータの事前分布と、データが得られる確率（尤度）を掛け合わせた事後分布を推定しています。GRMでは事後分布を式で解けないため、PyMCはNUTS（MCMCの一種）で事後分布から標本を引き、その平均や区間を推定値にします。考え方の土台はベイズの定理とは？公式の意味・具体例からナイーブベイズ実装までで解説しています。
弁別力（識別力）とは何ですか？
θの違いを回答の違いとしてどれだけ鋭く反映するかを表すパラメータaです。aが大きい項目は、境界付近でθが少し変わるだけで選ぶカテゴリが切り替わります。GRMでは項目ごとに1つ推定し、本記事の推定では0.82〜1.59でした。
PyMCのNUTSサンプラーで発散（divergences）が出たらどうしますか？
まず target_accept を既定の0.8から0.9〜0.95へ上げます。それでも残る場合は、事前分布が広すぎないか、aを正に縛っているか、θの尺度を固定しているかを見直し、必要なら非中心化などの再パラメータ化をします。発散が1件でも残る結果は、そのまま解釈しないのが原則です。
境界ごとの曲線とカテゴリごとの曲線は何が違いますか？
境界ごとの曲線は「k以上を選ぶ確率」P(Y≥k|θ) で、5件法なら4本の右上がりのS字になります。カテゴリごとの曲線は「ちょうどkを選ぶ確率」P(Y=k|θ) で、5本のうち両端は単調、中間カテゴリは山型です。中間カテゴリがどのθでも最頻カテゴリにならない場合は、選択肢の機能を点検する材料になります。ただし、それだけで機能していないとは断定せず、選択肢の意味、観測度数、測定情報も確認します。
