W21 실습 · 집계와 영향력 벡터

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import font_manager

# 한국어 글꼴을 후보 목록에서 있는 것으로 고른다
_cands = ['Apple SD Gothic Neo', 'AppleGothic', 'NanumGothic', 'Noto Sans CJK KR', 'Malgun Gothic']
_have = {f.name for f in font_manager.fontManager.ttflist}
for _f in _cands:
    if _f in _have:
        plt.rcParams['font.family'] = _f
        break
plt.rcParams['axes.unicode_minus'] = False

# 5부문 투입산출 몫 행렬을 직접 입력한다 — 행 = 구매자, 열 = 공급자
SECTORS = ['에너지', '소재', '제조', '물류', '서비스']
W = np.array([[0.4, 0.2, 0.2, 0.1, 0.1],
              [0.5, 0.2, 0.1, 0.1, 0.1],
              [0.2, 0.4, 0.1, 0.2, 0.1],
              [0.6, 0.0, 0.2, 0.1, 0.1],
              [0.2, 0.0, 0.2, 0.2, 0.4]])
ALPHA, N5 = 0.5, 5
B5 = np.full(N5, 1.0 / N5)
assert np.all(W >= 0) and np.allclose(W.sum(axis=1), 1.0)   # A1: 비음, 행합 1
d = W.sum(axis=0)                                           # 가중 출차수 d = W'1 (열합)


print('글꼴:', plt.rcParams['font.family'][0])
print('행합 :', W.sum(axis=1))
print('열합 d:', d, '  합 =', d.sum(), '= n')

# 영향력 벡터를 선형계 (I-(1-α)W')v = b 로 푼다 — 역행렬을 부르지 않는다
def influence(Wm, alpha, b):
    n = Wm.shape[0]
    return np.linalg.solve(np.eye(n) - (1.0 - alpha) * Wm.T, b)


v5 = influence(W, ALPHA, B5)
first = (1.0 + (1.0 - ALPHA) * d) / N5        # 급수를 k=0,1 에서 끊은 1차 근사
print(f"{'부문':<7}{'열합 d':>9}{'1차 근사':>11}{'v_i 정확':>11}{'차':>9}")
for i in range(N5):
    print(f'{SECTORS[i]:<7}{d[i]:9.2f}{first[i]:11.3f}{v5[i]:11.3f}{v5[i] - first[i]:9.3f}')
print(f'{"합":<7}{d.sum():9.2f}{first.sum():11.3f}{v5.sum():11.3f}')
print()
print(f'가장 큰 v_i : {SECTORS[int(np.argmax(v5))]} ({v5.max():.4f})')
print(f'합 1\'v      = {v5.sum():.6f}   1/α = {1.0 / ALPHA:.6f}')
print(f"고정점 확인  : max|v - b - (1-α)W'v| = {np.abs(v5 - B5 - (1 - ALPHA) * (W.T @ v5)).max():.2e}")

# 2-노름과 허핀달을 쌓고, 노이만 부분합을 반복문으로 회차마다 더한다
nv, low = np.linalg.norm(v5), 1.0 / (ALPHA * np.sqrt(N5))
herf = np.sum((v5 / v5.sum()) ** 2)
print(f'||v||_2 = {nv:.4f}   하한 1/(α√n) = {low:.4f}   상한 1/α = {1 / ALPHA:.4f}')
print(f'허핀달 = α²||v||² = {herf:.4f}   유효 부문 수 = {1 / herf:.3f}')
print()
term, acc, parts = B5.copy(), B5.copy(), [B5.copy()]
for k in range(8):
    term = (1.0 - ALPHA) * (W.T @ term)       # k 바퀴째 전파 (1-α)^k (W')^k b
    acc = acc + term
    parts.append(acc.copy())
parts = np.array(parts)
err = np.abs(parts - v5).max(axis=1)
ratio = err[1:] / err[:-1]
print(f"{'K':>3}{'부분합 합':>12}{'최대 오차':>12}{'직전 대비':>11}")
for k in range(9):
    tail = '' if k == 0 else f'{ratio[k - 1]:11.4f}'
    print(f'{k:3d}{parts[k].sum():12.4f}{err[k]:12.4f}{tail}')
print(f'오차비 평균 {ratio.mean():.4f}   1-α = {1 - ALPHA:.4f}   최대 편차 {np.abs(ratio - 0.5).max():.4f}')

# 부분합의 수렴과 영향력 벡터를 두 패널에 그린다
cols = ['#c0392b', '#3b6ea5', '#222222', '#3f8f5f', '#d98c3f']
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9.6, 3.8))
ks = np.arange(parts.shape[0])
for i in range(N5):
    ax1.plot(ks, parts[:, i], marker='o', ms=3.5, color=cols[i], label=SECTORS[i])
    ax1.axhline(v5[i], ls=':', lw=1.0, color=cols[i])
