20260313 DINA pymc - 🫶🏻 U1FAF6-project

In [ ]:

Copied!

pip install "nutpie[pymc]"
pip install "nutpie[pymc]"

Collecting nutpie[pymc]
  Downloading nutpie-0.16.8-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (6.5 kB)
Requirement already satisfied: pyarrow>=12.0.0 in /usr/local/lib/python3.12/dist-packages (from nutpie[pymc]) (18.1.0)
Collecting arro3-core>=0.6.0 (from nutpie[pymc])
  Downloading arro3_core-0.8.0-cp311-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (515 bytes)
Requirement already satisfied: pandas>=2.0 in /usr/local/lib/python3.12/dist-packages (from nutpie[pymc]) (2.2.2)
Requirement already satisfied: xarray>=2025.1.2 in /usr/local/lib/python3.12/dist-packages (from nutpie[pymc]) (2025.12.0)
Requirement already satisfied: arviz<1.0,>=0.20.0 in /usr/local/lib/python3.12/dist-packages (from nutpie[pymc]) (0.22.0)
Collecting obstore>=0.8.0 (from nutpie[pymc])
  Downloading obstore-0.9.2-cp311-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (3.9 kB)
Collecting zarr>=3.1.0 (from nutpie[pymc])
  Downloading zarr-3.1.5-py3-none-any.whl.metadata (10 kB)
Requirement already satisfied: pymc>=5.20.1 in /usr/local/lib/python3.12/dist-packages (from nutpie[pymc]) (5.28.1)
Requirement already satisfied: numba>=0.60.0 in /usr/local/lib/python3.12/dist-packages (from nutpie[pymc]) (0.60.0)
Requirement already satisfied: setuptools>=60.0.0 in /usr/local/lib/python3.12/dist-packages (from arviz<1.0,>=0.20.0->nutpie[pymc]) (75.2.0)
Requirement already satisfied: matplotlib>=3.8 in /usr/local/lib/python3.12/dist-packages (from arviz<1.0,>=0.20.0->nutpie[pymc]) (3.10.0)
Requirement already satisfied: numpy>=1.26.0 in /usr/local/lib/python3.12/dist-packages (from arviz<1.0,>=0.20.0->nutpie[pymc]) (2.0.2)
Requirement already satisfied: scipy>=1.11.0 in /usr/local/lib/python3.12/dist-packages (from arviz<1.0,>=0.20.0->nutpie[pymc]) (1.16.3)
Requirement already satisfied: packaging in /usr/local/lib/python3.12/dist-packages (from arviz<1.0,>=0.20.0->nutpie[pymc]) (26.0)
Requirement already satisfied: h5netcdf>=1.0.2 in /usr/local/lib/python3.12/dist-packages (from arviz<1.0,>=0.20.0->nutpie[pymc]) (1.8.1)
Requirement already satisfied: typing-extensions>=4.1.0 in /usr/local/lib/python3.12/dist-packages (from arviz<1.0,>=0.20.0->nutpie[pymc]) (4.15.0)
Requirement already satisfied: xarray-einstats>=0.3 in /usr/local/lib/python3.12/dist-packages (from arviz<1.0,>=0.20.0->nutpie[pymc]) (0.10.0)
Requirement already satisfied: llvmlite<0.44,>=0.43.0dev0 in /usr/local/lib/python3.12/dist-packages (from numba>=0.60.0->nutpie[pymc]) (0.43.0)
Requirement already satisfied: python-dateutil>=2.8.2 in /usr/local/lib/python3.12/dist-packages (from pandas>=2.0->nutpie[pymc]) (2.9.0.post0)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.12/dist-packages (from pandas>=2.0->nutpie[pymc]) (2025.2)
Requirement already satisfied: tzdata>=2022.7 in /usr/local/lib/python3.12/dist-packages (from pandas>=2.0->nutpie[pymc]) (2025.3)
Requirement already satisfied: cachetools<7,>=4.2.1 in /usr/local/lib/python3.12/dist-packages (from pymc>=5.20.1->nutpie[pymc]) (6.2.6)
Requirement already satisfied: cloudpickle in /usr/local/lib/python3.12/dist-packages (from pymc>=5.20.1->nutpie[pymc]) (3.1.2)
Requirement already satisfied: pytensor<2.39,>=2.38.0 in /usr/local/lib/python3.12/dist-packages (from pymc>=5.20.1->nutpie[pymc]) (2.38.2)
Requirement already satisfied: rich>=13.7.1 in /usr/local/lib/python3.12/dist-packages (from pymc>=5.20.1->nutpie[pymc]) (13.9.4)
Requirement already satisfied: threadpoolctl<4.0.0,>=3.1.0 in /usr/local/lib/python3.12/dist-packages (from pymc>=5.20.1->nutpie[pymc]) (3.6.0)
Collecting donfig>=0.8 (from zarr>=3.1.0->nutpie[pymc])
  Downloading donfig-0.8.1.post1-py3-none-any.whl.metadata (5.0 kB)
Requirement already satisfied: google-crc32c>=1.5 in /usr/local/lib/python3.12/dist-packages (from zarr>=3.1.0->nutpie[pymc]) (1.8.0)
Collecting numcodecs>=0.14 (from zarr>=3.1.0->nutpie[pymc])
  Downloading numcodecs-0.16.5-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl.metadata (3.4 kB)
