"""Case 013: same deterministic teaching data as R; NumPy only."""
import numpy as np

i = np.arange(1, 121)

def unit(a, m):
    return (((i * a) % m) - (m - 1) / 2) / ((m - 1) / 2)   # repeatable values in [-1, 1]

support = ((i * 37) % 11 - 5) / 2.5
trait = 0.8 * support + 0.6 * unit(53, 13)

def item(a, m):
    return np.clip(np.floor(3 + trait + 0.8 * unit(a, m) + 0.5), 1, 5)

# Four-item well-being scale (1-5). Items 2 and 4 are reverse-worded:
# a HIGH raw score on them means LOW well-being.
wb1, wb2, wb3, wb4 = item(17, 7), 6 - item(19, 7), item(23, 7), 6 - item(29, 7)

def slope(y, x):
    X = np.column_stack([np.ones(len(y)), x])
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    resid = y - X @ beta
    se = np.sqrt(resid @ resid / (len(y) - 2) * np.linalg.inv(X.T @ X)[1, 1])
    return beta[1], se, beta[1] / se

# Check: inter-item correlations. Negative values show items that point the other way.
print(np.round(np.corrcoef([wb1, wb2, wb3, wb4]), 2))

# The mistake: sum the raw items without reverse-scoring items 2 and 4.
wrong = slope(wb1 + wb2 + wb3 + wb4, support)

# The fix: reverse-score items 2 and 4 (minimum + maximum - response), then sum.
wb2r, wb4r = 6 - wb2, 6 - wb4
print(np.round(np.corrcoef([wb1, wb2r, wb3, wb4r]), 2))
right = slope(wb1 + wb2r + wb3 + wb4r, support)

print("Unreversed: estimate %.3f, SE %.3f, t %.3f" % wrong)
print("Reversed:   estimate %.3f, SE %.3f, t %.3f" % right)

assert len(i) == 120
assert np.allclose(wrong, [0.105, 0.079, 1.342], atol=5e-4, rtol=0)
assert np.allclose(right, [3.042, 0.147, 20.740], atol=5e-4, rtol=0)
print("All assertions passed; NumPy", np.__version__)