ax1.set_xlabel('전파 회차 $K$')
ax1.set_ylabel(r'영향력 부분합 $v_i^{(K)}$')
ax1.set_xticks(ks)
ax1.legend(loc='lower right', fontsize=8, ncol=2)
ax1.set_title('(a) 노이만 부분합 — 점선은 정확한 $v$', fontsize=10, loc='left')
ax2.bar(SECTORS, v5, 0.6, color=['#c0392b'] + ['#3f8f5f'] * 4)
ax2.axhline(1.0 / (ALPHA * N5), ls='--', lw=1.2, color='k')
ax2.text(4.4, 1.0 / (ALPHA * N5) + 0.015, r'균등 $1/(\alpha n)=0.4$', ha='right', fontsize=9)
ax2.text(4.4, v5.max() * 0.97, f'합 = {v5.sum():.2f} = 1/α', ha='right', va='top', fontsize=10)
ax2.set_xlabel('부문')
ax2.set_ylabel(r'Domar 가중치 $v_i$')
ax2.set_ylim(0, v5.max() * 1.15)
ax2.set_title(f'(b) 영향력 벡터 — $\\|v\\|_2$ = {nv:.3f} > {low:.3f}', fontsize=10, loc='left')
plt.tight_layout()
plt.show()

# 계수 1 족 W = 1π' 의 닫힌형과 전체 L 로 푼 값을 작은 n 에서 맞춰 본다
def rank_size_pi(n, zeta):
    S = np.arange(1, n + 1, dtype=float) ** (-1.0 / zeta)   # S_i = i^{-1/ζ}
    return S / S.sum()


def norm_v_closed(pi, alpha):
    c = (1.0 - alpha) / alpha                               # (eq-w21-13)
    return np.sqrt((1.0 + 2.0 * c) / pi.size + c ** 2 * np.sum(pi ** 2))


print(f"{'n':>5}{'π 종류':>10}{'닫힌형':>11}{'전체 L':>11}{'차':>11}")
for n in (5, 8, 20):
    for nm, pi in (('균등', np.full(n, 1.0 / n)), ('ζ=1.5', rank_size_pi(n, 1.5)),
                   ('별', np.eye(n)[0])):
        Wr = np.outer(np.ones(n), pi)                       # 모든 구매자가 같은 비율로 산다
        assert np.allclose(Wr.sum(axis=1), 1.0)
        v = influence(Wr, ALPHA, np.full(n, 1.0 / n))
        a, b = norm_v_closed(pi, ALPHA), np.linalg.norm(v)
        print(f'{n:5d}{nm:>10}{a:11.6f}{b:11.6f}{a - b:11.1e}')

# 표 2 — n을 열 자리씩 키우며 네 구조의 ||v||_2 를 쌓는다
# 계수 1 족의 영향력 벡터 v = b + ((1-α)/α)π — (eq-w21-13) 의 벡터형 (b 균등)
def v_rank1(pi, alpha):
    return np.full(pi.size, 1.0 / pi.size) + (1.0 - alpha) / alpha * pi


ns2 = [10 ** k for k in range(2, 7)]
tab = {}
print(f"{'n':>8}{'균등':>10}{'ζ=1.5':>10}{'ζ=1.1':>10}{'별':>10}{'합 1′v':>10}")
for n in ns2:
    pis = [np.full(n, 1.0 / n), rank_size_pi(n, 1.5),
           rank_size_pi(n, 1.1), np.eye(1, n, 0)[0]]
    row = [norm_v_closed(pi, ALPHA) for pi in pis]
    sums = [v_rank1(pi, ALPHA).sum() for pi in pis]          # 네 벡터의 합을 실제로 잰다
    assert abs(row[0] - 1.0 / (ALPHA * np.sqrt(n))) < 1e-12  # 균등이면 정확히 1/(α√n)
    assert max(abs(s - 1.0 / ALPHA) for s in sums) < 1e-9    # 넷 다 1'v = 1/α (eq-w21-12)
    tab[n] = row
    print(f'{n:8d}' + ''.join(f'{x:10.4f}' for x in row) + f'{max(sums):10.4f}')
print()
lg = np.log(10.0 ** 3)
for j, nm in enumerate(['균등', 'ζ=1.5', 'ζ=1.1', '별']):
    fac = tab[10 ** 3][j] / tab[10 ** 6][j]
    print(f'{nm:>6}: 10^3 → 10^6 에서 {fac:8.2f} 분의 1  '
          f'(log-log 기울기 {-np.log(fac) / lg:+.3f})')