Requirement already satisfied: pyyaml in /usr/local/lib/python3.12/dist-packages (from donfig>=0.8->zarr>=3.1.0->nutpie[pymc]) (6.0.3)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.8->arviz<1.0,>=0.20.0->nutpie[pymc]) (1.3.3)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.8->arviz<1.0,>=0.20.0->nutpie[pymc]) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.8->arviz<1.0,>=0.20.0->nutpie[pymc]) (4.62.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.8->arviz<1.0,>=0.20.0->nutpie[pymc]) (1.5.0)
Requirement already satisfied: pillow>=8 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.8->arviz<1.0,>=0.20.0->nutpie[pymc]) (11.3.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.8->arviz<1.0,>=0.20.0->nutpie[pymc]) (3.3.2)
Requirement already satisfied: filelock>=3.15 in /usr/local/lib/python3.12/dist-packages (from pytensor<2.39,>=2.38.0->pymc>=5.20.1->nutpie[pymc]) (3.25.1)
Requirement already satisfied: etuples in /usr/local/lib/python3.12/dist-packages (from pytensor<2.39,>=2.38.0->pymc>=5.20.1->nutpie[pymc]) (0.3.10)
Requirement already satisfied: logical-unification in /usr/local/lib/python3.12/dist-packages (from pytensor<2.39,>=2.38.0->pymc>=5.20.1->nutpie[pymc]) (0.4.7)
Requirement already satisfied: miniKanren in /usr/local/lib/python3.12/dist-packages (from pytensor<2.39,>=2.38.0->pymc>=5.20.1->nutpie[pymc]) (1.0.5)
Requirement already satisfied: cons in /usr/local/lib/python3.12/dist-packages (from pytensor<2.39,>=2.38.0->pymc>=5.20.1->nutpie[pymc]) (0.4.7)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.12/dist-packages (from python-dateutil>=2.8.2->pandas>=2.0->nutpie[pymc]) (1.17.0)
Requirement already satisfied: markdown-it-py>=2.2.0 in /usr/local/lib/python3.12/dist-packages (from rich>=13.7.1->pymc>=5.20.1->nutpie[pymc]) (4.0.0)
Requirement already satisfied: pygments<3.0.0,>=2.13.0 in /usr/local/lib/python3.12/dist-packages (from rich>=13.7.1->pymc>=5.20.1->nutpie[pymc]) (2.19.2)
Requirement already satisfied: mdurl~=0.1 in /usr/local/lib/python3.12/dist-packages (from markdown-it-py>=2.2.0->rich>=13.7.1->pymc>=5.20.1->nutpie[pymc]) (0.1.2)
Requirement already satisfied: toolz in /usr/local/lib/python3.12/dist-packages (from logical-unification->pytensor<2.39,>=2.38.0->pymc>=5.20.1->nutpie[pymc]) (0.12.1)
Requirement already satisfied: multipledispatch in /usr/local/lib/python3.12/dist-packages (from logical-unification->pytensor<2.39,>=2.38.0->pymc>=5.20.1->nutpie[pymc]) (1.0.0)
Downloading arro3_core-0.8.0-cp311-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (3.0 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.0/3.0 MB 31.6 MB/s eta 0:00:00
Downloading obstore-0.9.2-cp311-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (4.2 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 4.2/4.2 MB 80.0 MB/s eta 0:00:00
Downloading zarr-3.1.5-py3-none-any.whl (284 kB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 284.1/284.1 kB 9.7 MB/s eta 0:00:00
Downloading nutpie-0.16.8-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (8.7 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 8.7/8.7 MB 64.5 MB/s eta 0:00:00
Downloading donfig-0.8.1.post1-py3-none-any.whl (21 kB)
Downloading numcodecs-0.16.5-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl (9.2 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 9.2/9.2 MB 36.3 MB/s eta 0:00:00
Installing collected packages: obstore, numcodecs, donfig, arro3-core, zarr, nutpie
Successfully installed arro3-core-0.8.0 donfig-0.8.1.post1 numcodecs-0.16.5 nutpie-0.16.8 obstore-0.9.2 zarr-3.1.5

DINA model in PyMC¶

2026-03-13
jw

本文先退而以 DINA model 實做。

DINA (Deterministic Input, Noisy AND gate). 假設該題目需要所對應的所有屬性都掌握 (AND)，才能答對。（屬性的乘法關係，交互效果部份。）適合認知測驗。
DINO (Deterministic Input, Noisy OR gate). 假設該題目需所對應的屬性掌握一種，就能答對 (OR)。（屬性的加法關係，主效果部份。）適合非認知測驗。
GDINA. 納入所有可能的主效果和交互效果的飽和模式。讓資料決定哪些效果重要，不預設結構。

考量點。¶

GDINA 雖然比較彈性，但每個題目的參數量也較多。因此考量到參數數量和樣本數量的權衡，要使用飽和模式 (GDINA, LCDM, ...) 還是縮減模式(DINA, DINO, ...) 需要研究者選擇。以 GDINA 而言，每個題目需要 $2^{K_j}$ 個參數，${K_j}$ 是看某題目 $j$ 對應幾個屬性 $K$. 而 DINA 模式每個題目就對應2個參數 (slip, guess)。

以 20 題，5 個屬性為例，

DINA. $20 \times 2 = 40$ 個參數。
GDINA. 最低情況每題只對應 1 個屬性（但這會有識別性問題），$20 \times 2^1 = 40$ 個參數，最高情況每題對應 5 個屬性（這也有識別性問題），$20 \times 2^5 = 640$ 個參數。正常情況應該介於 40--640 個參數之間。仍然是非常大的參數量。在 Q-matrix 設計得當的情況下，也需要有足夠的樣本才能估計得好。

另外， GDINA 的估計方法也更為複雜。DINA 是相對容易且穩定的選項之一。

DCM 與 mixture modeling。¶

診斷分類模式 (diagnostic classification model, DCM) (這邊所說的 DINA, DINO, GDINA 都屬於此類) 可視為一種「驗證性」（或說限制性）的分類模式，也就是說要分幾類是確定好的。即 $2^K$ 類，K 是屬性數量。例如 3 個屬性，就有 $2^3=8$ 個精熟組型，從 {0,0,0}, {0,0,1},...,{1,1,1}, 只是有些組可能沒有人。詳見本分文件的最後。

相較於「探索性」的 latent class analysis (LCA)(可視為固定潛在變異數為 0 的 GMM), 或是 k-means (用距離決定的) ，或其他各式 mixture model（例如 Gaussian mixture model），要幾類往往是「探索」出來的。也就是研究者透過模式比較，決定要分成幾類。

因此，DCM 有別於其他的分類模式。一般的分類模式只是單純根據變項把人做分類，但 DCM 更重要的目標是把人分進「設定好的精熟類型」中，藉此知道每位學生在這些屬性（或技能等等）是否精熟/不精熟。達到認知診斷的效果。

MCMC¶

也可以使用 EM 估計。但需要另外檢驗。

在貝氏的架構下可以直覺的寫出統計模式，並透過自動抽樣的方式得到參數估計。因為 DCM 涉及到離散的潛在變數，因此過去通常用 Gibbs sampler 估計（例如 jags, 見下方文獻）

但考量到在 python 上實做，比較成熟的貝氏統計軟體為 PyMC . 但 PyMC 是用 HMC/NUTS 抽樣器，該方法比較快，但對於離散潛在變數不友善。因此需要特殊處理（logsumexp）（見下）。

Ref¶

(Bayesian-Gibbs)

Zhan, P., Jiao, H., Man, K., & Wang, L. (2019). Using JAGS for Bayesian cognitive diagnosis modeling: A tutorial. Journal of Educational and Behavioral Statistics, 44(4), 473-503. LINK

(EM-DINA)

De La Torre, J. (2009). DINA model and parameter estimation: A didactic. Journal of educational and behavioral statistics, 34(1), 115-130. LINK

(EM-GDINA)

De La Torre, J. (2011). The generalized DINA model framework. Psychometrika, 76(2), 179-199. LINK

In [ ]:

Copied!





"""
DINA + PyMC 貝氏估計
========================
模型：Deterministic Input, Noisy And Gate (DINA)
  - 每題兩個參數：slip (s_j) 與 guess (g_j)
  - 理想反應：η_ij = ∏_k α_ik^q_jk  （必須掌握所有必要屬性）
  - P(Y_ij=1 | α_i) = (1 - s_j)^η_ij * g_j^(1 - η_ij)
模擬資料：500人、20題、3屬性
"""

import numpy as np
import pandas as pd
import itertools
import warnings
warnings.filterwarnings("ignore")

import pymc as pm
import pytensor.tensor as pt
from scipy.special import logsumexp
from scipy.stats import pearsonr

import nutpie
"""
DINA + PyMC 貝氏估計
========================
模型：Deterministic Input, Noisy And Gate (DINA)
  - 每題兩個參數：slip (s_j) 與 guess (g_j)
  - 理想反應：η_ij = ∏_k α_ik^q_jk  （必須掌握所有必要屬性）
  - P(Y_ij=1 | α_i) = (1 - s_j)^η_ij * g_j^(1 - η_ij)
模擬資料：500人、20題、3屬性
"""

import numpy as np
import pandas as pd
import itertools
import warnings
warnings.filterwarnings("ignore")

import pymc as pm
import pytensor.tensor as pt
from scipy.special import logsumexp
from scipy.stats import pearsonr

import nutpie

In [ ]:

Copied!





# ══════════════════════════════════════════════════════════════
# PART 1：模擬資料
# ══════════════════════════════════════════════════════════════

np.random.seed(1234)

N, J, K = 1500, 20, 3
all_alpha = np.array(list(itertools.product([0, 1], repeat=K)))  # (8, 3)
L = len(all_alpha)  # 2^K = 8 種屬性組合

# Q 矩陣：指定每題測量哪些屬性
Q = np.array([
    [1, 0, 0], [1, 0, 0], [1, 0, 0], [1, 0, 0],   # 題 1-4：僅測屬性1
    [0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0],   # 題 5-8：僅測屬性2
    [0, 0, 1], [0, 0, 1], [0, 0, 1], [0, 0, 1],   # 題 9-12：僅測屬性3
    [1, 1, 0], [1, 1, 0], [1, 1, 0],               # 題 13-15：測屬性1+2
    [1, 0, 1], [1, 0, 1], [1, 0, 1],               # 題 16-18：測屬性1+3
    [1, 1, 1], [1, 1, 1],                           # 題 19-20：測全部屬性
])

# 屬性組合的先驗比例
pi_true = np.array([0.05, 0.08, 0.08, 0.12, 0.10, 0.12, 0.15, 0.30])

# DINA 真實參數：slip 與 guess
# slip：掌握所有屬性但答錯的概率；guess：未掌握屬性但答對的概率
slip_true  = np.array([0.10, 0.15, 0.10, 0.15,   # 題 1-4
                        0.10, 0.15, 0.10, 0.15,   # 題 5-8
                        0.10, 0.15, 0.10, 0.15,   # 題 9-12
                        0.10, 0.15, 0.10,           # 題 13-15
                        0.10, 0.15, 0.10,           # 題 16-18
                        0.10, 0.15])               # 題 19-20

guess_true = np.array([0.15, 0.20, 0.10, 0.25,   # 題 1-4
                        0.15, 0.20, 0.10, 0.25,   # 題 5-8
                        0.15, 0.20, 0.10, 0.25,   # 題 9-12
                        0.10, 0.15, 0.10,           # 題 13-15
                        0.10, 0.15, 0.10,           # 題 16-18
                        0.10, 0.15])               # 題 19-20

# 計算各屬性組合下，每題的理想反應 η_lj = ∏_k α_lk^q_jk
# eta_true[l, j] = 1 iff 屬性組合 l 掌握題 j 的所有必要屬性
eta_true = np.array([
    [int(np.all(all_alpha[l] >= Q[j])) for j in range(J)]
    for l in range(L)
])  # (L, J)

# DINA 成功概率：P(Y=1|α_l) = (1-s_j)^η_lj * g_j^(1-η_lj)
P_true = np.zeros((J, L))
for j in range(J):
    for l in range(L):
        if eta_true[l, j] == 1:
            P_true[j, l] = 1 - slip_true[j]    # 掌握 → 成功概率高
        else:
            P_true[j, l] = guess_true[j]        # 未掌握 → 猜對概率低

# 模擬受試者屬性組合
alpha_idx  = np.random.choice(L, size=N, p=pi_true)
Alpha_true = all_alpha[alpha_idx]

# 模擬作答資料
Y = np.array([
    [np.random.binomial(1, P_true[j, alpha_idx[i]]) for j in range(J)]
    for i in range(N)
])

print(f"模擬資料完成：Y={Y.shape}, Q={Q.shape}")
print(f"slip 範圍：{slip_true.min():.2f} ~ {slip_true.max():.2f}")
print(f"guess 範圍：{guess_true.min():.2f} ~ {guess_true.max():.2f}")
# ══════════════════════════════════════════════════════════════
# PART 1：模擬資料
# ══════════════════════════════════════════════════════════════

np.random.seed(1234)

N, J, K = 1500, 20, 3
all_alpha = np.array(list(itertools.product([0, 1], repeat=K)))  # (8, 3)
L = len(all_alpha)  # 2^K = 8 種屬性組合

# Q 矩陣：指定每題測量哪些屬性
Q = np.array([
    [1, 0, 0], [1, 0, 0], [1, 0, 0], [1, 0, 0],   # 題 1-4：僅測屬性1
    [0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0],   # 題 5-8：僅測屬性2
    [0, 0, 1], [0, 0, 1], [0, 0, 1], [0, 0, 1],   # 題 9-12：僅測屬性3
    [1, 1, 0], [1, 1, 0], [1, 1, 0],               # 題 13-15：測屬性1+2
    [1, 0, 1], [1, 0, 1], [1, 0, 1],               # 題 16-18：測屬性1+3
    [1, 1, 1], [1, 1, 1],                           # 題 19-20：測全部屬性
])

# 屬性組合的先驗比例
pi_true = np.array([0.05, 0.08, 0.08, 0.12, 0.10, 0.12, 0.15, 0.30])

# DINA 真實參數：slip 與 guess
# slip：掌握所有屬性但答錯的概率；guess：未掌握屬性但答對的概率
slip_true  = np.array([0.10, 0.15, 0.10, 0.15,   # 題 1-4
                        0.10, 0.15, 0.10, 0.15,   # 題 5-8
                        0.10, 0.15, 0.10, 0.15,   # 題 9-12
                        0.10, 0.15, 0.10,           # 題 13-15
                        0.10, 0.15, 0.10,           # 題 16-18
                        0.10, 0.15])               # 題 19-20

guess_true = np.array([0.15, 0.20, 0.10, 0.25,   # 題 1-4
                        0.15, 0.20, 0.10, 0.25,   # 題 5-8
                        0.15, 0.20, 0.10, 0.25,   # 題 9-12
                        0.10, 0.15, 0.10,           # 題 13-15
                        0.10, 0.15, 0.10,           # 題 16-18
                        0.10, 0.15])               # 題 19-20

# 計算各屬性組合下，每題的理想反應 η_lj = ∏_k α_lk^q_jk
# eta_true[l, j] = 1 iff 屬性組合 l 掌握題 j 的所有必要屬性
eta_true = np.array([
    [int(np.all(all_alpha[l] >= Q[j])) for j in range(J)]
    for l in range(L)
])  # (L, J)

# DINA 成功概率：P(Y=1|α_l) = (1-s_j)^η_lj * g_j^(1-η_lj)
P_true = np.zeros((J, L))
for j in range(J):
    for l in range(L):
        if eta_true[l, j] == 1:
            P_true[j, l] = 1 - slip_true[j]    # 掌握 → 成功概率高
        else:
            P_true[j, l] = guess_true[j]        # 未掌握 → 猜對概率低

# 模擬受試者屬性組合
alpha_idx  = np.random.choice(L, size=N, p=pi_true)
Alpha_true = all_alpha[alpha_idx]

# 模擬作答資料
Y = np.array([
    [np.random.binomial(1, P_true[j, alpha_idx[i]]) for j in range(J)]
    for i in range(N)
])

print(f"模擬資料完成：Y={Y.shape}, Q={Q.shape}")
print(f"slip 範圍：{slip_true.min():.2f} ~ {slip_true.max():.2f}")
print(f"guess 範圍：{guess_true.min():.2f} ~ {guess_true.max():.2f}")

模擬資料完成：Y=(1500, 20), Q=(20, 3)
slip 範圍：0.10 ~ 0.15
guess 範圍：0.10 ~ 0.25

In [ ]:

Copied!





# ══════════════════════════════════════════════════════════════
# PART 2：預計算 η 矩陣（PyMC 外部）
# ══════════════════════════════════════════════════════════════
# η_lj = 1 iff 屬性組合 l 掌握了題 j 的所有必要屬性
# 此為固定的結構矩陣，由 Q 矩陣決定，不需要估計

# eta[l, j]: 屬性組合 l 對題 j 的理想反應
eta = np.array([
    [int(np.all(all_alpha[l] >= Q[j])) for j in range(J)]
    for l in range(L)
], dtype=np.float32)  # (L, J)

print("η 矩陣（前 4 題 × 全部 8 個屬性組合）：")
df_eta = pd.DataFrame(eta[:, :4],
                       index=["".join(map(str, a)) for a in all_alpha],
                       columns=[f"Q{j+1}" for j in range(4)])
print(df_eta.to_string())
print("\n（η=1 表示該屬性組合掌握了題目所有必要屬性）")
# ══════════════════════════════════════════════════════════════
# PART 2：預計算 η 矩陣（PyMC 外部）
# ══════════════════════════════════════════════════════════════
# η_lj = 1 iff 屬性組合 l 掌握了題 j 的所有必要屬性
# 此為固定的結構矩陣，由 Q 矩陣決定，不需要估計

# eta[l, j]: 屬性組合 l 對題 j 的理想反應
eta = np.array([
    [int(np.all(all_alpha[l] >= Q[j])) for j in range(J)]
    for l in range(L)
], dtype=np.float32)  # (L, J)

print("η 矩陣（前 4 題 × 全部 8 個屬性組合）：")
df_eta = pd.DataFrame(eta[:, :4],
                       index=["".join(map(str, a)) for a in all_alpha],
                       columns=[f"Q{j+1}" for j in range(4)])
print(df_eta.to_string())
print("\n（η=1 表示該屬性組合掌握了題目所有必要屬性）")

η 矩陣（前 4 題 × 全部 8 個屬性組合）：
      Q1   Q2   Q3   Q4
000  0.0  0.0  0.0  0.0
001  0.0  0.0  0.0  0.0
010  0.0  0.0  0.0  0.0
011  0.0  0.0  0.0  0.0
100  1.0  1.0  1.0  1.0
101  1.0  1.0  1.0  1.0
110  1.0  1.0  1.0  1.0
111  1.0  1.0  1.0  1.0

（η=1 表示該屬性組合掌握了題目所有必要屬性）

這邊用 VI 估計比較快，但是參數回復性比較不好。用 MCMC 估計比較慢（在 colab 上約須要 1-3 分鐘，到本地上可能可加速），但參數回復性較好。

In [ ]:

Copied!





# ══════════════════════════════════════════════════════════════
# PART 3：PyMC DINA（MCMC 估計）
# ══════════════════════════════════════════════════════════════
#
# DINA 模型：
#   P(Y_ij=1 | α_i) = (1 - s_j)^η_ij * g_j^(1 - η_ij)
#
# 先驗設定（參考論文 de la Torre 2009; Junker & Sijtsma 2001）：
#   slip_j ~ Beta(1, 4)   → 偏向小值（slip 通常較小）
#   guess_j ~ Beta(1, 4)  → 偏向小值（guess 通常較小）
#   π ~ Dirichlet(1, ..., 1)  → 均勻先驗
#
# 可識別性約束：s_j + g_j < 1  （透過先驗形狀隱式控制）

eta_pt = pt.constant(eta)  # (L, J)
Y_pt   = pt.constant(Y.astype(np.float32))  # (N, J)

with pm.Model() as dina_model:

    # ── 先驗 ────────────────────────────────────────────────
    # 屬性組合的混合比例
    pi = pm.Dirichlet("pi", a=np.ones(L))  # (L,)

    # 每題的 slip 與 guess 參數
    # Beta(1,4) 先驗讓參數偏向 0，符合 DINA 的識別條件
    slip  = pm.Beta("slip",  alpha=1, beta=4, shape=(J,))  # (J,)
    guess = pm.Beta("guess", alpha=1, beta=4, shape=(J,))  # (J,)

    # ── DINA 成功概率矩陣 P_mat (J, L) ─────────────────────
    # P_mat[j, l] = (1 - slip[j])^η[l,j] * guess[j]^(1 - η[l,j])
    # 利用 log-sum trick 計算
    #   log P = η[l,j] * log(1-slip[j]) + (1-η[l,j]) * log(guess[j])

    log_1_minus_slip = pt.log(1 - slip + 1e-8)  # (J,)
    log_guess        = pt.log(guess + 1e-8)      # (J,)

    # eta_pt: (L, J) → 轉置成 (J, L) 以便廣播
    eta_T = eta_pt.T  # (J, L)

    # log P_mat[j, l] = η[l,j]*log(1-s_j) + (1-η[l,j])*log(g_j)
    log_P_mat = (
        eta_T * log_1_minus_slip[:, None] +
        (1 - eta_T) * log_guess[:, None]
    )  # (J, L)

    P_mat = pt.exp(log_P_mat)  # (J, L)

    # ── Mixture likelihood（對 L 個屬性組合做 marginalization）─
    # log p(Y_i | π, s, g) = log Σ_l π_l * ∏_j P_mat[j,l]^Y_ij * (1-P_mat[j,l])^(1-Y_ij)

    log_likes = pt.stack([
        pt.log(pi[l] + 1e-8) +
        pt.dot(Y_pt,       pt.log(P_mat[:, l] + 1e-8)) +
        pt.dot(1 - Y_pt,   pt.log(1 - P_mat[:, l] + 1e-8))
        for l in range(L)
    ], axis=1)  # (N, L)

    pm.Potential("obs", pt.logsumexp(log_likes, axis=1).sum())

    # ── VI 估計（ADVI）──────────────────────────────────────
    """
    approx = pm.fit(20000,
                    #method="advi",
                    method="fullrank_advi",
                    callbacks=[pm.callbacks.CheckParametersConvergence(tolerance=1e-4)])

    trace = approx.sample(2000)
    """
    #trace = pm.sample(cores=2, target_accept=0.85, init="advi+adapt_diag")


    # -- nutpie

    compiled = nutpie.compile_pymc_model(dina_model)
    trace = nutpie.sample(compiled, chains=2, cores=2, target_accept=0.85)
# ══════════════════════════════════════════════════════════════
# PART 3：PyMC DINA（MCMC 估計）
# ══════════════════════════════════════════════════════════════
#
# DINA 模型：
#   P(Y_ij=1 | α_i) = (1 - s_j)^η_ij * g_j^(1 - η_ij)
#
# 先驗設定（參考論文 de la Torre 2009; Junker & Sijtsma 2001）：
#   slip_j ~ Beta(1, 4)   → 偏向小值（slip 通常較小）
#   guess_j ~ Beta(1, 4)  → 偏向小值（guess 通常較小）
#   π ~ Dirichlet(1, ..., 1)  → 均勻先驗
#
# 可識別性約束：s_j + g_j < 1  （透過先驗形狀隱式控制）

eta_pt = pt.constant(eta)  # (L, J)
Y_pt   = pt.constant(Y.astype(np.float32))  # (N, J)

with pm.Model() as dina_model:

    # ── 先驗 ────────────────────────────────────────────────
    # 屬性組合的混合比例
    pi = pm.Dirichlet("pi", a=np.ones(L))  # (L,)

    # 每題的 slip 與 guess 參數
    # Beta(1,4) 先驗讓參數偏向 0，符合 DINA 的識別條件
    slip  = pm.Beta("slip",  alpha=1, beta=4, shape=(J,))  # (J,)
    guess = pm.Beta("guess", alpha=1, beta=4, shape=(J,))  # (J,)

    # ── DINA 成功概率矩陣 P_mat (J, L) ─────────────────────
    # P_mat[j, l] = (1 - slip[j])^η[l,j] * guess[j]^(1 - η[l,j])
    # 利用 log-sum trick 計算
    #   log P = η[l,j] * log(1-slip[j]) + (1-η[l,j]) * log(guess[j])

    log_1_minus_slip = pt.log(1 - slip + 1e-8)  # (J,)
    log_guess        = pt.log(guess + 1e-8)      # (J,)

    # eta_pt: (L, J) → 轉置成 (J, L) 以便廣播
    eta_T = eta_pt.T  # (J, L)

    # log P_mat[j, l] = η[l,j]*log(1-s_j) + (1-η[l,j])*log(g_j)
    log_P_mat = (
        eta_T * log_1_minus_slip[:, None] +
        (1 - eta_T) * log_guess[:, None]
    )  # (J, L)

    P_mat = pt.exp(log_P_mat)  # (J, L)

    # ── Mixture likelihood（對 L 個屬性組合做 marginalization）─
    # log p(Y_i | π, s, g) = log Σ_l π_l * ∏_j P_mat[j,l]^Y_ij * (1-P_mat[j,l])^(1-Y_ij)

    log_likes = pt.stack([
        pt.log(pi[l] + 1e-8) +
        pt.dot(Y_pt,       pt.log(P_mat[:, l] + 1e-8)) +
        pt.dot(1 - Y_pt,   pt.log(1 - P_mat[:, l] + 1e-8))
        for l in range(L)
    ], axis=1)  # (N, L)

    pm.Potential("obs", pt.logsumexp(log_likes, axis=1).sum())

    # ── VI 估計（ADVI）──────────────────────────────────────
    """
    approx = pm.fit(20000,
                    #method="advi",
                    method="fullrank_advi",
                    callbacks=[pm.callbacks.CheckParametersConvergence(tolerance=1e-4)])

    trace = approx.sample(2000)
    """
    #trace = pm.sample(cores=2, target_accept=0.85, init="advi+adapt_diag")


    # -- nutpie

    compiled = nutpie.compile_pymc_model(dina_model)
    trace = nutpie.sample(compiled, chains=2, cores=2, target_accept=0.85)

Sampler Progress

Total Chains: 2

Active Chains: 0

Finished Chains: 2

Sampling for a minute

Estimated Time to Completion: now

Progress	Draws	Divergences	Step Size	Gradients/Draw
	1400	0	0.44	7
	1400	0	0.35	15

In [ ]:

Copied!





# ══════════════════════════════════════════════════════════════
# PART 4：基本參數回收結果
# ══════════════════════════════════════════════════════════════

# 提取後驗均值
slip_est  = trace.posterior["slip"].values.reshape(-1, J).mean(axis=0)   # (J,)
guess_est = trace.posterior["guess"].values.reshape(-1, J).mean(axis=0)  # (J,)
pi_est    = trace.posterior["pi"].values.reshape(-1, L).mean(axis=0)     # (L,)

# Pearson r
r_slip,  _ = pearsonr(slip_true,  slip_est)
r_guess, _ = pearsonr(guess_true, guess_est)
r_pi,    _ = pearsonr(pi_true,    pi_est)

print("參數回收（Pearson r）：")
print(f"  slip   r = {r_slip:.4f}")
print(f"  guess  r = {r_guess:.4f}")
print(f"  π      r = {r_pi:.4f}")
# ══════════════════════════════════════════════════════════════
# PART 4：基本參數回收結果
# ══════════════════════════════════════════════════════════════

# 提取後驗均值
slip_est  = trace.posterior["slip"].values.reshape(-1, J).mean(axis=0)   # (J,)
guess_est = trace.posterior["guess"].values.reshape(-1, J).mean(axis=0)  # (J,)
pi_est    = trace.posterior["pi"].values.reshape(-1, L).mean(axis=0)     # (L,)

# Pearson r
r_slip,  _ = pearsonr(slip_true,  slip_est)
r_guess, _ = pearsonr(guess_true, guess_est)
r_pi,    _ = pearsonr(pi_true,    pi_est)

print("參數回收（Pearson r）：")
print(f"  slip   r = {r_slip:.4f}")
print(f"  guess  r = {r_guess:.4f}")
print(f"  π      r = {r_pi:.4f}")

參數回收（Pearson r）：
  slip   r = 0.8583
  guess  r = 0.9777
  π      r = 0.9914

In [ ]:

Copied!





# ══════════════════════════════════════════════════════════════
# PART 5：各參數估計表
# ══════════════════════════════════════════════════════════════

profile_labels = ["".join(map(str, a)) for a in all_alpha]

# ── 5a. slip 與 guess：真實 vs 估計 ─────────────────────────
slip_sd  = trace.posterior["slip"].values.reshape(-1, J).std(axis=0)
guess_sd = trace.posterior["guess"].values.reshape(-1, J).std(axis=0)

df_sg = pd.DataFrame({
    "題目":       [f"Q{j+1}" for j in range(J)],
    "slip 真實":  slip_true.round(3),
    "slip 估計":  slip_est.round(3),
    "slip SD":    slip_sd.round(3),
    "slip 偏差":  (slip_est - slip_true).round(3),
    "guess 真實": guess_true.round(3),
    "guess 估計": guess_est.round(3),
    "guess SD":   guess_sd.round(3),
    "guess 偏差": (guess_est - guess_true).round(3),
})
print("\n" + "="*70)
print("【DINA 參數：slip & guess — 真實 vs 估計】")
print("="*70)
print(df_sg.to_string(index=False))

# 平均絕對偏差
mab_slip  = np.abs(slip_est  - slip_true).mean()
mab_guess = np.abs(guess_est - guess_true).mean()
print(f"\n  slip  MAB = {mab_slip:.4f}")
print(f"  guess MAB = {mab_guess:.4f}")
# ══════════════════════════════════════════════════════════════
# PART 5：各參數估計表
# ══════════════════════════════════════════════════════════════

profile_labels = ["".join(map(str, a)) for a in all_alpha]

# ── 5a. slip 與 guess：真實 vs 估計 ─────────────────────────
slip_sd  = trace.posterior["slip"].values.reshape(-1, J).std(axis=0)
guess_sd = trace.posterior["guess"].values.reshape(-1, J).std(axis=0)

df_sg = pd.DataFrame({
    "題目":       [f"Q{j+1}" for j in range(J)],
    "slip 真實":  slip_true.round(3),
    "slip 估計":  slip_est.round(3),
    "slip SD":    slip_sd.round(3),
    "slip 偏差":  (slip_est - slip_true).round(3),
    "guess 真實": guess_true.round(3),
    "guess 估計": guess_est.round(3),
    "guess SD":   guess_sd.round(3),
    "guess 偏差": (guess_est - guess_true).round(3),
})
print("\n" + "="*70)
print("【DINA 參數：slip & guess — 真實 vs 估計】")
print("="*70)
print(df_sg.to_string(index=False))

# 平均絕對偏差
mab_slip  = np.abs(slip_est  - slip_true).mean()
mab_guess = np.abs(guess_est - guess_true).mean()
print(f"\n  slip  MAB = {mab_slip:.4f}")
print(f"  guess MAB = {mab_guess:.4f}")

======================================================================
【DINA 參數：slip & guess — 真實 vs 估計】
======================================================================
 題目  slip 真實  slip 估計  slip SD  slip 偏差  guess 真實  guess 估計  guess SD  guess 偏差
 Q1     0.10    0.110    0.010    0.010      0.15     0.167     0.017     0.017
 Q2     0.15    0.170    0.012    0.020      0.20     0.205     0.019     0.005
 Q3     0.10    0.112    0.010    0.012      0.10     0.128     0.015     0.028
 Q4     0.15    0.147    0.011   -0.003      0.25     0.254     0.020     0.004
 Q5     0.10    0.100    0.010   -0.000      0.15     0.165     0.018     0.015
 Q6     0.15    0.152    0.012    0.002      0.20     0.192     0.018    -0.008
 Q7     0.10    0.101    0.010    0.001      0.10     0.120     0.016     0.020
 Q8     0.15    0.144    0.012   -0.006      0.25     0.248     0.020    -0.002
 Q9     0.10    0.105    0.011    0.005      0.15     0.147     0.016    -0.003
Q10     0.15    0.140    0.012   -0.010      0.20     0.199     0.017    -0.001
Q11     0.10    0.088    0.010   -0.012      0.10     0.105     0.013     0.005
Q12     0.15    0.141    0.011   -0.009      0.25     0.258     0.019     0.008
Q13     0.10    0.105    0.012    0.005      0.10     0.100     0.011     0.000
Q14     0.15    0.127    0.013   -0.023      0.15     0.135     0.012    -0.015
Q15     0.10    0.091    0.011   -0.009      0.10     0.112     0.012     0.012
Q16     0.10    0.118    0.013    0.018      0.10     0.087     0.010    -0.013
Q17     0.15    0.141    0.014   -0.009      0.15     0.151     0.012     0.001
Q18     0.10    0.080    0.010   -0.020      0.10     0.100     0.010     0.000
Q19     0.10    0.131    0.015    0.031      0.10     0.087     0.009    -0.013
Q20     0.15    0.166    0.017    0.016      0.15     0.162     0.011     0.012

  slip  MAB = 0.0109
  guess MAB = 0.0092

In [ ]:

Copied!





# ── 5b. 由 slip/guess 重建 P 矩陣 ───────────────────────────
# P_est[j, l] = (1 - slip_est[j])^η[l,j] * guess_est[j]^(1 - η[l,j])

P_est = np.zeros((J, L))
for j in range(J):
    for l in range(L):
        if eta[l, j] == 1:
            P_est[j, l] = 1 - slip_est[j]
        else:
            P_est[j, l] = guess_est[j]

df_P_true = pd.DataFrame(P_true,
                          index=[f"Q{j+1}" for j in range(J)],
                          columns=profile_labels)
df_P_est  = pd.DataFrame(P_est,
                          index=[f"Q{j+1}" for j in range(J)],
                          columns=profile_labels)

print("\n" + "="*60)
print("【各題成功概率 P — 真實值】")
print("  （每題只有 2 種不同值：guess vs 1-slip）")
print("="*60)
print(df_P_true.round(3).to_string())

print("\n" + "="*60)
print("【各題成功概率 P — PyMC DINA 估計值（後驗均值）】")
print("="*60)
print(df_P_est.round(3).to_string())

df_P_diff = df_P_est - df_P_true
mab = df_P_diff.abs().values.mean()
rmse = np.sqrt((df_P_diff.values**2).mean())
print(f"\n  P 矩陣 MAB  = {mab:.4f}")
print(f"  P 矩陣 RMSE = {rmse:.4f}")
# ── 5b. 由 slip/guess 重建 P 矩陣 ───────────────────────────
# P_est[j, l] = (1 - slip_est[j])^η[l,j] * guess_est[j]^(1 - η[l,j])

P_est = np.zeros((J, L))
for j in range(J):
    for l in range(L):
        if eta[l, j] == 1:
            P_est[j, l] = 1 - slip_est[j]
        else:
            P_est[j, l] = guess_est[j]

df_P_true = pd.DataFrame(P_true,
                          index=[f"Q{j+1}" for j in range(J)],
                          columns=profile_labels)
df_P_est  = pd.DataFrame(P_est,
                          index=[f"Q{j+1}" for j in range(J)],
                          columns=profile_labels)

print("\n" + "="*60)
print("【各題成功概率 P — 真實值】")
print("  （每題只有 2 種不同值：guess vs 1-slip）")
print("="*60)
print(df_P_true.round(3).to_string())

print("\n" + "="*60)
print("【各題成功概率 P — PyMC DINA 估計值（後驗均值）】")
print("="*60)
print(df_P_est.round(3).to_string())

df_P_diff = df_P_est - df_P_true
mab = df_P_diff.abs().values.mean()
rmse = np.sqrt((df_P_diff.values**2).mean())
print(f"\n  P 矩陣 MAB  = {mab:.4f}")
print(f"  P 矩陣 RMSE = {rmse:.4f}")

============================================================
【各題成功概率 P — 真實值】
  （每題只有 2 種不同值：guess vs 1-slip）
============================================================
      000   001   010   011   100   101   110   111
Q1   0.15  0.15  0.15  0.15  0.90  0.90  0.90  0.90
Q2   0.20  0.20  0.20  0.20  0.85  0.85  0.85  0.85
Q3   0.10  0.10  0.10  0.10  0.90  0.90  0.90  0.90
Q4   0.25  0.25  0.25  0.25  0.85  0.85  0.85  0.85
Q5   0.15  0.15  0.90  0.90  0.15  0.15  0.90  0.90
Q6   0.20  0.20  0.85  0.85  0.20  0.20  0.85  0.85
Q7   0.10  0.10  0.90  0.90  0.10  0.10  0.90  0.90
Q8   0.25  0.25  0.85  0.85  0.25  0.25  0.85  0.85
Q9   0.15  0.90  0.15  0.90  0.15  0.90  0.15  0.90
Q10  0.20  0.85  0.20  0.85  0.20  0.85  0.20  0.85
Q11  0.10  0.90  0.10  0.90  0.10  0.90  0.10  0.90
Q12  0.25  0.85  0.25  0.85  0.25  0.85  0.25  0.85
Q13  0.10  0.10  0.10  0.10  0.10  0.10  0.90  0.90
Q14  0.15  0.15  0.15  0.15  0.15  0.15  0.85  0.85
Q15  0.10  0.10  0.10  0.10  0.10  0.10  0.90  0.90
Q16  0.10  0.10  0.10  0.10  0.10  0.90  0.10  0.90
Q17  0.15  0.15  0.15  0.15  0.15  0.85  0.15  0.85
Q18  0.10  0.10  0.10  0.10  0.10  0.90  0.10  0.90
Q19  0.10  0.10  0.10  0.10  0.10  0.10  0.10  0.90
Q20  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.85

============================================================
【各題成功概率 P — PyMC DINA 估計值（後驗均值）】
============================================================
       000    001    010    011    100    101    110    111
Q1   0.167  0.167  0.167  0.167  0.890  0.890  0.890  0.890
Q2   0.205  0.205  0.205  0.205  0.830  0.830  0.830  0.830
Q3   0.128  0.128  0.128  0.128  0.888  0.888  0.888  0.888
Q4   0.254  0.254  0.254  0.254  0.853  0.853  0.853  0.853
Q5   0.165  0.165  0.900  0.900  0.165  0.165  0.900  0.900
Q6   0.192  0.192  0.848  0.848  0.192  0.192  0.848  0.848
Q7   0.120  0.120  0.899  0.899  0.120  0.120  0.899  0.899
Q8   0.248  0.248  0.856  0.856  0.248  0.248  0.856  0.856
Q9   0.147  0.895  0.147  0.895  0.147  0.895  0.147  0.895
Q10  0.199  0.860  0.199  0.860  0.199  0.860  0.199  0.860
Q11  0.105  0.912  0.105  0.912  0.105  0.912  0.105  0.912
Q12  0.258  0.859  0.258  0.859  0.258  0.859  0.258  0.859
Q13  0.100  0.100  0.100  0.100  0.100  0.100  0.895  0.895
Q14  0.135  0.135  0.135  0.135  0.135  0.135  0.873  0.873
Q15  0.112  0.112  0.112  0.112  0.112  0.112  0.909  0.909
Q16  0.087  0.087  0.087  0.087  0.087  0.882  0.087  0.882
Q17  0.151  0.151  0.151  0.151  0.151  0.859  0.151  0.859
Q18  0.100  0.100  0.100  0.100  0.100  0.920  0.100  0.920
Q19  0.087  0.087  0.087  0.087  0.087  0.087  0.087  0.869
Q20  0.162  0.162  0.162  0.162  0.162  0.162  0.162  0.834

  P 矩陣 MAB  = 0.0092
  P 矩陣 RMSE = 0.0116

In [ ]:

Copied!





# ── 5c. π：真實 vs 估計（含 95% HDI）───────────────────────
pi_samples  = trace.posterior["pi"].values.reshape(-1, L)
pi_hdi_low  = np.percentile(pi_samples, 2.5,  axis=0)
pi_hdi_high = np.percentile(pi_samples, 97.5, axis=0)

df_pi = pd.DataFrame({
    "屬性組合":   profile_labels,
    "真實 π":     pi_true.round(4),
    "估計 π":     pi_est.round(4),
    "SD":         pi_samples.std(axis=0).round(4),
    "95% HDI 下": pi_hdi_low.round(4),
    "95% HDI 上": pi_hdi_high.round(4),
    "偏差":       (pi_est - pi_true).round(4),
})
print("\n" + "="*60)
print("【屬性組合比例 π — 真實 vs 估計】")
print("="*60)
print(df_pi.to_string(index=False))
# ── 5c. π：真實 vs 估計（含 95% HDI）───────────────────────
pi_samples  = trace.posterior["pi"].values.reshape(-1, L)
pi_hdi_low  = np.percentile(pi_samples, 2.5,  axis=0)
pi_hdi_high = np.percentile(pi_samples, 97.5, axis=0)

df_pi = pd.DataFrame({
    "屬性組合":   profile_labels,
    "真實 π":     pi_true.round(4),
    "估計 π":     pi_est.round(4),
    "SD":         pi_samples.std(axis=0).round(4),
    "95% HDI 下": pi_hdi_low.round(4),
    "95% HDI 上": pi_hdi_high.round(4),
    "偏差":       (pi_est - pi_true).round(4),
})
print("\n" + "="*60)
print("【屬性組合比例 π — 真實 vs 估計】")
print("="*60)
print(df_pi.to_string(index=False))

============================================================
【屬性組合比例 π — 真實 vs 估計】
============================================================
屬性組合  真實 π   估計 π     SD  95% HDI 下  95% HDI 上      偏差
 000  0.05 0.0606 0.0070     0.0475     0.0746  0.0106
 001  0.08 0.0837 0.0080     0.0692     0.1005  0.0037
 010  0.08 0.0880 0.0083     0.0726     0.1049  0.0080
 011  0.12 0.1000 0.0085     0.0839     0.1169 -0.0200
 100  0.10 0.0993 0.0081     0.0834     0.1152 -0.0007
 101  0.12 0.1184 0.0083     0.1031     0.1349 -0.0016
 110  0.15 0.1420 0.0092     0.1247     0.1608 -0.0080
 111  0.30 0.3081 0.0117     0.2849     0.3314  0.0081

In [ ]:

Copied!





# ── 5d. slip/guess 後驗分佈（95% HDI）──────────────────────
slip_samples  = trace.posterior["slip"].values.reshape(-1, J)
guess_samples = trace.posterior["guess"].values.reshape(-1, J)

slip_hdi_low   = np.percentile(slip_samples,  2.5,  axis=0)
slip_hdi_high  = np.percentile(slip_samples,  97.5, axis=0)
guess_hdi_low  = np.percentile(guess_samples, 2.5,  axis=0)
guess_hdi_high = np.percentile(guess_samples, 97.5, axis=0)

df_hdi = pd.DataFrame({
    "題目":          [f"Q{j+1}" for j in range(J)],
    "slip 真實":     slip_true.round(3),
    "slip 估計":     slip_est.round(3),
    "slip 95%下":   slip_hdi_low.round(3),
    "slip 95%上":   slip_hdi_high.round(3),
    "guess 真實":    guess_true.round(3),
    "guess 估計":    guess_est.round(3),
    "guess 95%下":  guess_hdi_low.round(3),
    "guess 95%上":  guess_hdi_high.round(3),
})
print("\n" + "="*80)
print("【slip & guess — 後驗均值與 95% HDI】")
print("="*80)
print(df_hdi.to_string(index=False))
# ── 5d. slip/guess 後驗分佈（95% HDI）──────────────────────
slip_samples  = trace.posterior["slip"].values.reshape(-1, J)
guess_samples = trace.posterior["guess"].values.reshape(-1, J)

slip_hdi_low   = np.percentile(slip_samples,  2.5,  axis=0)
slip_hdi_high  = np.percentile(slip_samples,  97.5, axis=0)
guess_hdi_low  = np.percentile(guess_samples, 2.5,  axis=0)
guess_hdi_high = np.percentile(guess_samples, 97.5, axis=0)

df_hdi = pd.DataFrame({
    "題目":          [f"Q{j+1}" for j in range(J)],
    "slip 真實":     slip_true.round(3),
    "slip 估計":     slip_est.round(3),
    "slip 95%下":   slip_hdi_low.round(3),
    "slip 95%上":   slip_hdi_high.round(3),
    "guess 真實":    guess_true.round(3),
    "guess 估計":    guess_est.round(3),
    "guess 95%下":  guess_hdi_low.round(3),
    "guess 95%上":  guess_hdi_high.round(3),
})
print("\n" + "="*80)
print("【slip & guess — 後驗均值與 95% HDI】")
print("="*80)
print(df_hdi.to_string(index=False))

================================================================================
【slip & guess — 後驗均值與 95% HDI】
================================================================================
 題目  slip 真實  slip 估計  slip 95%下  slip 95%上  guess 真實  guess 估計  guess 95%下  guess 95%上
 Q1     0.10    0.110      0.090      0.131      0.15     0.167       0.134       0.201
 Q2     0.15    0.170      0.147      0.196      0.20     0.205       0.169       0.243
 Q3     0.10    0.112      0.092      0.132      0.10     0.128       0.098       0.158
 Q4     0.15    0.147      0.124      0.170      0.25     0.254       0.214       0.295
 Q5     0.10    0.100      0.081      0.120      0.15     0.165       0.132       0.204
 Q6     0.15    0.152      0.129      0.177      0.20     0.192       0.158       0.230
 Q7     0.10    0.101      0.082      0.121      0.10     0.120       0.091       0.150
 Q8     0.15    0.144      0.121      0.166      0.25     0.248       0.210       0.287
 Q9     0.10    0.105      0.086      0.127      0.15     0.147       0.117       0.180
Q10     0.15    0.140      0.117      0.165      0.20     0.199       0.167       0.234
Q11     0.10    0.088      0.070      0.108      0.10     0.105       0.080       0.132
Q12     0.15    0.141      0.120      0.163      0.25     0.258       0.223       0.295
Q13     0.10    0.105      0.082      0.130      0.10     0.100       0.081       0.122
Q14     0.15    0.127      0.102      0.155      0.15     0.135       0.112       0.160
Q15     0.10    0.091      0.069      0.115      0.10     0.112       0.091       0.135
Q16     0.10    0.118      0.094      0.145      0.10     0.087       0.069       0.107
Q17     0.15    0.141      0.115      0.170      0.15     0.151       0.128       0.176
Q18     0.10    0.080      0.061      0.102      0.10     0.100       0.081       0.122
Q19     0.10    0.131      0.100      0.163      0.10     0.087       0.070       0.105
Q20     0.15    0.166      0.133      0.202      0.15     0.162       0.142       0.185

In [ ]:

Copied!





# ══════════════════════════════════════════════════════════════
# PART 6：各人分類結果（MAP 與 EAP）
# ══════════════════════════════════════════════════════════════

# ── 6a. 計算後驗 γ（各人對每個屬性組合的後驗概率）──────────
# 使用估計的 slip_est、guess_est、pi_est 做前向計算
log_gamma = np.zeros((N, L))
for l in range(L):
    # P(Y_i | α_l)
    log_p_given_l = np.zeros(N)
    for j in range(J):
        if eta[l, j] == 1:
            p_jl = 1 - slip_est[j]
        else:
            p_jl = guess_est[j]
        log_p_given_l += (
            Y[:, j] * np.log(p_jl + 1e-10) +
            (1 - Y[:, j]) * np.log(1 - p_jl + 1e-10)
        )
    log_gamma[:, l] = np.log(pi_est[l] + 1e-10) + log_p_given_l

log_norm = logsumexp(log_gamma, axis=1, keepdims=True)
gamma    = np.exp(log_gamma - log_norm)  # (N, L) 後驗概率

# ── 6b. MAP 分類 ─────────────────────────────────────────────
map_idx      = np.argmax(gamma, axis=1)
map_profiles = all_alpha[map_idx]

# ── 6c. EAP（各屬性的期望掌握概率）─────────────────────────
eap = gamma @ all_alpha  # (N, K)

# ── 6d. 組合分類結果 DataFrame ───────────────────────────────
df_class = pd.DataFrame({
    "真實屬性組合":  ["".join(map(str, all_alpha[i])) for i in alpha_idx],
    "MAP預測組合":   ["".join(map(str, p)) for p in map_profiles],
    "分類正確":      alpha_idx == map_idx,
    "MAP後驗概率":   gamma[np.arange(N), map_idx].round(3),
    "EAP_α1":       eap[:, 0].round(3),
    "EAP_α2":       eap[:, 1].round(3),
    "EAP_α3":       eap[:, 2].round(3),
})
for l, label in enumerate(profile_labels):
    df_class[f"P({label})"] = gamma[:, l].round(3)

print("\n" + "="*60)
print("【各人分類結果（前 15 人）】")
print("="*60)
print(df_class.head(15).to_string(index=True))

# ── 6e. 整體分類準確率 ───────────────────────────────────────
acc = df_class["分類正確"].mean()
print(f"\n  MAP 整體分類準確率 = {acc:.4f} ({acc*100:.1f}%)")

# ── 6f. 各屬性組合的分類準確率 ──────────────────────────────
print("\n  各屬性組合分類準確率：")
summary = df_class.groupby("真實屬性組合").agg(
    人數=("分類正確", "count"),
    正確數=("分類正確", "sum"),
    準確率=("分類正確", "mean"),
    平均MAP後驗=("MAP後驗概率", "mean"),
).round(3)
print(summary.to_string())

# ── 6g. 混淆矩陣 ────────────────────────────────────────────
conf = pd.crosstab(
    df_class["真實屬性組合"],
    df_class["MAP預測組合"],
    rownames=["真實"], colnames=["預測"]
)
print("\n" + "="*60)
print("【混淆矩陣（行=真實，列=預測）】")
print("="*60)
print(conf.to_string())
# ══════════════════════════════════════════════════════════════
# PART 6：各人分類結果（MAP 與 EAP）
# ══════════════════════════════════════════════════════════════

# ── 6a. 計算後驗 γ（各人對每個屬性組合的後驗概率）──────────
# 使用估計的 slip_est、guess_est、pi_est 做前向計算
log_gamma = np.zeros((N, L))
for l in range(L):
    # P(Y_i | α_l)
    log_p_given_l = np.zeros(N)
    for j in range(J):
        if eta[l, j] == 1:
            p_jl = 1 - slip_est[j]
        else:
            p_jl = guess_est[j]
        log_p_given_l += (
            Y[:, j] * np.log(p_jl + 1e-10) +
            (1 - Y[:, j]) * np.log(1 - p_jl + 1e-10)
        )
    log_gamma[:, l] = np.log(pi_est[l] + 1e-10) + log_p_given_l

log_norm = logsumexp(log_gamma, axis=1, keepdims=True)
gamma    = np.exp(log_gamma - log_norm)  # (N, L) 後驗概率

# ── 6b. MAP 分類 ─────────────────────────────────────────────
map_idx      = np.argmax(gamma, axis=1)
map_profiles = all_alpha[map_idx]

# ── 6c. EAP（各屬性的期望掌握概率）─────────────────────────
eap = gamma @ all_alpha  # (N, K)

# ── 6d. 組合分類結果 DataFrame ───────────────────────────────
df_class = pd.DataFrame({
    "真實屬性組合":  ["".join(map(str, all_alpha[i])) for i in alpha_idx],
    "MAP預測組合":   ["".join(map(str, p)) for p in map_profiles],
    "分類正確":      alpha_idx == map_idx,
    "MAP後驗概率":   gamma[np.arange(N), map_idx].round(3),
    "EAP_α1":       eap[:, 0].round(3),
    "EAP_α2":       eap[:, 1].round(3),
    "EAP_α3":       eap[:, 2].round(3),
})
for l, label in enumerate(profile_labels):
    df_class[f"P({label})"] = gamma[:, l].round(3)

print("\n" + "="*60)
print("【各人分類結果（前 15 人）】")
print("="*60)
print(df_class.head(15).to_string(index=True))

# ── 6e. 整體分類準確率 ───────────────────────────────────────
acc = df_class["分類正確"].mean()
print(f"\n  MAP 整體分類準確率 = {acc:.4f} ({acc*100:.1f}%)")

# ── 6f. 各屬性組合的分類準確率 ──────────────────────────────
print("\n  各屬性組合分類準確率：")
summary = df_class.groupby("真實屬性組合").agg(
    人數=("分類正確", "count"),
    正確數=("分類正確", "sum"),
    準確率=("分類正確", "mean"),
    平均MAP後驗=("MAP後驗概率", "mean"),
).round(3)
print(summary.to_string())

# ── 6g. 混淆矩陣 ────────────────────────────────────────────
conf = pd.crosstab(
    df_class["真實屬性組合"],
    df_class["MAP預測組合"],
    rownames=["真實"], colnames=["預測"]
)
print("\n" + "="*60)
print("【混淆矩陣（行=真實，列=預測）】")
print("="*60)
print(conf.to_string())

============================================================
【各人分類結果（前 15 人）】
============================================================
   真實屬性組合 MAP預測組合  分類正確  MAP後驗概率  EAP_α1  EAP_α2  EAP_α3  P(000)  P(001)  P(010)  P(011)  P(100)  P(101)  P(110)  P(111)
0     010     010  True    0.972   0.000   0.980   0.008   0.019   0.000   0.972   0.008   0.000   0.000     0.0     0.0
1     110     110  True    1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000     1.0     0.0
2     101     101  True    0.984   0.991   0.000   0.993   0.000   0.009   0.000   0.000   0.007   0.984     0.0     0.0
3     111     111  True    1.000   1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000     0.0     1.0
4     111     111  True    1.000   1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000     0.0     1.0
5     011     011  True    0.906   0.000   0.999   0.908   0.000   0.001   0.092   0.906   0.000   0.000     0.0     0.0
6     011     011  True    0.941   0.000   0.942   0.999   0.000   0.058   0.001   0.941   0.000   0.000     0.0     0.0
7     111     111  True    1.000   1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000     0.0     1.0
8     111     111  True    1.000   1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000     0.0     1.0
9     111     111  True    1.000   1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000     0.0     1.0
10    100     100  True    0.951   0.951   0.000   0.002   0.047   0.002   0.000   0.000   0.951   0.000     0.0     0.0
11    101     101  True    1.000   1.000   0.000   1.000   0.000   0.000   0.000   0.000   0.000   1.000     0.0     0.0
12    110     110  True    1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000     1.0     0.0
13    111     111  True    1.000   1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000     0.0     1.0
14    100     100  True    0.979   0.979   0.000   0.000   0.020   0.000   0.000   0.000   0.979   0.000     0.0     0.0

  MAP 整體分類準確率 = 0.9560 (95.6%)

  各屬性組合分類準確率：
         人數  正確數    準確率  平均MAP後驗
真實屬性組合                          
000      88   77  0.875    0.866
001     127  112  0.882    0.883
010     127  119  0.937    0.922
011     153  135  0.882    0.910
100     145  139  0.959    0.938
101     180  176  0.978    0.969
110     218  214  0.982    0.972
111     462  462  1.000    0.996

============================================================
【混淆矩陣（行=真實，列=預測）】
============================================================
預測   000  001  010  011  100  101  110  111
真實                                         
000   77    2    3    0    5    1    0    0
001    7  112    0    6    1    1    0    0
010    4    0  119    2    0    0    2    0
011    0   11    6  135    1    0    0    0
100    2    1    0    0  139    2    1    0
101    0    2    0    0    0  176    0    2
110    0    1    2    0    1    0  214    0
111    0    0    0    0    0    0    0  462