AI003. 缺失值恢复

# 岭回归法方程 (X'X+λR)β=X'y,R 对角阵且偏置位为 0(偏置不正则);高斯消元选列主元。
# 反推缺失特征 x_k=(y−b−Σ其余)/w_k;|w_k|<1e-10 输出 NA。
import sys


def gauss(a, b):
    n = len(b)
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(a[r][col]))
        if abs(a[pivot][col]) < 1e-12:
            raise RuntimeError("singular")
        a[col], a[pivot] = a[pivot], a[col]
        b[col], b[pivot] = b[pivot], b[col]
        div = a[col][col]
        for j in range(col, n):
            a[col][j] /= div
        b[col] /= div
        for r in range(n):
            if r == col:
                continue
            factor = a[r][col]
            if factor == 0:
                continue
            for j in range(col, n):
                a[r][j] -= factor * a[col][j]
            b[r] -= factor * b[col]
    return b


def solve() -> None:
    tokens = sys.stdin.read().split()
    if not tokens:
        return
    at = 0
    n = int(tokens[at]); at += 1
    d = int(tokens[at]); at += 1
    m = int(tokens[at]); at += 1
    lam = float(tokens[at]); at += 1
    size = d + 1
    a = [[0.0] * size for _ in range(size)]
    b = [0.0] * size
    for _ in range(n):
        z = [1.0]
        for _ in range(d):
            z.append(float(tokens[at])); at += 1
        y = float(tokens[at]); at += 1
        for r in range(size):
            b[r] += z[r] * y
            for c in range(size):
                a[r][c] += z[r] * z[c]
    for j in range(1, size):
        a[j][j] += lam
    w = gauss(a, b)
    ans = []
    for _ in range(m):
        row = []
        missing = -1
        for j in range(d):
            token = tokens[at]; at += 1
            if token == "?":
                missing = j
                row.append(0.0)
            else:
                row.append(float(token))
        target = float(tokens[at]); at += 1
        denom = w[missing + 1]
        if abs(denom) < 1e-10:
            ans.append("NA")
            continue
        known = w[0]
        for j in range(d):
            if j != missing:
                known += w[j + 1] * row[j]
        value = (target - known) / denom
        if abs(value) < 0.00005:
            value = 0.0
        ans.append(f"{value:.4f}")
    print("\n".join(ans))


if __name__ == "__main__":
    solve()