pn = {n: np.linalg.norm(rank_size_pi(n, 1.5)) for n in ns2}
print('ζ=1.5 의 n^{1/3}||π||_2 :', ' '.join(f'{n ** (1 / 3) * pn[n]:.3f}' for n in ns2))
print(f'||π||_2 만의 기울기 {-np.log(pn[10 ** 3] / pn[10 ** 6]) / lg:+.3f}   (극한 -1/3)')

# 같은 계산을 촘촘한 n 격자에서 되풀이해 log-log 에 찍는다
# (로그 축 눈금은 직접 붙인다 — 기본 서식의 유니코드 마이너스가 한국어 글꼴에 없다)
from matplotlib.ticker import FuncFormatter, NullFormatter


def log_ticks(ax, axes='xy'):
    for nm in axes:
        a = ax.xaxis if nm == 'x' else ax.yaxis
        a.set_major_formatter(FuncFormatter(lambda v, _: f'$10^{{{round(np.log10(v))}}}$'))
        a.set_minor_formatter(NullFormatter())


ns = np.unique(np.concatenate([10.0 ** np.arange(1, 7), 3 * 10.0 ** np.arange(1, 6)])).astype(int)
curves = {nm: [] for nm in ['균등', 'ζ=1.5', 'ζ=1.1', '별']}
for n in ns:
    curves['균등'].append(norm_v_closed(np.full(n, 1.0 / n), ALPHA))
    curves['ζ=1.5'].append(norm_v_closed(rank_size_pi(n, 1.5), ALPHA))
    curves['ζ=1.1'].append(norm_v_closed(rank_size_pi(n, 1.1), ALPHA))
    curves['별'].append(norm_v_closed(np.eye(1, n, 0)[0], ALPHA))
fig, ax = plt.subplots(figsize=(6.4, 4.2))
for nm, col in zip(curves, ['#3f8f5f', '#3b6ea5', '#d98c3f', '#c0392b']):
    ax.loglog(ns, curves[nm], marker='o', ms=3.5, color=col, label=nm)
ax.loglog([1e3, 1e6], 0.045 * np.array([1.0, 10.0 ** -1.5]), ls='--', lw=1.0, color='k')
ax.text(3e4, 0.0050, '기울기 $-1/2$', fontsize=9, ha='center', va='top')
ax.loglog([1e4, 1e6], 0.06 * np.array([1.0, 1e2 ** (-1 / 3)]), ls='--', lw=1.0, color='k')
ax.text(3e5, 0.028, '기울기 $-1/3$', fontsize=9, ha='center', va='bottom')
ax.set_xlabel('부문 수 $n$')
ax.set_ylabel(r'$\|v\|_2$ = 집계 변동성 / 개별 변동성')
log_ticks(ax)
ax.legend(loc='lower left', fontsize=9)
ax.set_title('같은 합, 다른 길이 — 구조가 감쇠율을 정한다', fontsize=10, loc='left')
plt.tight_layout()
plt.show()

# iid Pareto 를 역함수법으로 뽑아 ||π||_2 의 중앙값과 10–90% 분위를 쌓는다
rng = np.random.default_rng(21)
REPS = 300
nsm = [10, 30, 100, 300, 1000, 3000, 10 ** 4, 3 * 10 ** 4, 10 ** 5]   # 본문 3절과 같은 격자
mc = {}
for zeta in (1.5, 3.0):
    med, q10, q90 = [], [], []
    for n in nsm:
        vals = np.empty(REPS)
        for r in range(REPS):
            S = rng.random(n) ** (-1.0 / zeta)          # S = U^{-1/ζ}, S_min = 1
            vals[r] = np.sqrt((S ** 2).sum()) / S.sum()
        med.append(np.median(vals))
        q10.append(np.quantile(vals, 0.10))
        q90.append(np.quantile(vals, 0.90))
    mc[zeta] = tuple(np.array(a) for a in (med, q10, q90))
for zeta in (1.5, 3.0):
    med, q10, q90 = mc[zeta]
    sc = np.array(nsm, dtype=float) ** (1.0 - 1.0 / zeta) if zeta < 2 else np.sqrt(nsm)
    lab = 'n^{1-1/ζ}' if zeta < 2 else '√n'
    print(f'--- ζ = {zeta}   척도 {lab} ---')
    print(f"{'n':>8}{'중앙값':>10}{'10%':>10}{'90%':>10}{'분위비':>10}")
    for i, n in enumerate(nsm):
        print(f'{n:8d}{sc[i] * med[i]:10.3f}{sc[i] * q10[i]:10.3f}'
              f'{sc[i] * q90[i]:10.3f}{q90[i] / q10[i]:10.2f}')
ZETA3 = 3.0
lim3 = (ZETA3 - 1.0) / np.sqrt(ZETA3 * (ZETA3 - 2.0))       # (eq-w21-18) 을 그대로 계산
print(f'ζ={ZETA3:g} 의 극한 (ζ-1)/√(ζ(ζ-2)) = {lim3:.4f}')

# 척도를 맞추기 전과 맞춘 뒤를 나란히 그린다 — 띠가 줄어드는가만 본다
# 기준선 lim3 은 앞 셀에서 (eq-w21-18) 로 계산한 값을 그대로 쓴다
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9.6, 3.8))
for zeta, col in ((1.5, '#3b6ea5'), (3.0, '#222222')):
    med, q10, q90 = mc[zeta]
    ax1.loglog(nsm, med, color=col, marker='o', ms=3.5, label=f'ζ = {zeta}')
    ax1.fill_between(nsm, q10, q90, color=col, alpha=0.2, lw=0)
    sc = np.array(nsm, dtype=float) ** (1.0 - 1.0 / zeta) if zeta < 2 else np.sqrt(nsm)
    ax2.semilogx(nsm, sc * med, color=col, marker='o', ms=3.5, label=f'ζ = {zeta}')
    ax2.fill_between(nsm, sc * q10, sc * q90, color=col, alpha=0.2, lw=0)
ax1.loglog(nsm, lim3 * np.array(nsm, dtype=float) ** -0.5, ls='--', lw=1.2, color='#3f8f5f',
           label=rf'${lim3:.3f}\,n^{{-1/2}}$')
ax1.set_xlabel('부문 수 $n$')
ax1.set_ylabel(r'$\|\pi\|_2$')
log_ticks(ax1)
ax1.legend(loc='lower left', fontsize=9)
ax1.set_title('(a) 중앙값과 10–90% 띠', fontsize=10, loc='left')
ax2.axhline(lim3, ls='--', lw=1.2, color='#3f8f5f')
ax2.text(1.2e1, lim3 + 0.22, rf'$(\zeta-1)/\sqrt{{\zeta(\zeta-2)}}={lim3:.3f}$', fontsize=9)
ax2.set_xlabel('부문 수 $n$')
ax2.set_ylabel(r'$n^{1-1/\zeta}\|\pi\|_2$   ($\zeta=3$: $\sqrt{n}\,\|\pi\|_2$)')
ax2.set_ylim(0, 2.4)
log_ticks(ax2, 'x')
ax2.legend(loc='upper right', fontsize=9)
ax2.set_title('(b) 척도를 맞춘 뒤 — 조이는가', fontsize=10, loc='left')
plt.tight_layout()
plt.show()

# 대조 표에 옮길 수치를 예측 항목 순서로 모은다
print(f'[예측 1] 최대 v_i = {SECTORS[int(np.argmax(v5))]} {v5.max():.4f} (열합 {d.max():.1f}), '
      f"합 1'v = {v5.sum():.4f} = 1/α, "
      f'||v||_2 = {nv:.4f} > {low:.4f} = 1/(α√n), '
      f'노이만 오차비 {ratio.min():.4f}–{ratio.max():.4f} (약 1-α = 0.5)')
r15 = tab[10 ** 3][1] / tab[10 ** 6][1]
r_un = tab[10 ** 3][0] / tab[10 ** 6][0]
print(f'[예측 2] ζ=1.5: {tab[10 ** 3][1]:.4f} → {tab[10 ** 6][1]:.4f} = {r15:.1f}분의 1 '
      f'(기울기 {-np.log(r15) / lg:+.3f}), '
      f'균등: {tab[10 ** 3][0]:.4f} → {tab[10 ** 6][0]:.4f} = {r_un:.1f}분의 1 '
      f'(기울기 {-np.log(r_un) / lg:+.3f})')
sel = [nsm.index(k) for k in (10 ** 2, 10 ** 3, 10 ** 4, 10 ** 5)]   # 본문 3절이 싣는 n
for zeta in (1.5, 3.0):
    med, q10, q90 = mc[zeta]
    sc = np.array(nsm, dtype=float) ** (1.0 - 1.0 / zeta) if zeta < 2 else np.sqrt(nsm)
    print(f'[예측 3] ζ={zeta} (n = 10^2, 10^3, 10^4, 10^5): '
          f'분위비 ' + ', '.join(f'{q90[i] / q10[i]:.2f}' for i in sel)
          + ', 척도 맞춘 중앙값 ' + ', '.join(f'{sc[i] * med[i]:.2f}' for i in sel)
          + (f' (극한 {lim3:.3f})' if zeta > 2 else ' (극한 상수 없음)'))